How to program a Monte Carlo simulation of Bertrand's box paradox?

Question

The following problem has been posted on Mensa International Facebook Page:

$\quad\quad\quad\quad\quad\quad\quad\quad$

The post itself received 1000+ comments but I won't go into details about the debate there since I know this is Bertrand's box paradox and the answer is $\frac23$. What make me interested here is how does one answer this problem using a Monte Carlo approach? How is the algorithm to solve this problem?

Here is my attempt:

1. Generate $N$ uniformly distributed random numbers between $0$ and $1$.
2. Let the event of box contains 2 gold balls (box 1) selected be less than half.
3. Count the numbers that less than $0.5$ and call the result as $S$.
4. Since it's a certainty to get a gold ball if the box 1 is selected and it's only 50% chance of getting a gold ball if the box 2 is selected, hence the probability of getting a sequence GG is $$P(B2=G|B1=G)=\frac{S}{S+0.5(N-S)}$$

Implementing the algorithm above in R:

N <- 10000
S <- sum(runif(N)<0.5)
S/(S+0.5*(N-S))

The output of program above is around $0.67$ which almost match with the correct answer but I'm not sure this is the correct way. Is there a proper way to solve this problem programmatically?

Answer 1

Like @Henry, I don't really feel that your solution is a Monte Carlo. Surely, you sample from the distribution, but it doesn't have much to do with imitating the data generating process. If you wanted to use Monte Carlo to convince someone that the theoretical solution is correct, you would need to use solution that imitates the data generating process. I would imagine it to be something like below:

boxes <- list(
  c(0, 0),
  c(0, 1),
  c(1, 1)
)

count_successes = 0
count_valid_samples = 0

for (i in 1:5000) {
  sampled_box <- unlist(sample(boxes, 1)) # sample box
  sampled_balls <- sample(sampled_box)    # shuffle balls in the box

  if (sampled_balls[1] == 1) {            # if first ball is golden
    if (sampled_balls[2] == 1) {          # if second ball is golden
      count_successes = count_successes + 1
    }
    count_valid_samples = count_valid_samples + 1
  }
}
count_successes / count_valid_samples

or using "vectorized" code:

mean(replicate(5000, {       # repeat 5000 times, next calculate empirical probability
  x <- boxes[[sample(3, 1)]] # pick a box
  if (x[sample(2, 1)] == 1)  # pick a ball, check if it is golden
    return(sum(x) == 2)      # check if you have two golden balls in the box
  else
    return(NA)               # ignore if sampled ball is silver
  }), na.rm = TRUE)          # not count if silver

Notice that since you condition of the fact that first ball was already drawn and it is golden, so the above code could simply use two boxes boxes <- list(c(0, 1), c(1, 1)) and then sample from them x <- boxes[[sample(2, 1)]], so the code would be faster since it wouldn't make the 1/3 empty runs that we discount. However since the problem is simple and the code runs fast, we could afford to explicitly simulate the whole data generating process "to be sure" that the result is correct. Besides that, this step is not needed since it would yield same results in both cases.

Answer 2

I feel your S/(S+0.5*(N-S)) calculation is not really simulation

Try something like this

N <- 10^6
ballsinboxes <- c("G","G", "G","S", "S","S")
selectedbox <- sample(c(1,2,3), N, replace=TRUE)
selectedball <- sample(c(1,2), N, replace=TRUE)
selectedcolour <- ballsinboxes[(selectedbox-1)*2 + selectedball ]
othercolour <- ballsinboxes[(selectedbox-1)*2 + 3 - selectedball ]
sum(selectedcolour == "G" & othercolour == "G") / sum(selectedcolour == "G")

Answer 3

Why not simply list the cases?

Here: G is for "gold", S is for "silver", capital is for the initial extraction:

Gg

gG

Gs

... all other cases invove an initial silver (S) extraction and don't satisfy the conditional (initial extraction G).

Such, P(g|G) = 2/3.