Search code examples
rmontecarlo

R: Unfair Coin Toss with Monte Carlo Method. Goal is Probability of N Exact Heads


I have a task to use the Monte Carlo method to evaluate an unfair coin flip and determine the probability of obtaining n heads out of n flips within n simulations. The guide my professor provided (at bottom) isn't making sense for how to adjust it for an unfair coin. I've included what I started coding, but it's nowhere near complete. How can I adjust Monte Carlo to work for an unfair coin?

Problem:

Assume you have a coin that is not fair, where the probability of having a heads is p. We are interested in the probability of seeing exactly nheads heads out of a total of nflips flips. Write a function that takes these arguments to evaluate this probability over nsim simulations (a total of four arguments), and test it out on a few different values. Hint: use the prob argument in the sample function.

nsim <- 10000
nheads <-
nflips <- 
p <- NULL
unfairfunc <- function(x){
  for(nheads in 1:nsim)
  nheads <-
  nflips <-
    }

Below is the base code my professor provided for Monte Carlo probabilities for a fair coin toss.

 flip_function <- function(n) {
  flips <- sample(c("heads", "tails"), n, replace=TRUE) 
  percent_heads <- length(which(flips== "heads")) / n
  return(percent_heads)
}

Solution

  • sample function has prob argument, which can be used to assign a probability weight value for each possible outcome. If none is assigned to prob, the sample function will take simple random sampling in which each possible outcome has the same probability weight, so each will have the same probability to be selected, which is a fair condition.

    For example:

    # Case 1:
    sample(c("a", "b"), 1) 
    
    #Case 2:
    sample(c("a", "b"), 1, prob = c(0.2, 0.8)) 
    

    In Case #1, the probability of "a" being sampled equals the probability of "b" being sampled, which is 1/2 = 50%.

    In Case #2, the probability of "a" being sampled is 20%, but the probability of "b" being sampled is 80%.

    If you repeat the sample function N times, it can be expected that "a" will be selected about 50% N times in Case #1, and about 20% N times in Case #2.

    Try many times:

    Case #1 
    s <- numeric(10)
    for(k in 1:10) s[k] <- sample(c("a", "b"), 1)
    table(s)
    
    Case #2 
    s <- numeric(10)
    for(k in 1:10) s[k] <- sample(c("a", "b"), 1, prob = c(0.2, 0.8))
    table(s)
    

    Note that :

    • if prob is to be set, a weight value must be assigned to each possible outcome, and
    • the sum of the weights assigned to prob must equal 1.