Asumsikan saya memiliki fungsi yang ingin saya integrasikan
Karena , saya dapat mengganti dalam f ( x ) = U ( x ) / g ( x ) untuk membatalkan g dari integral, menghasilkan ekspresi bentuk 1 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
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:
1rnorm
Ini seharusnya dalam teori saya mengevaluasi ke1/√
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 ) = { 1 Saya tidak harus menggunakanU(x)sebagai kerapatan seragam ini. Saya bisa menggunakan beberapa kerapatan lain yang terintegrasi ke1, misalnya kerapatan probabilitas P(x)=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.
Jawaban:
Ini adalah pertanyaan yang paling menarik, yang berkaitan dengan masalah perkiraan konstanta normalisasi kepadatang 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
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 to1/π−−√ :
We discuss this method in details in two papers with Darren Wraith and with Jean-Michel Marin.
sumber