Wie programmiere ich eine Monte-Carlo-Simulation von Bertrands Box-Paradoxon?

12

Das folgende Problem wurde auf der Mensa International Facebook-Seite veröffentlicht:

Bildbeschreibung hier eingeben

Der Post selbst hat über 1000 Kommentare erhalten, aber ich werde hier nicht näher auf die Debatte eingehen, da ich weiß, dass dies Bertrands Box-Paradoxon ist und die Antwort lautet . Mich interessiert hier, wie man dieses Problem mit einem Monte-Carlo-Ansatz beantwortet. Wie löst der Algorithmus dieses Problem?23

Hier ist mein Versuch:

  1. Generiere gleichmäßig verteilte Zufallszahlen zwischen und .N01
  2. Lassen Sie das Ereignis der Box enthält 2 Goldkugeln (Box 1) weniger als die Hälfte ausgewählt.
  3. Zählen Sie die Zahlen , die weniger als und nenne das Ergebnis als .0.5S
  4. Da es eine Gewissheit ist, eine Goldkugel zu erhalten, wenn die Box 1 ausgewählt ist, und es nur eine 50% ige Chance gibt, eine Goldkugel zu erhalten, wenn die Box 2 ausgewählt ist, ist die Wahrscheinlichkeit, eine Folge GG zu erhalten,
    P(B2=G|B1=G)=SS+0.5(NS)

Implementierung des obigen Algorithmus in R:

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

Die Ausgabe des obigen Programms liegt bei was fast der richtigen Antwort entspricht, aber ich bin mir nicht sicher, ob dies der richtige Weg ist. Gibt es einen geeigneten Weg, um dieses Problem programmgesteuert zu lösen?0.67

Anastasiya-Romanova 秀
quelle

Antworten:

14

Wie bei @Henry glaube ich nicht wirklich, dass Ihre Lösung ein Monte Carlo ist. Sicher, Sie nehmen ein Beispiel aus der Distribution, aber es hat nicht viel damit zu tun, den Datenerzeugungsprozess zu imitieren. Wenn Sie Monte Carlo verwenden möchten, um jemanden davon zu überzeugen, dass die theoretische Lösung korrekt ist, müssen Sie eine Lösung verwenden, die den Datenerzeugungsprozess imitiert. Ich würde es mir wie folgt vorstellen:

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

oder mit "vektorisiertem" 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

Beachten Sie, dass, da Sie davon ausgehen, dass die erste Kugel bereits gezogen wurde und sie golden ist, der obige Code einfach zwei Kästchen verwenden boxes <- list(c(0, 1), c(1, 1))und sie dann abtasten könnte x <- boxes[[sample(2, 1)]], der Code also schneller wäre, da er das 1/3 nicht ergeben würde Leerfahrten, die wir rabattieren. Da das Problem jedoch einfach ist und der Code schnell ausgeführt wird, können wir es uns leisten, den gesamten Datenerzeugungsprozess explizit zu simulieren, "um sicherzustellen", dass das Ergebnis korrekt ist. Außerdem ist dieser Schritt nicht erforderlich, da er in beiden Fällen zu den gleichen Ergebnissen führen würde.

Tim
quelle
Bedeutet das, x <- boxes[[sample(3, 1)]]dass Sie eine Schachtel aus 3 Schachteln nehmen? Wenn ja, warum ist es notwendig, seit wir wissen, dass Sie bereits eine goldene Kugel ausgewählt haben?
Anastasiya-Romanova 秀
7
@ Anastasiya-Romanova 秀 Sie könnten stattdessen aus zwei Kästchen abtasten boxes <- list(c(0, 1), c(1, 1))und dann x <- boxes[[sample(2, 1)]], aber da dies fast dieselbe Rechenzeit ist, warum nicht den zusätzlichen Schritt verwenden, der genau dem Abtastvorgang ähnelt? Es ändert nichts am Ergebnis, macht aber die Simulation explizit.
Tim
Alles klar Tim, danke für deine Antwort. Geben Sie mir zunächst Zeit, um Ihre Antwort zu verstehen, da ich in R ziemlich neu bin. Im Moment +1 für Sie und @ Henry.
Anastasiya-Romanova 秀
1
@ Anastasiya-Romanova 秀 ja, genau. Der Code tastet eine Schachtel ab, tastet dann einen Ball aus der Schachtel ab, ist er golden (= 1), prüft er, ob der andere Ball aus der Schachtel ebenfalls golden ist (1 + 1 = 2), wenn ja, zählt er und am Ende wird die Summe durch die Gesamtzahl geteilt (dh verwendet mean).
Tim
1
@ Anastasiya-Romanova 秀return(NA)gibt den fehlenden Wert zurück und wird dann mean(, na.rm = TRUE)verwendet, wobei das na.rm = TRUEArgument die Funktion anweist, die fehlenden Werte zu ignorieren. In anderen Programmiersprachen kann dies anders erfolgen, z . B. mit continueoder passSchlüsselwörtern.
Tim
4

Ich glaube, Ihre S/(S+0.5*(N-S))Berechnung ist nicht wirklich eine Simulation

Versuchen Sie so etwas

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")
Henry
quelle
-2

Warum nicht einfach die Fälle auflisten?

Hier: G steht für "Gold", S steht für "Silber", Kapital für die anfängliche Gewinnung:

Gg

gG

Gs

... in allen anderen Fällen wird eine anfängliche Silber (S) -Extraktion vorausgesetzt und die Bedingung (anfängliche Extraktion G) nicht erfüllt.

So ist P (g | G) = 2/3.

ghuezt
quelle
7
Die Frage fragt nach der Monte-Carlo-Lösung.
Tim
Nun, ALLE Möglichkeiten aufzulisten ist ein extremer Fall von Monte Carlo.
Ghuezt
4
Nein, Monte Carlo handelt von Simulation / Randomisierung.
Tim
Tim, mach deine Rechnung richtig. Mit unendlich vielen Ziehungen erhalten Sie genau eine Liste aller Fälle mit gleichen Wahrscheinlichkeiten. Ich traurig "Extremfall" und meinte Grenze.
Ghuezt
1
Sicher, aber alle Fälle aufzuzählen ist nicht Monte Carlo. Quadrat ist ein Rechteck, aber Rechteck ist kein Quadrat.
Tim