Perhitungan kemungkinan marginal dari sampel MCMC

24

Ini adalah pertanyaan berulang (lihat posting ini , posting ini dan posting ini ), tetapi saya memiliki putaran berbeda.

Misalkan saya memiliki banyak sampel dari sampler MCMC generik. Untuk setiap sampel , saya tahu nilai kemungkinanlog f ( x | θ ) log f ( θ ) log f ( x i | θ )θlogf(x|θ) dan dari prior . Jika ini membantu, saya juga tahu nilai kemungkinan log per titik data, (informasi ini membantu dengan metode tertentu, seperti WAIC dan PSIS-LOO).logf(θ)logf(xi|θ)

Saya ingin mendapatkan perkiraan (kasar) dari kemungkinan marjinal, hanya dengan sampel yang saya miliki, dan mungkin beberapa evaluasi fungsi lainnya (tetapi tanpa menjalankan kembali MCMC ad hoc ).

Pertama-tama, mari kita membersihkan meja. Kita semua tahu bahwa estimator harmonik adalah estimator terburuk yang pernah ada . Mari kita lanjutkan. Jika Anda melakukan pengambilan sampel Gibbs dengan prior dan posterior dalam bentuk tertutup, Anda dapat menggunakan metode Chib ; tapi saya tidak yakin bagaimana menggeneralisasi di luar kasus-kasus itu. Ada juga metode yang mengharuskan Anda untuk memodifikasi prosedur pengambilan sampel (seperti via tempered posteriors ), tetapi saya tidak tertarik dengan itu di sini.

Pendekatan yang saya pikirkan terdiri dari perkiraan distribusi yang mendasari dengan bentuk parametrik (atau nonparametrik) g(θ) , dan kemudian mencari tahu konstanta normalisasi Z sebagai masalah optimasi 1-D (yaitu, Z yang meminimalkan beberapa kesalahan antara Zg(θ) dan f(x|θ)f(θ) , dievaluasi pada sampel). Dalam kasus paling sederhana, misalkan posterior kira-kira multivariat normal, saya dapat memasukkan g(θ) sebagai multivariat normal dan mendapatkan sesuatu yang mirip dengan perkiraan Laplace (saya mungkin ingin menggunakan beberapa evaluasi fungsi tambahan untuk memperbaiki posisi mode). Namun, saya dapat menggunakan sebagai g(θ)keluarga yang lebih fleksibel seperti campuran variational dari multivariat t distribusi.

Saya menghargai bahwa metode ini hanya berfungsi jika Zg(θ) adalah perkiraan yang masuk akal untuk f(x|θ)f(θ) , tetapi alasan atau dongeng peringatan mengapa sangat tidak bijaksana untuk lakukan? Adakah bacaan yang akan Anda rekomendasikan?

Pendekatan sepenuhnya nonparametrik menggunakan beberapa keluarga nonparametrik, seperti proses Gaussian (GP), untuk memperkirakan logf(x|θ)+logf(θ) (atau transformasi nonlinier lainnya, seperti sebagai akar kuadrat), dan quadrature Bayesian untuk secara implisit mengintegrasikan target yang mendasarinya (lihat di sini dan di sini ). Ini tampaknya menjadi pendekatan alternatif yang menarik, tetapi analog dalam semangat (juga, perhatikan bahwa dokter akan sulit dalam kasus saya).

Lacerbi
sumber
6
Saya pikir Chib, S. dan Jeliazkov, I. 2001 "Kemungkinan marjinal dari keluaran Metropolis - Hastings" digeneralisasi ke keluaran MCMC normal - akan tertarik untuk mendengar pengalaman dengan pendekatan ini. Adapun GP - pada dasarnya, ini bermuara pada persaingan posterior, yang Anda juga dapat mempertimbangkan untuk masalah lain. Saya kira masalahnya adalah bahwa Anda tidak pernah yakin tentang kualitas aproksimasi. Yang saya juga ingin tahu adalah apakah sampel MCMC ideal untuk model GP, atau apakah Anda harus berinvestasi lebih banyak di bagian ekor.
Florian Hartig
2
(+1) Terima kasih untuk referensi, terlihat tepat - Saya akan memeriksanya. Saya setuju bahwa semua pendekatan berbasis model bisa bermasalah (hal yang baik dengan quadrature Bayesian adalah Anda mendapatkan perkiraan ketidakpastian, meskipun tidak yakin seberapa dikalibrasi itu). Untuk saat ini, tujuan sederhana saya adalah melakukan sesuatu yang "lebih baik daripada perkiraan Laplace".
lacerbi

Jawaban:

26

Perpanjangan oleh Chib dan Jeliazkov (2001) sayangnya menjadi cepat mahal atau sangat bervariasi, yang merupakan alasan mengapa itu tidak banyak digunakan di luar kasus sampling Gibbs.

Walaupun ada banyak cara dan pendekatan untuk masalah estimasi konstanta normalisasi (seperti yang diilustrasikan oleh pembicaraan yang cukup beragam dalam lokakarya Estimating Constant yang kami jalankan minggu lalu di University of Warwick, slide tersedia di sana ), beberapa solusi memang mengeksploitasi langsung output MCMC .Z

  1. Seperti yang Anda sebutkan, penaksir rata-rata harmonik dari Newton dan Raftery (1994) hampir selalu buruk karena memiliki varian yang tak terbatas. Namun, ada cara untuk menghindari kutukan varians tak terbatas dengan menggunakan sebagai gantinya target dukungan yang terbatas dalam identitas rata-rata harmonik dengan memilihαsebagai indikator wilayah HPD untuk posterior. Ini memastikan varians yang terbatas dengan menghilangkan ekor di rata-rata harmonik. (Perinciannya dapat ditemukan dalammakalah yang saya tulis dengan Darren Wraithdan dalambab tentang normalisasi konstanta yangditulis dengan Jean-Michel Marin.) Singkatnya, metode ini mendaur ulang output MCMCθ1,,θMdengan mengidentifikasiβ( 20% mengatakan) nilai terbesar dari targetπ(θ)f(x|θ)dan menciptakanα

    α(θ)π(θ)f(x|θ)dπ(θ|x)=1Z
    αθ1,,θMβπ(θ)f(x|θ)αsebagai seragam atas persatuan bola berpusat pada orang-orang kepadatan terbesar (HPD) simulasi dan dengan radius ρ , berarti estimasi normalisasi konstan Z diberikan oleh Z - 1 = 1θi0ρZ jikadadalah dimensiθ(koreksi berlaku untuk bola berpotongan) dan jikaρcukup kecil untuk bola tidak pernah berpotongan (berarti bahwa hanya satu indikator terbaik pada bola adalah berbeda dari nol). Penjelasan untukpenyebutαM2adalah bahwa ini adalah jumlah ganda dariistilahβM2: 1
    Z^1=1βM2m=1Mdouble sum overβM ball centres θi0and M simulations θmI(0,ρ)(mini||θmθi0||){π(θm)f(x|θm)}1/πd/2ρdΓ(d/2+1)1volume of ball with radius ρβMα(θm)π(θm)f(x|θm)
    dθραM2βM2 dengan masing-masing istilah dalamθmmengintegrasikan keZ-1.
    1βMi=1βM1Mm=1MU(θi0,ρ)(θm)same as with min×1π(θm)f(x|θm)
    θmZ1
  2. Pendekatan lain adalah mengubah konstanta normalisasi menjadi parameter. Ini kedengarannya seperti bid'ah statistik tetapi makalah oleh Guttmann dan Hyvärinen (2012) meyakinkan saya akan hal sebaliknya. Tanpa terlalu banyak ke rincian, ide rapi di dalamnya adalah untuk mengubah diamati log-likelihood n Σ i = 1 f ( x i | q ) - n log exp f ( x | q ) d x menjadi log-likelihood bersama n i = 1 [ fZ

    i=1nf(xi|θ)nlogexpf(x|θ)dx
    yang merupakan log-kemungkinan proses titik Poisson dengan fungsi intensitas exp { f ( x | θ ) + ν + log n }
    i=1n[f(xi|θ)+ν]nexp[f(x|θ)+ν]dx
    exp{f(x|θ)+ν+logn}
    Ini adalah model alternatif di mana kemungkinan asli tidak muncul sebagai marginal di atas. Hanya mode yang bertepatan, dengan mode bersyarat di ν yang menyediakan konstanta normalisasi. Dalam praktiknya, kemungkinan proses Poisson di atas tidak tersedia dan Guttmann dan Hyvärinen (2012) menawarkan perkiraan dengan menggunakan regresi logistik. Untuk menghubungkan lebih baik dengan pertanyaan Anda, perkiraan Geyer adalah MLE, maka solusi untuk masalah maksimalisasi.
  3. π(θ|x)π(θ|x)g(θ)π(θ|x)g(θ)). Dengan regresi menjadi nilai-nilai dari kedua kepadatan, dinormalisasi atau tidak. Ini kebetulan terkait langsung dengan jembatan sampel Gelman dan Meng (1997), yang juga mendaur ulang sampel dari target yang berbeda. Dan versi selanjutnya, seperti Meng's MLE.
  4. Pendekatan berbeda yang memaksa seseorang untuk menjalankan sampler MCMC spesifik adalah sampling bersarang Skilling . Sementara saya [dan yang lain] memiliki beberapa keraguan mengenai efisiensi metode ini, metode ini cukup populer dalam astrostatistik dan kosmologi, dengan perangkat lunak yang tersedia seperti multinest .
  5. H0:θ=θ0ξπ1(θ)π2(ξ)H0
    B01(x)=πθ(θ0|x)π1(θ0)
    πθ(θ0|x)θθ0H0:θ=θ0
    m0(x)=Ξf(x|θ0,ξ)π2(ξ)dξ
    is available in closed form, one can derive the marginal density for the unconstrained model
    ma(x)=Θ×Ξf(x|θ,ξ)π1(θ)π2(ξ)dθdξ
    from the Bayes factor. (This Savage-Dickey representation relies on specific versions of three different densities and so is fraught with danger, not even mentioning the computational challenge of producing the marginal posterior.)

[Here is a set of slides I wrote about estimating normalising constants for a NIPS workshop last December.]

Xi'an
sumber
2
(+1) Incredibly rich answer, thank you. This will be useful to me and, I suppose, many other people. It will take me some time to have a look at the various approaches, and then I might come back with specific questions.
lacerbi
2
Starting from point (1)... I read the relevant articles. The "corrected" harmonic mean estimator seems exactly what I was looking for. It's neat and easy to compute given a MCMC output. So... what's the catch? It doesn't look like the method is being widely used, judging from a quick search on Google Scholar. What are its limitations? (besides the need to identify the HPD regions, which I imagine might become an issue for very complicated posteriors in high dimension). I am definitely going to give it a try -- but I wonder if there is something I need to be wary of.
lacerbi
2
Saya menambahkan beberapa rincian lebih lanjut: masalah dalam menerapkan seragam HPD adalah untuk mengetahui pendekatan kompak yang tepat untuk wilayah HPD. Cembung hull poin dengan nilai posterior tinggi (NP?) Sulit untuk menentukan sementara bola yang berpusat pada titik-titik tersebut dapat berpotongan, yang menciptakan masalah normalisasi konstan sekunder.
Xi'an
2
@Xi'an : very helpful, thanks! Can I ask: of all the mentioned approaches, what would currently be your recommendation if one looks for a general approach that tends to work out of the box (i.e. no tuning / checking required from the user)? I would be especially interested in the case of models with a low (< 50) number of parameters, non-normal posteriors, and strong correlations between parameters.
Florian Hartig
1
@FlorianHartig: the fact that a generic software like BUGS does not return a generic estimate of Zadalah semacam mengungkapkan sejauh mana masalah. Banyak solusi yang dapat ditemukan dalam literatur khusus belum menghasilkan perkiraan konsensus. Oleh karena itu, rekomendasi saya adalah memilih solusi regresi logistik Geyer, yang agak tidak sensitif terhadap dimensi.
Xi'an