Jak zaprogramować symulację Monte Carlo paradoksu pudełkowego Bertranda?

Następujący problem został opublikowany na stronie Mensa International na Facebooku:

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

Sam post otrzymał ponad 1000 komentarzy, ale nie będę wchodził w szczegóły na temat debaty, ponieważ wiem, że jest to paradoks skrzynki Bertranda, a odpowiedź brzmi . Zainteresowało mnie to, w jaki sposób można rozwiązać ten problem, stosując podejście Monte Carlo? Jak algorytm rozwiązuje ten problem?$\frac23$

Oto moja próba:

1. Wygeneruj równomiernie rozmieszczonych liczb losowych od do .$N$$0$$1$
2. Niech wydarzenie w pudełku zawiera 2 wybrane złote kule (pole 1), jest mniejsze niż połowa.
3. Policz liczby mniejsze niż i nazwij wynik jako$0.5$$S$ .
4. Ponieważ jest pewne, że zdobędziesz złotą piłkę, jeśli zostanie wybrane pudełko 1, a szansa na zdobycie złotej piłki jest tylko 50%, jeśli zostanie wybrane pudełko 2, stąd prawdopodobieństwo uzyskania sekwencji GG wynosi $$P(B2=G|B1=G)=\frac{S}{S+0.5(N-S)}$$

Implementowanie powyższego algorytmu w języku R:

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

Wyjście powyższego programu wynosi około co prawie odpowiada poprawnej odpowiedzi, ale nie jestem pewien, czy to jest właściwy sposób. Czy istnieje właściwy sposób programowego rozwiązania tego problemu?$0.67$

Podobnie jak @Henry , tak naprawdę nie uważam, że twoim rozwiązaniem jest Monte Carlo. Z pewnością próbujesz z dystrybucji, ale nie ma to wiele wspólnego z naśladowaniem procesu generowania danych. Jeśli chcesz użyć Monte Carlo, aby przekonać kogoś, że teoretyczne rozwiązanie jest poprawne, musisz użyć rozwiązania imitującego proces generowania danych. Wyobrażam sobie, że jest to coś takiego:

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

lub używając kodu „wektoryzowanego”:

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

Zauważ, że skoro warunek jest taki, że pierwsza kula została już narysowana i jest złota, więc powyższy kod może po prostu użyć dwóch pól, boxes <- list(c(0, 1), c(1, 1))a następnie pobrać z nich próbki x <- boxes[[sample(2, 1)]], więc kod byłby szybszy, ponieważ nie zrobiłby 1/3 puste przebiegi, które pomijamy. Ponieważ jednak problem jest prosty, a kod działa szybko, możemy sobie pozwolić na jawną symulację całego procesu generowania danych, aby „upewnić się”, że wynik jest poprawny. Poza tym krok ten nie jest potrzebny, ponieważ w obu przypadkach przyniosłby takie same wyniki.

czuję twój S/(S+0.5*(N-S)) obliczenia nie są tak naprawdę symulacją

Spróbuj czegoś takiego

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")

Dlaczego nie po prostu wymienić skrzyń?

Tutaj: G oznacza „złoto”, S oznacza „srebro”, kapitał służy do początkowej ekstrakcji:

Gg

gG

Gs

... wszystkie pozostałe przypadki wymagają wstępnej ekstrakcji srebra (S) i nie spełniają warunkowego (początkowej ekstrakcji G).

Taki P (g | G) = 2/3.