Saya memiliki model untuk mencapai perkiraan Bayesian ukuran populasi dan probabilitas deteksi dalam distribusi binomial semata-mata berdasarkan pada jumlah obyek yang diamati yang diamati : untuk . Untuk kesederhanaan, kita asumsikan bahwa N ditetapkan pada nilai yang sama dan tidak diketahui untuk setiap y_i . Dalam contoh ini, y = 53,57,66,67,73 .
Model ini, ketika diperkirakan dalam rstan
, menyimpang dari hasil yang diperoleh dari pendekatan grid posterior. Saya mencoba menjelaskan mengapa. (Pembaca yang tertarik mungkin menemukan bahwa pertanyaan ini adalah kelanjutan dari jawaban saya di sini .)
rstan
Perkiraan
Untuk referensi, ini adalah kode rstan.
raftery.model <- "
data{
int I;
int y[I];
}
parameters{
real<lower=max(y)> N;
simplex[2] theta;
}
transformed parameters{
}
model{
vector[I] Pr_y;
for(i in 1:I){
Pr_y[i] <- binomial_coefficient_log(N, y[i])
+multiply_log(y[i], theta[1])
+multiply_log((N-y[i]), theta[2]);
}
increment_log_prob(sum(Pr_y));
increment_log_prob(-log(N));
}
"
raft.data <- list(y=c(53,57,66,67,72), I=5)
system.time(fit.test <- stan(model_code=raftery.model, data=raft.data,iter=10))
system.time(fit <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))
Perhatikan bahwa saya berperan theta
sebagai 2-simpleks. Ini hanya untuk kesederhanaan. Jumlah bunga adalah theta[1]
; jelas theta[2]
informasi yang berlebihan.
Selain itu, adalah nilai nyata ( rstan
hanya menerima parameter bernilai riil karena merupakan metode gradien), jadi saya menulis distribusi binomial bernilai nyata.
Hasil Rstan
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
N 1078.75 256.72 15159.79 94.44 148.28 230.61 461.63 4575.49 3487 1
theta[1] 0.29 0.00 0.19 0.01 0.14 0.27 0.42 0.67 2519 1
theta[2] 0.71 0.00 0.19 0.33 0.58 0.73 0.86 0.99 2519 1
lp__ -19.88 0.02 1.11 -22.89 -20.31 -19.54 -19.09 -18.82 3339 1
Perkiraan Kisi
Perkiraan kisi diproduksi seperti di bawah ini. Batasan memori membuat saya tidak membuat kisi yang lebih bagus di laptop saya.
theta <- seq(0+1e-10,1-1e-10, len=1e3)
N <- round(seq(72, 5000, len=1e3)); N[2]-N[1]
grid <- expand.grid(N,theta)
y <- c(53,57,66,67,72)
raftery.prob <- function(x, z=y){
N <- x[1]
theta <- x[2]
exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N
}
post <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F)
post.norm <- post/sum(post)
Saya menggunakan pendekatan grid untuk menghasilkan tampilan kepadatan posterior ini. Kita dapat melihat bahwa posterior berbentuk pisang; posterior semacam ini dapat menjadi masalah untuk metrik HMC euclidian. (Tingkat keparahan bentuk pisang sebenarnya ditekan di sini karena pada skala log.) Jika Anda berpikir tentang bentuk pisang selama satu menit, Anda akan menyadari bahwa itu harus terletak pada garis . (Selain itu, perkiraan kisi-kisi yang ditampilkan dalam grafik ini tidak dinormalisasi untuk alasan kejelasan - selain itu pisang agak terlalu sempit untuk terlihat jelas.)
Hasil perkiraan grid
do.call(cbind, lapply(c(0.025, .25, .5, .75, .975), function(quantile){
approx(y=N, x=cumsum(rowSums(post.norm))/sum(post.norm), xout=quantile)
}))
[,1] [,2] [,3] [,4] [,5]
x 0.025 0.25 0.5 0.75 0.975
y 92.55068 144.7091 226.7845 443.6359 2475.398
Diskusi
Kuantil 97,5% untuk jauh lebih besar dalam model saya daripada untuk pendekatan grid, tetapi kuantilasinya mirip dengan perkiraan grid. Saya menafsirkan ini sebagai menunjukkan bahwa kedua metode umumnya sepakat. Saya tidak tahu bagaimana menafsirkan perbedaan dalam kuantil 97,5%, meskipun.rstan
Saya telah mengembangkan beberapa penjelasan yang mungkin untuk apa yang mungkin menjelaskan perbedaan antara perkiraan grid dan hasil dari rstan
sampel HMC-NUTS, tetapi saya tidak yakin bagaimana memahami apakah satu, keduanya atau tidak ada penjelasan yang benar.
- Rstan salah dan kisi sudah benar. Kepadatan berbentuk pisang bermasalah untuk
rstan
, terutama karena melayang ke arah , sehingga jumlah ekor ini tidak dapat dipercaya. Kita bisa melihat dari plot posterior atas grid yang ekor sangat tajam pada nilai-nilai yang lebih besar . - Rstan benar dan grid salah. Grid membuat dua perkiraan yang dapat merusak hasil. Pertama, kisi-kisi itu hanya sekumpulan titik hingga subruang posterior, jadi ini merupakan perkiraan kasar. Kedua, karena itu adalah ruang bagian yang terbatas, kita palsu menyatakan ada menjadi 0 probabilitas posterior atas nilai-nilai lebih besar dari nilai jaringan terbesar kami untuk . Demikian juga, lebih baik masuk ke ekor kisi-kisi, sehingga quanitles ekornya benar.
rstan
Saya membutuhkan lebih banyak ruang untuk mengklarifikasi poin dari jawaban Juho. Jika saya mengerti dengan benar, kita dapat mengintegrasikan dari posterior untuk mendapatkan distribusi beta-binomial:
Dalam kasus kami, dan karena kami memiliki seragam sebelum . Saya percaya bahwa posterior seharusnya mana karena . Tapi ini tampaknya sangat menyimpang dari jawaban Juho. Di mana saya salah?
sumber