Integrasi Metropolis-Hastings - mengapa strategi saya tidak berhasil?

16

Asumsikan saya memiliki fungsi g(x) yang ingin saya integrasikan

g(x)dx.
Tentu saja dengan asumsi g(x) menjadi nol pada titik akhir, tidak ada semburan, fungsi yang bagus. Salah satu cara yang saya telah mengutak-atik adalah dengan menggunakan algoritma Metropolis-Hastings untuk menghasilkan daftar sampel x1,x2,,xn dari distribusi proporsional untuk g(x) , yang hilang konstanta normalisasi
N=g(x)dx
yang akan saya sebutp(x) , dan kemudian menghitung beberapa statistikf(x) pada inix:
1ni=0nf(xi)f(x)p(x)dx.

Karena , saya dapat mengganti dalam f ( x ) = U ( x ) / g ( x ) untuk membatalkan g dari integral, menghasilkan ekspresi bentuk 1p(x)=g(x)/Nf(x)=U(x)/g(x)g Jadi asalkanU(x)terintegrasi ke1 disepanjang wilayah itu, saya harus mendapatkan hasil1/N, yang bisa saya ambil secara timbal balik untuk mendapatkan jawaban yang saya inginkan. Oleh karena itu saya dapat mengambil rentang sampel saya (untuk menggunakan titik-titik yang paling efektif)r=xmax-xmindan biarkanU(x)=1/runtuk setiap sampel yang telah saya gambar. Dengan begituU

1NU(x)g(x)g(x)dx=1NU(x)dx.
U(x)11/Nr=xmaxxminU(x)=1/rU(x)mengevaluasi nol di luar wilayah di mana sampel saya tidak, tetapi terintegrasi ke di wilayah itu. Jadi jika sekarang saya mengambil nilai yang diharapkan, saya harus mendapatkan: E [ U ( x )1
E[U(x)g(x)]=1N1ni=0nU(x)g(x).

Saya mencoba menguji ini dalam R untuk fungsi sampel . Dalam hal ini saya tidak menggunakan Metropolis-Hastings untuk menghasilkan sampel tetapi menggunakan probabilitas aktual untuk menghasilkan sampel (hanya untuk menguji). Saya tidak begitu mendapatkan hasil yang saya cari. Pada dasarnya ekspresi penuh dari apa yang saya hitung adalah: 1g(x)=ex2rnorm Ini seharusnya dalam teori saya mengevaluasi ke1/

1n(xmaxxmin)i=0n1exi2.
1/π . Mendekat tapi tentu saja tidak menyatu dengan cara yang diharapkan, apakah saya melakukan sesuatu yang salah?
ys = rnorm(1000000, 0, 1/sqrt(2))
r = max(ys) - min(ys)
sum(sapply(ys, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.6019741. 1/sqrt(pi) = 0.5641896

Edit untuk CliffAB

Alasan saya menggunakan rentang hanya untuk dengan mudah mendefinisikan fungsi yang bukan nol di atas wilayah di mana poin saya berada, tetapi itu terintegrasi ke pada rentang [ - , ] . Spesifikasi lengkap fungsi adalah: U ( x ) = { 11[,] Saya tidak harus menggunakanU(x)sebagai kerapatan seragam ini. Saya bisa menggunakan beberapa kerapatan lain yang terintegrasi ke1, misalnya kerapatan probabilitas P(x)=1

U(x)={1xmaxxminxmax>x>xmin0otherwise.
U(x)1 Namun ini akan membuat menjumlahkan sampel individu sepele yaitu 1
P(x)=1πex2.
1ni=0nP(x)g(x)=1ni=0nexi2/πexi2=1ni=0n1π=1π.

Saya bisa mencoba teknik ini untuk distribusi lain yang terintegrasi ke . Namun, saya masih ingin tahu mengapa itu tidak berhasil untuk distribusi yang seragam.1

Mike Flynn
sumber
Hanya dengan cepat melihat ini, jadi saya tidak yakin persis mengapa Anda memutuskan untuk menggunakan rentang (x). Syaratnya itu valid, itu sangat tidak efisien! Kisaran sampel dengan ukuran itu hanya sekitar statistik paling tidak stabil yang bisa Anda ambil.
Cliff AB
@CliffAB Tidak ada yang khusus tentang saya menggunakan rentang, selain mendefinisikan distribusi seragam pada interval di mana poin saya berada. Lihat hasil edit.
Mike Flynn
1
Saya akan melihat ini nanti secara lebih rinci. Tetapi sesuatu yang perlu dipertimbangkan adalah seolah-olah x adalah seperangkat RV seragam, maka sebagain, jarak(x)1. Tetapi jika x adalah seperangkat RV normal non-degenarate, maka sebagain, jarak(x).
Cliff AB
@CliffAB you might have been right, I think the reason was that the bounds of the integral were not fixed, and so the variance of the estimator will never converge...
Mike Flynn

Jawaban:

13

Ini adalah pertanyaan yang paling menarik, yang berkaitan dengan masalah perkiraan konstanta normalisasi kepadatan g berdasarkan pada output MCMC dari kepadatan yang sama g. (Komentar sampingan adalah bahwa asumsi yang benar untuk dibuat adalah itug dapat diintegrasikan, pergi ke nol tanpa batas tidak cukup.)

Menurut pendapat saya, entri yang paling relevan tentang topik ini sehubungan dengan saran Anda adalah makalah oleh Gelfand dan Dey (1994, JRSS B ), di mana penulis mengembangkan pendekatan yang sangat mirip untuk menemukan

Xg(x)dx
saat menghasilkan dari hal(x)g(x). Salah satu hasil dalam makalah ini adalah, untuk setiap probabilitas kepadatanα(x) [this is equivalent to your U(x)] such that
{x;α(x)>0}{x;g(x)>0}
the following identity
Xα(x)g(x)p(x)dx=Xα(x)Ndx=1N
shows that a sample from p can produce an unbiased evaluation of 1/N by the importance sampling estimator
η^=1ni=1nα(xi)g(xi)xiiidp(x)
Obviously, the performances (convergence speed, existence of a variance, &tc.) of the estimator η^ do depend on the choice of α [even though its expectation does not]. In a Bayesian framework, a choice advocated by Gelfand and Dey is to take α=π, the prior density. This leads to
α(x)g(x)=1(x)
where (x) is the likelihood function, since g(x)=π(x)(x). Unfortunately, the resulting estimator
N^=ni=1n1/(xi)
is the harmonic mean estimator, also called the worst Monte Carlo estimator ever by Radford Neal, from the University of Toronto. So it does not always work out nicely. Or even hardly ever.

Your idea of using the range of your sample (min(xi),max(xi)) and the uniform over that range is connected with the harmonic mean issue: this estimator does not have a variance if only because because of the exp{x2} appearing in the numerator (I suspect it could always be the case for an unbounded support!) and it thus converges very slowly to the normalising constant. For instance, if you rerun your code several times, you get very different numerical values after 10⁶ iterations. This means you cannot even trust the magnitude of the answer.

A generic fix to this infinite variance issue is to use for α a more concentrated density, using for instance the quartiles of your sample (q.25(xi),q.75(xi)), because g then remains lower-bounded over this interval.

When adapting your code to this new density, the approximation is much closer to 1/π:

ys = rnorm(1e6, 0, 1/sqrt(2))
r = quantile(ys,.75) - quantile(ys,.25)
yc=ys[(ys>quantile(ys,.25))&(ys<quantile(ys,.75))]
sum(sapply(yc, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.5649015. 1/sqrt(pi) = 0.5641896

We discuss this method in details in two papers with Darren Wraith and with Jean-Michel Marin.

Xi'an
sumber