Analisis Bayesian dengan histogram sebelumnya. Mengapa menggambar simulasi dari posterior?

Ini adalah pertanyaan pemula tentang latihan dalam “Bayesian Computation with R” karya Jim Albert. Perhatikan bahwa sementara ini mungkin pekerjaan rumah, dalam kasus saya tidak, karena saya belajar metode Bayesian di R karena saya pikir saya mungkin menggunakannya dalam analisis masa depan saya.

Bagaimanapun, sementara ini adalah pertanyaan khusus, itu mungkin melibatkan pemahaman dasar metode Bayesian.

Jadi, dalam latihan 2.2, Jim Albert meminta kami untuk menganalisis percobaan lemparan sen. Lihat disini. Kita harus menggunakan histogram sebelumnya, yaitu, membagi ruang pnilai yang mungkin dalam 10 interval panjang .1dan menetapkan probabilitas sebelumnya untuk ini.

Karena saya tahu bahwa probabilitas yang sebenarnya adalah .5, dan saya pikir sangat tidak mungkin bahwa alam semesta telah mengubah hukum probabilitas atau penny yang kasar, prioritas saya adalah:

prior <- c(1,5,20,100,5000,5000,100,20,5,1)
prior <- prior/sum(prior)

di sepanjang titik tengah interval

midpt <- seq(0.05, 0.95, by=0.1)

Sejauh ini bagus. Selanjutnya, kami memutar sen 20 kali dan mencatat jumlah keberhasilan (kepala) dan kegagalan (ekor). Mudah dilakukan:

y <- rbinom(n=20,p=.5,size=1)
s <- sum(y==1)
f <- sum(y==0)

Di expeniment saya, s == 7dan f == 13. Berikutnya adalah bagian yang saya tidak mengerti:

Simulasi dari distribusi posterior dengan (1) menghitung kepadatan posterior p pada kisi nilai pada (0,1) dan (2) mengambil sampel simulasi dengan penggantian dari kisi. (Fungsi histpriordan samplesangat membantu dalam perhitungan ini). Bagaimana probabilitas interval berubah berdasarkan data Anda?

Beginilah caranya:

p <- seq(0,1, length=500)
post <- histprior(p,midpt,prior) * dbeta(p,s+1,f+1)
post <- post/sum(post)

ps <- sample(p, replace=TRUE, prob = post)

Tetapi mengapa kita melakukan itu ?

Kita dapat dengan mudah memperoleh kerapatan posterior dengan mengalikan yang sebelumnya dengan kemungkinan yang sesuai, seperti yang dilakukan pada baris dua dari blok di atas. Ini adalah sebidang distribusi posterior:

Ketika distribusi posterior dipesan, kita dapat memperoleh hasil untuk interval yang diperkenalkan dalam histogram sebelumnya dengan merangkum elemen-elemen kepadatan posterior:

post.vector <- vector()
post.vector[1] <- sum(post[p < 0.1])
post.vector[2] <- sum(post[p > 0.1 & p <= 0.2])
post.vector[3] <- sum(post[p > 0.2 & p <= 0.3])
post.vector[4] <- sum(post[p > 0.3 & p <= 0.4])
post.vector[5] <- sum(post[p > 0.4 & p <= 0.5])
post.vector[6] <- sum(post[p > 0.5 & p <= 0.6])
post.vector[7] <- sum(post[p > 0.6 & p <= 0.7])
post.vector[8] <- sum(post[p > 0.7 & p <= 0.8])
post.vector[9] <- sum(post[p > 0.8 & p <= 0.9])
post.vector[10] <- sum(post[p > 0.9 & p <= 1])

(Ahli R mungkin menemukan cara yang lebih baik untuk membuat vektor itu. Saya kira itu mungkin ada hubungannya dengan sweep?)

round(cbind(midpt,prior,post.vector),3)

      midpt prior post.vector
 [1,]  0.05 0.000       0.000
 [2,]  0.15 0.000       0.000
 [3,]  0.25 0.002       0.003
 [4,]  0.35 0.010       0.022
 [5,]  0.45 0.488       0.737
 [6,]  0.55 0.488       0.238
 [7,]  0.65 0.010       0.001
 [8,]  0.75 0.002       0.000
 [9,]  0.85 0.000       0.000
[10,]  0.95 0.000       0.000

Selain itu, kami memiliki 500 undian dari distribusi posterior yang memberi tahu kami tidak ada yang berbeda. Berikut adalah plot kepadatan undian yang disimulasikan:

Sekarang kami menggunakan data yang disimulasikan untuk mendapatkan probabilitas untuk interval kami dengan menghitung proporsi simulasi dalam interval:

sim.vector <- vector()
sim.vector[1] <- length(ps[ps < 0.1])/length(ps)
sim.vector[2] <- length(ps[ps > 0.1 & ps <= 0.2])/length(ps)
sim.vector[3] <- length(ps[ps > 0.2 & ps <= 0.3])/length(ps)
sim.vector[4] <- length(ps[ps > 0.3 & ps <= 0.4])/length(ps)
sim.vector[5] <- length(ps[ps > 0.4 & ps <= 0.5])/length(ps)
sim.vector[6] <- length(ps[ps > 0.5 & ps <= 0.6])/length(ps)
sim.vector[7] <- length(ps[ps > 0.6 & ps <= 0.7])/length(ps)
sim.vector[8] <- length(ps[ps > 0.7 & ps <= 0.8])/length(ps)
sim.vector[9] <- length(ps[ps > 0.8 & ps <= 0.9])/length(ps)
sim.vector[10] <- length(ps[ps > 0.9 & ps <= 1])/length(ps)

(Lagi: Apakah ada cara yang lebih efisien untuk melakukan ini?)

Ringkas hasil:

round(cbind(midpt,prior,post.vector,sim.vector),3)

      midpt prior post.vector sim.vector
 [1,]  0.05 0.000       0.000      0.000
 [2,]  0.15 0.000       0.000      0.000
 [3,]  0.25 0.002       0.003      0.000
 [4,]  0.35 0.010       0.022      0.026
 [5,]  0.45 0.488       0.737      0.738
 [6,]  0.55 0.488       0.238      0.236
 [7,]  0.65 0.010       0.001      0.000
 [8,]  0.75 0.002       0.000      0.000
 [9,]  0.85 0.000       0.000      0.000
[10,]  0.95 0.000       0.000      0.000

Tidak mengherankan bahwa simultion tidak menghasilkan hasil lain selain posterior, yang menjadi basisnya. Jadi, mengapa kita menggambar simulasi itu sejak awal ?

Untuk menjawab pertanyaan Anda: Bagaimana cara mengikuti yang lebih elegan?

post.vector <- vector()
post.vector[1] <- sum(post[p < 0.1])
post.vector[2] <- sum(post[p > 0.1 & p <= 0.2])
post.vector[3] <- sum(post[p > 0.2 & p <= 0.3])
post.vector[4] <- sum(post[p > 0.3 & p <= 0.4])
post.vector[5] <- sum(post[p > 0.4 & p <= 0.5])
post.vector[6] <- sum(post[p > 0.5 & p <= 0.6])
post.vector[7] <- sum(post[p > 0.6 & p <= 0.7])
post.vector[8] <- sum(post[p > 0.7 & p <= 0.8])
post.vector[9] <- sum(post[p > 0.8 & p <= 0.9])
post.vector[10] <- sum(post[p > 0.9 & p <= 1])

Cara termudah untuk melakukannya menggunakan basis R adalah:

group <- cut(p, breaks=seq(0,1,0.1), include.lowest = T)
post.vector.alt <- aggregate(post, FUN=sum, by=list(group))

Perhatikan bahwa jeda mulai dari 0 hingga 1. Ini menghasilkan:

     Group.1            x
1    [0,0.1] 3.030528e-13
2  (0.1,0.2] 1.251849e-08
3  (0.2,0.3] 6.385088e-06
4  (0.3,0.4] 6.732672e-04
5  (0.4,0.5] 2.376448e-01
6  (0.5,0.6] 7.372805e-01
7  (0.6,0.7] 2.158296e-02
8  (0.7,0.8] 2.691182e-03
9  (0.8,0.9] 1.205200e-04
10   (0.9,1] 3.345072e-07

Dan kita mempunyai:

> all.equal (post.vector.alt$x, post.vector)
[1] TRUE

Pemahaman saya adalah bahwa karena kerapatan posterior yang diperoleh dari produk dengan kerapatan sebelumnya dan kemungkinan hanya merupakan PERKAYAAN dari kerapatan posterior, maka kami tidak dapat membuat kesimpulan EXACT dari itu secara langsung.

Akibatnya, kita perlu mengambil sampel acak dari itu, dan melakukan inferensi dari sampel, seperti metode simulasi untuk posterior dari keluarga beta.