Bagaimana memprogram simulasi Monte Carlo dari paradoks kotak Bertrand?

Masalah berikut telah diposting di Halaman Facebook Mensa International:

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

Posting itu sendiri menerima 1000+ komentar tetapi saya tidak akan membahas detail tentang perdebatan di sana karena saya tahu ini adalah paradoks kotak Bertrand dan jawabannya adalah . Apa yang membuat saya tertarik di sini adalah bagaimana seseorang menjawab masalah ini menggunakan pendekatan Monte Carlo? Bagaimana algoritma untuk mengatasi masalah ini?$\frac23$

Ini usaha saya:

1. Hasilkan angka acak yang terdistribusi secara seragam antara 0 dan 1 .$N$$0$$1$
2. Biarkan acara kotak berisi 2 bola emas (kotak 1) yang dipilih kurang dari setengah.
3. Menghitung angka yang kurang dari dan memanggil hasil sebagai S .$0.5$$S$
4. Karena kepastian untuk mendapatkan bola emas jika kotak 1 dipilih dan itu hanya peluang 50% untuk mendapatkan bola emas jika kotak 2 dipilih, maka kemungkinan mendapatkan urutan GG adalah $$P(B2=G|B1=G)=\frac{S}{S+0.5(N-S)}$$

Menerapkan algoritma di atas dalam R:

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

Output dari program di atas sekitar yang hampir cocok dengan jawaban yang benar tetapi saya tidak yakin ini adalah cara yang benar. Apakah ada cara yang tepat untuk menyelesaikan masalah ini secara terprogram?$0.67$

Seperti @Henry , saya tidak merasa bahwa solusi Anda adalah Monte Carlo. Tentunya, Anda mengambil sampel dari distribusi, tetapi tidak ada hubungannya dengan meniru proses pembuatan data. Jika Anda ingin menggunakan Monte Carlo untuk meyakinkan seseorang bahwa solusi teoretisnya benar, Anda perlu menggunakan solusi yang meniru proses pembuatan data. Saya akan membayangkannya menjadi seperti di bawah ini:

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

atau menggunakan kode "vektor":

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

Perhatikan bahwa karena Anda mengkondisikan fakta bahwa bola pertama sudah ditarik dan itu adalah emas, maka kode di atas hanya dapat menggunakan dua kotak boxes <- list(c(0, 1), c(1, 1))dan kemudian sampel dari mereka x <- boxes[[sample(2, 1)]], sehingga kode akan lebih cepat karena tidak akan membuat 1/3 berjalan kosong yang kami diskon. Namun karena masalahnya sederhana dan kodenya berjalan cepat, kami dapat mensimulasikan secara eksplisit seluruh proses pembuatan data "untuk memastikan" bahwa hasilnya benar. Selain itu, langkah ini tidak diperlukan karena akan menghasilkan hasil yang sama dalam kedua kasus.

Saya merasa S/(S+0.5*(N-S))perhitungan Anda tidak benar-benar simulasi

Coba sesuatu seperti ini

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

Mengapa tidak daftar kasus saja?

Di sini: G adalah untuk "emas", S adalah untuk "perak", modal untuk ekstraksi awal:

Gg

gG

Gs

... semua kasing lainnya mengekstraksi perak awal (S) dan tidak memenuhi syarat (ekstraksi awal G).

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