A Tricky Dice Problem

Step 1.

One replication of the experiment consists of rolling a die a very large number of times. Create a function roll(nroll) that rolls a fair die nroll times and returns TRUE if there are more low than high numbers AND ALSO there are more odd than even numbers, and returns FALSE otherwise.

roll <- function(nroll) {
  rolls <- sample(1:6, nroll, replace = T)
  low_count <- sum(rolls <= 3)
  odd_count <- sum(rolls %% 2 == 1)
  if ((low_count > nroll/2) && (odd_count > nroll/2)) {
    return(TRUE)
  }
  return(FALSE)
}

We will consider 10,000 large enough and use nroll=10000.

Step 2.

The roll function just constructed gives us a single replication of the experiment. To find the probability, we must replicate the experiment a large number of times and find the proportion of TRUE replications.

Replicate the experiment 500 times, with 10,000 rolls each time. The result should be a logical vector of length 500. Assign this result to a variable called replicates.

replicates <- replicate(500, roll(10000))

Step 3.

What is the estimated answer to the dice problem?

Our estimate is the proportion of TRUE replications

mean(replicates)
## [1] 0.31

Step 4.

We can visualize our approximation by plotting the running mean; we should expect that our approximation gets better with more replicates.

First, create a function running_mean(m) that gives the average of the first m elements of replicates. Check that running_mean(500) gives your answer from Step 3.

running_mean <- function(m) {
  mean(replicates[1:m])
}
running_mean(500)
## [1] 0.31

Next, find the running mean with the following line of code:

running <- sapply(1:500, running_mean)

running now contains the running mean over all 500 replicates.

Finally, make the following plot:

plot(1:500, running, type="l",
     xlab = "replication", ylab = "estimate", main = "Approximation")
p <- 1/4 + asin(1/3) / (2 * pi)
abline(p, 0, col="red")

The red horizontal line is the true value given in \((*)\). You can find the true value with asin, which is the \(\arcsin\) function.

Your plot will of course be different due to randomness. Remember to set a seed so that your results are replicable.