Dapatkah saya melakukan semiotomatis diagnostik konvergensi MCMC untuk mengatur panjang burn-in?

Saya ingin mengotomatiskan pilihan burn-in untuk rantai MCMC, misalnya dengan menghapus baris pertama berdasarkan pada diagnostik konvergensi.

Sejauh mana langkah ini dapat diotomatisasi dengan aman? Bahkan jika saya masih mengecek autokorelasi, jejak mcmc, dan pdf, alangkah baiknya jika memiliki pilihan burn-in length otomatis.

Pertanyaan saya bersifat umum, tetapi akan lebih bagus jika Anda bisa memberikan spesifik untuk berurusan dengan R mcmc.object; Saya menggunakan paket rjags dan coda di R.

Berikut ini satu pendekatan di otomatisasi. Umpan balik sangat dihargai. Ini adalah upaya untuk menggantikan inspeksi visual awal dengan perhitungan, diikuti dengan inspeksi visual berikutnya, sesuai dengan praktik standar.

Solusi ini sebenarnya menggabungkan dua solusi potensial, pertama, menghitung burn-in untuk menghapus panjang rantai sebelum beberapa ambang batas tercapai, dan kemudian menggunakan matriks autokorelasi untuk menghitung interval penjarangan.

1. menghitung vektor median maksimum faktor penyusutan diagnostik konvergensi Gelman-Rubin (grsf) untuk semua variabel di
2. menemukan jumlah sampel minimum di mana grsf di semua variabel berada di bawah ambang tertentu, misalnya 1,1 dalam contoh, mungkin lebih rendah dalam praktiknya
3. sub sampel rantai dari titik ini ke ujung rantai
4. menipiskan rantai menggunakan autokorelasi rantai yang paling berkorelasi otomatis
5. secara visual mengkonfirmasi konvergensi dengan jejak, autokorelasi, dan plot kerapatan

Objek mcmc dapat diunduh di sini: jags.out.Rdata

# jags.out is the mcmc.object with m variables
library(coda)    
load('jags.out.Rdata')
# 1. calculate max.gd.vec, 
# max.gd.vec is a vector of the maximum shrink factor
max.gd.vec     <- apply(gelman.plot(jags.out)$shrink[, ,'median'], 1, max)
# 2. will use window() to subsample the jags.out mcmc.object
# 3. start window at min(where max.gd.vec < 1.1, 100) 
window.start   <- max(100, min(as.numeric(names(which(max.gd.vec - 1.1 < 0)))))
jags.out.trunc <- window(jags.out, start = window.start)
# 4. calculate thinning interval
# thin.int is the chain thin interval
# step is very slow 
# 4.1 find n most autocorrelated variables
n = min(3, ncol(acm))
acm             <- autocorr.diag(jags.out.trunc)
acm.subset      <- colnames(acm)[rank(-colSums(acm))][1:n]
jags.out.subset <- jags.out.trunc[,acm.subset]
# 4.2 calculate the thinning interval
# ac.int is the time step interval for autocorrelation matrix
ac.int          <- 500 #set high to reduce computation time
thin.int        <- max(apply(acm2 < 0, 2, function(x) match(T,x)) * ac.int, 50)
# 4.3 thin the chain 
jags.out.thin   <- window(jags.out.trunc, thin = thin.int)
# 5. plots for visual diagnostics
plot(jags.out.thin)
autocorr.plot(jags.win.out.thin)

--memperbarui--

Seperti yang diterapkan dalam R, perhitungan matriks autokorelasi lebih lambat daripada yang diinginkan (> 15 menit dalam beberapa kasus), pada tingkat yang lebih rendah, demikian juga perhitungan faktor penyusutan GR. Ada pertanyaan tentang bagaimana mempercepat langkah 4 di stackoverflow di sini

--perbarui bagian 2--

jawaban tambahan:

1. Tidak mungkin untuk mendiagnosis konvergensi, hanya untuk mendiagnosis kurangnya konvergensi (Brooks, Giudici, dan Philippe, 2003)

2. Fungsi autorun.jags dari paket runjags mengotomatiskan perhitungan run run dan diagnosa konvergensi. Itu tidak mulai memantau rantai sampai diagnostik Gelman rubin di bawah 1,05; itu menghitung panjang rantai menggunakan diagnostik Raftery dan Lewis.

3. Gelman et al (Gelman 2004 Bayesian Data Analysis, hal. 295, Gelman dan Shirley, 2010 ) menyatakan bahwa mereka menggunakan pendekatan konservatif dengan membuang paruh pertama rantai. Meskipun merupakan solusi yang relatif sederhana, dalam praktiknya ini cukup untuk menyelesaikan masalah bagi serangkaian model dan data saya.

#code for answer 3
chain.length <- summary(jags.out)$end
jags.out.trunc <- window(jags.out, start = chain.length / 2)
# thin based on autocorrelation if < 50, otherwise ignore
acm <- autocorr.diag(jags.out.trunc, lags = c(1, 5, 10, 15, 25))
# require visual inspection, check acceptance rate
if (acm == 50) stop('check acceptance rate, inspect diagnostic figures') 
thin.int <- min(apply(acm2 < 0, 2, function(x) match(TRUE, x)), 50)
jags.out.thin <- window(jags.out.trunc, thin = thin.int)