Comment programmer une simulation Monte Carlo du paradoxe de la boîte de Bertrand?

12

Le problème suivant a été publié sur la page Facebook de Mensa International:

entrez la description de l'image ici

Le message lui-même a reçu plus de 1000 commentaires, mais je n'entrerai pas dans les détails du débat là-bas, car je sais que c'est le paradoxe de la boîte de Bertrand et la réponse est . Ce qui m'intéresse ici, c'est comment répondre à ce problème en utilisant une approche Monte Carlo? Comment est l'algorithme pour résoudre ce problème?23

Voici ma tentative:

  1. Générez nombres aléatoires uniformément répartis entre et .0 1N01
  2. Que l'événement de la boîte contient 2 boules d'or (boîte 1) sélectionnées inférieures à la moitié.
  3. Compter le nombre que moins de et de communiquer avec le résultat que .S0.5S
  4. Puisqu'il est certain d'obtenir une boule d'or si la case 1 est sélectionnée et que c'est seulement 50% de chances d'obtenir une boule d'or si la case 2 est sélectionnée, d'où la probabilité d'obtenir une séquence GG est
    P(B2=G|B1=G)=SS+0.5(NS)

Implémentation de l'algorithme ci-dessus dans R:

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

La sortie du programme ci-dessus est d'environ ce qui correspond presque à la bonne réponse, mais je ne suis pas sûr que ce soit la bonne façon. Existe-t-il un moyen approprié de résoudre ce problème par programme?0.67

Anastasiya-Romanova 秀
la source

Réponses:

14

Comme @Henry , je ne pense pas vraiment que votre solution soit un Monte Carlo. Certes, vous échantillonnez à partir de la distribution, mais cela n'a pas grand-chose à voir avec l'imitation du processus de génération de données. Si vous souhaitez utiliser Monte Carlo pour convaincre quelqu'un que la solution théorique est correcte, vous devez utiliser une solution qui imite le processus de génération de données. J'imagine que c'est quelque chose comme ci-dessous:

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

ou en utilisant du code "vectorisé":

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

Notez que puisque vous conditionnez le fait que la première balle a déjà été tirée et qu'elle est dorée, le code ci-dessus pourrait simplement utiliser deux cases boxes <- list(c(0, 1), c(1, 1)), puis en échantillonner x <- boxes[[sample(2, 1)]], de sorte que le code serait plus rapide car il ne ferait pas le 1/3 pistes vides que nous actualisons. Cependant, comme le problème est simple et que le code s'exécute rapidement, nous pourrions nous permettre de simuler explicitement l'ensemble du processus de génération de données "pour être sûr" que le résultat est correct. En plus de cela, cette étape n'est pas nécessaire car elle donnerait les mêmes résultats dans les deux cas.

Tim
la source
Cela x <- boxes[[sample(3, 1)]]signifie- t - il que vous prenez une boîte de 3 boîtes? Si oui, pourquoi est-ce nécessaire puisque nous savons que vous avez déjà choisi une boule d'or?
Anastasiya-Romanova
7
@ Anastasiya-Romanova 秀 vous pourriez à la place échantillonner à partir de deux boîtes boxes <- list(c(0, 1), c(1, 1))et ensuite x <- boxes[[sample(2, 1)]], mais comme il s'agit presque du même temps de calcul, pourquoi ne pas utiliser l'étape supplémentaire qui ressemble exactement au processus d'échantillonnage? Cela ne change rien au résultat, mais rend la simulation explicite.
Tim
Très bien Tim, merci pour votre réponse. Donnez-moi le temps de comprendre votre réponse en premier car je suis assez nouveau dans R. Pour l'instant, +1 pour vous et @Henry.
Anastasiya-Romanova
1
@ Anastasiya-Romanova 秀 oui, exactement. Le code échantillonne une boîte, puis échantillonne une balle de la boîte, si elle est dorée (= 1) puis il vérifie si l'autre balle de la boîte est également dorée (1 + 1 = 2), si oui alors il la compte et à la fin, il divise la somme par le nombre total (c'est-à-dire les utilisations mean).
Tim
1
@ Anastasiya-Romanova 秀return(NA)renvoie la valeur manquante, puis mean(, na.rm = TRUE)est utilisée, où l' na.rm = TRUEargument indique à la fonction d'ignorer les valeurs manquantes. Dans d'autres langages de programmation, cela pourrait être fait différemment, par exemple en utilisant continueou des passmots clés.
Tim
4

Je sens que ton S/(S+0.5*(N-S))calcul n'est pas vraiment de la simulation

Essayez quelque chose comme ça

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")
Henri
la source
-2

Pourquoi ne pas simplement lister les cas?

Ici: G est pour "or", S est pour "argent", le capital est pour l'extraction initiale:

Gg

gG

Gs

... tous les autres cas impliquent une extraction initiale d'argent (S) et ne satisfont pas au conditionnel (extraction initiale G).

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

ghuezt
la source
7
La question porte sur la solution de Monte Carlo.
Tim
Eh bien, lister TOUTES les possibilités est un cas extrême de Monte Carlo.
ghuezt
4
Non, ce n'est pas le cas, Monte Carlo concerne la simulation / randomisation.
Tim
Tim, fais bien tes calculs. Avec une infinité de tirages, vous obtenez exactement une liste de tous les cas avec des probabilités égales. J'ai triste "cas extrême" et voulait dire limite.
ghuezt
1
Bien sûr, mais répertorier tous les cas n'est pas Monte Carlo. Le carré est un rectangle, mais le rectangle n'est pas un carré.
Tim