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

12

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

wprowadź opis zdjęcia tutaj

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?2)3)

Oto moja próba:

  1. Wygeneruj równomiernie rozmieszczonych liczb losowych od do .N.101
  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 jakoS.0,5S. .
  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)=sol|b1=sol)=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

Anastasiya-Romanova 秀
źródło

Odpowiedzi:

14

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.

Tim
źródło
Czy to x <- boxes[[sample(3, 1)]]znaczy, że bierzesz pudełko z 3 pudełek? Jeśli tak, dlaczego jest to konieczne, skoro wiemy, że już wybrałeś złotą piłkę?
Anastasiya-Romanova 秀
7
@ Anastasiya-Romanowa秀mógłbyś zamiast próbki z dwóch pól boxes <- list(c(0, 1), c(1, 1)), a następnie x <- boxes[[sample(2, 1)]], ale ponieważ jest to prawie taki sam czas obliczeń, dlaczego nie wykorzystać dodatkowy krok, który dokładnie podobny proces próbkowania? Nie zmienia to niczego w wyniku, ale wyraźnie określa symulację.
Tim
Dobra Tim, dzięki za odpowiedź. Daj mi czas, aby najpierw zrozumieć twoją odpowiedź, ponieważ jestem dość nowy w R. Na razie +1 dla ciebie i @Henry.
Anastasiya-Romanova 秀
1
@ Anastasiya-Romanova 秀 tak, dokładnie. Kod pobiera próbkę z pudełka, a następnie pobiera próbkę z pudełka, jeśli jest złota (= 1), następnie sprawdza, czy druga piłka z pudełka jest również złota (1 + 1 = 2), jeśli tak, to liczy i na koniec dzieli sumę przez całkowitą liczbę (tj. wykorzystanie mean).
Tim
1
@ Anastasiya-Romanova 秀return(NA)zwraca brakującą wartość, a następnie mean(, na.rm = TRUE)jest używana, gdzie na.rm = TRUEargument mówi funkcji, aby zignorowała brakujące wartości. W innych językach programowania można to zrobić inaczej, np. Za pomocą continuelub passsłów kluczowych.
Tim
4

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")
Henz
źródło
-2

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.

Ghuezt
źródło
7
Pytanie dotyczy rozwiązania Monte Carlo.
Tim
Cóż, zestawienie WSZYSTKICH możliwości jest skrajnym przypadkiem Monte Carlo.
ghuezt
4
Nie, nie jest, Monte Carlo dotyczy symulacji / randomizacji.
Tim
Tim, popraw swoją matematykę. Przy nieskończenie wielu losowaniach otrzymujesz dokładnie listę wszystkich przypadków z jednakowymi prawdopodobieństwami. Smutno mi „ekstremalny przypadek” i miałem na myśli limit.
ghuezt
1
Jasne, ale wymienianie wszystkich przypadków nie jest Monte Carlo. Kwadrat jest prostokątem, ale prostokąt nie jest kwadratem.
Tim