Mensimulasikan variabel Bernoulli dengan probabilitas menggunakan koin bias

Dapatkah seseorang memberi tahu saya cara mensimulasikan , di mana , menggunakan lemparan koin (sebanyak yang Anda perlukan) dengan ?$\mathrm{Bernoulli}\left({a\over b}\right)$$a,b\in \mathbb{N}$$P(H)=p$

Saya sedang berpikir untuk menggunakan sampel penolakan, tetapi tidak bisa mengatasinya.

Karena ada banyak solusi yang tidak terhitung jumlahnya, mari kita cari yang efisien .

Ide di balik ini satu dimulai dengan cara standar untuk menerapkan variabel Bernoulli: membandingkan variabel acak seragam untuk parameter . Saat , kembalikan ; jika tidak, kembalikan .$U$$a/b$$U \lt a/b$$1$$0$

Kita dapat menggunakan -coin sebagai penghasil angka acak yang seragam$p$ . Untuk menghasilkan angka secara seragam dalam interval apa pun , balikkan koin. Ketika head, secara rekursif menghasilkan nilai seragam di bagian pertama dari interval; ketika itu buntut, secara menghasilkan secara rekursif dari bagian terakhir dari interval. Pada titik tertentu, interval target akan menjadi sangat kecil sehingga tidak masalah bagaimana Anda memilih nomor dari itu: itulah cara rekursi dimulai. Jelas prosedur ini menghasilkan variasi yang seragam (hingga presisi yang diinginkan), seperti yang mudah dibuktikan dengan induksi.$U$$[x, y)$$X$$p$$X$$1-p$

Ide ini tidak efisien, tetapi mengarah pada metode yang efisien. Karena pada setiap tahap Anda akan menggambar angka dari beberapa interval yang diberikan , mengapa tidak pertama-tama memeriksa apakah Anda perlu menggambar sama sekali? Jika nilai target Anda berada di luar interval ini, Anda sudah tahu hasil perbandingan antara nilai acak dan target. Dengan demikian, algoritma ini cenderung berakhir dengan cepat. (Ini dapat ditafsirkan sebagai prosedur pengambilan sampel penolakan yang diminta dalam pertanyaan.)$[x,y)$

Kami dapat mengoptimalkan algoritme ini lebih lanjut. Pada tahap apa pun, kita sebenarnya memiliki dua koin yang dapat kita gunakan: dengan memberi label ulang koin kita, kita dapat membuatnya menjadi satu yang berkepala . Oleh karena itu, sebagai perhitungan awal kami dapat secara rekursif memilih mana yang menandai ulang mengarah ke jumlah flips yang diharapkan lebih rendah yang diperlukan untuk penghentian. (Perhitungan ini bisa menjadi langkah mahal.)$1-p$

Misalnya, tidak efisien untuk menggunakan koin dengan untuk mengemulasi variabel Bernoulli secara langsung: dibutuhkan rata-rata hampir sepuluh flips. Tetapi jika kita menggunakan koin, maka hanya dalam dua flips kita pasti akan selesai dan jumlah flips yang diharapkan hanya .$p=0.9$$(0.01)$$p=1-0.0=0.1$$1.2$

Berikut detailnya.

Partisi setiap interval setengah terbuka ke dalam interval$I = [x, y)$

Ini mendefinisikan dua transformasi dan yang beroperasi pada interval setengah terbuka.$s(*,H)$$s(*,T)$

Sebagai masalah terminologi, jika adalah kumpulan bilangan real biarkan ungkapan$I$

berarti bahwa adalah batas bawah untuk : untuk semua . Demikian pula, berarti adalah batas atas untuk .$t$$I$$t \lt x$$x \in I$$t \gt I$$t$$I$

Tulis . (Sebenarnya, tidak ada bedanya jika nyata dan bukan rasional; kita hanya memerlukan )$a/b = t$$t$$0 \le t \le 1$

Berikut adalah algoritma untuk menghasilkan varian dengan parameter Bernoulli yang diinginkan:$Z$

1. Set dan I n = I 0 = [ 0 , 1 ) .$n=0$$I_n = I_0 = [0,1)$

2. While {Aduk koin untuk menghasilkan X n + 1 . Set I n + 1 = S ( I n , X n + 1 ) . Penambahan n .}$(t\in I_{n})$$X_{n+1}$$I_{n+1} = S(I_n, X_{n+1}).$$n$

3. Jika maka atur Z = 1 . Jika tidak, atur Z = 0 .$t \gt I_{n+1}$$Z=1$$Z=0$

Penerapan

Sebagai ilustrasi, berikut ini adalah Rimplementasi dari aloritma sebagai fungsinya draw. Argumennya adalah nilai target dan interval [ x , y ) , awalnya [ 0 , 1 ) . Ini menggunakan fungsi bantu menerapkan . Meskipun tidak perlu, itu juga melacak jumlah lemparan koin. Ini mengembalikan variabel acak, jumlah lemparan, dan interval terakhir yang diperiksa.$t$$[x,y)$$[0,1)$s$s$

s <- function(x, ab, p) {
  d <- diff(ab) * p
  if (x == 1) c(ab[1], ab[1] + d) else c(ab[1] + d, ab[2])
}
draw <- function(target, p) {
  between <- function(z, ab) prod(z - ab) <= 0
  ab <- c(0,1)
  n <- 0
  while(between(target, ab)) {
    n <- n+1; ab <- s(runif(1) < p, ab, p)
  }
  return(c(target > ab[2], n, ab))
}

Sebagai contoh penggunaan dan uji keakuratannya, ambil case dan . Mari menggambar nilai menggunakan algoritma, melaporkan rata-rata (dan kesalahan standarnya), dan menunjukkan jumlah rata-rata flips yang digunakan.$t=1/100$$p=0.9$$10,000$

target <- 0.01
p <- 0.9
set.seed(17)
sim <- replicate(1e4, draw(target, p))

(m <- mean(sim[1, ]))                           # The mean
(m - target) / (sd(sim[1, ]) / sqrt(ncol(sim))) # A Z-score to compare to `target`
mean(sim[2, ])                                  # Average number of flips

Dalam simulasi ini flips adalah kepala. Meskipun lebih rendah dari target , skor-Z tidak signifikan: penyimpangan ini dapat dikaitkan dengan kebetulan. Jumlah rata-rata flips adalah - sedikit kurang dari sepuluh. Jika kita menggunakan koin , nilai rata-rata adalah masih tidak jauh berbeda dari target, tetapi hanya flips yang diperlukan rata-rata.$0.0095$$0.01$$-0.5154$$9.886$$1-p$$0.0094$$1.177$

Ini solusinya (agak berantakan, tapi ini tikaman pertamaku). Anda benar-benar dapat mengabaikan dan WLOG menganggap . Mengapa? Ada algoritma pintar untuk menghasilkan flip koin yang tidak bias dari dua membalik koin bias. Jadi kita dapat mengasumsikan .$P(H) = p$$P(H)=1/2$$P(H)=1/2$

Untuk menghasilkan , saya dapat memikirkan dua solusi (yang pertama bukan milik saya, tetapi yang kedua adalah generalisasi):$\text{Bernoulli}(\frac{a}{b})$

Solusi 1

Balikkan koin yang tidak bias kali. Jika kepala tidak hadir, mulai dari awal. Jika kepala yang hadir, kembali apakah koin pertama adalah kepala atau tidak (karena )$b$$a$$a$$P(\text{first coin is heads | $a$ heads in $b$ coins}) = \frac{a}{b}$

Solusi 2

Ini dapat diperluas ke nilai . Tulis dalam bentuk biner. Misalnya,$\text{Bernoulli}(p)$$p$$0.1 = 0.0001100110011001100110011... \text{base 2}$

Kami akan membuat nomor biner baru menggunakan koin membalik. Mulai dengan , dan tambahkan digit tergantung pada apakah kepala (1) atau ekor (0) muncul. Pada setiap flip, bandingkan angka biner baru Anda dengan representasi biner hingga digit yang sama . Akhirnya keduanya akan berbeda, dan kembali jika lebih besar dari nomor biner Anda.$0.$$p$ $bin(p)$

Dengan Python:

def simulate(p):
    binary_p = float_to_binary(p)
    binary_string = '0.'
    index = 3
    while True:
        binary_string += '0' if random.random() < 0.5 else '1'
        if binary_string != binary_p[:index]:
            return binary_string < binary_p[:index]
        index += 1

Beberapa bukti:

np.mean([simulate(0.4) for i in range(10000)])

sekitar 0,4 (namun tidak cepat)

Saya melihat solusi sederhana, tetapi tidak diragukan lagi ada banyak cara untuk melakukannya, beberapa mungkin lebih sederhana dari ini. Pendekatan ini dapat dipecah menjadi dua langkah:

1. Menghasilkan dari dua peristiwa dengan probabilitas yang sama diberikan prosedur melempar koin yang tidak adil (kombinasi koin tertentu dan metode yang digunakan untuk menghasilkan kepala dengan probabilitas ). Kita dapat menyebut dua peristiwa yang sama-sama berpeluang ini , dan . [Ada pendekatan sederhana untuk ini yang mengharuskan mengambil pasangan lemparan dan untuk menghasilkan dua hasil yang kemungkinannya sama, dengan semua hasil lainnya mengarah pada menghasilkan pasangan baru gulungan untuk mencoba lagi.]$p$$H^*$$T^*$$H^*=(H,T)$$T^*=(T,H)$

2. Sekarang Anda menghasilkan jalan acak dengan dua negara menyerap menggunakan koin adil yang disimulasikan. Dengan memilih jarak kondisi serapan dari titik asal (satu di atas dan satu di bawahnya), Anda dapat mengatur peluang penyerapan dengan mengatakan kondisi serapan atas menjadi rasio bilangan bulat yang diinginkan. Khususnya, jika Anda menempatkan penghalang penyerap atas di dan yang lebih rendah di (dan memulai proses dari asal), dan menjalankan jalan acak hingga penyerapan, probabilitas penyerapan di penghalang atas adalah .$a$$-(b-a)$$\frac{a}{a+(b-a)} = \frac{a}{b}$

   (Ada beberapa perhitungan yang harus dilakukan di sini untuk menunjukkannya, tetapi Anda bisa mengeluarkan probabilitas dengan cukup mudah dengan bekerja dengan relasi perulangan ... atau Anda dapat melakukannya dengan menjumlahkan seri tak terbatas ... atau ada cara lain.)