¿Cómo programar una simulación de Monte Carlo de la paradoja de la caja de Bertrand?

Question

El siguiente problema ha sido publicado en la página de Facebook de Mensa International:

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

El post en sí recibió más de 1000 comentarios, pero no entraré en detalles sobre el debate allí, ya que sé que esto es La paradoja de la caja de Bertrand y la respuesta es $ \frac23 $ . Lo que me interesa aquí es cómo se responde a este problema usando un enfoque de Montecarlo. ¿Cómo es el algoritmo para resolver este problema?

Este es mi intento:

1. Generar $N$ números aleatorios uniformemente distribuidos entre $0$ y $1$ .
2. Que el evento de la caja contiene 2 bolas de oro (caja 1) seleccionadas sea menos de la mitad.
3. Cuente los números que menos de $0.5$ y llamar al resultado como $S$ .
4. Ya que es una certeza obtener una bola de oro si la caja 1 es seleccionada y es sólo el 50% de probabilidad de obtener una bola de oro si la caja 2 es seleccionada, por lo tanto la probabilidad de obtener una secuencia GG es $$P(B2=G|B1=G)= \frac {S}{S+0.5(N-S)}$$

Implementando el algoritmo de arriba en R:

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

El resultado del programa anterior es alrededor de $0.67$ que casi coinciden con la respuesta correcta, pero no estoy seguro de que esta sea la forma correcta. ¿Hay una manera apropiada de resolver este problema programáticamente?

Answer 1

Como @Henry no me parece que su solución sea un Monte Carlo. Seguramente, tomas muestras de la distribución, pero no tiene mucho que ver con imitar el proceso de generación de datos. Si quisieras usar Monte Carlo para convencer a alguien de que la solución teórica es correcta, tendrías que usar una solución que imite el proceso de generación de datos. Me imagino que sería algo como lo siguiente:

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

o utilizando código "vectorizado":

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

Fíjate que como condicionas el hecho de que la primera bola ya fue dibujada y es dorada, por lo que el código anterior podría usar simplemente dos cajas boxes <- list(c(0, 1), c(1, 1)) y luego tomar una muestra de ellos x <- boxes[[sample(2, 1)]] Así, el código sería más rápido ya que no haría los 1/3 vacíos que descontamos. Sin embargo, como el problema es sencillo y el código se ejecuta rápidamente, podríamos permitirnos simular explícitamente todo el proceso de generación de datos "para estar seguros" de que el resultado es correcto. Además, este paso no es necesario ya que daría los mismos resultados en ambos casos.

Answer 2

Siento su S/(S+0.5*(N-S)) el cálculo no es realmente una simulación

Intenta algo como esto

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

¿Por qué no enumerar simplemente los casos?

Aquí: G es para "oro", S es para "plata", capital es para la extracción inicial:

Gg

gG

Gs

... todos los demás casos implican una extracción inicial de plata (S) y no satisfacen el condicional (extracción inicial G).

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