Memecahkan persamaan integral sederhana dengan pengambilan sampel acak

8

Biarkan menjadi fungsi non-negatif. Saya tertarik menemukan sedemikian rupa sehingga Peringatan itu : yang bisa saya lakukan adalah sampel pada titik di . Aku bisa, bagaimanapun, memilih lokasi di mana saya sampel secara acak, jika saya menginginkannya. fz[0,1]

0zf(x)dx=1201f(x)dx
f[0,1]f

Pertanyaan:

  1. Apakah mungkin untuk mendapatkan estimasi tidak bias setelah banyak sampel? Jika demikian, berapakah varian terkecil yang dapat dimiliki dari estimasi tersebut setelah sampel ?zk
  2. Jika tidak, prosedur apa yang tersedia untuk memperkirakan , dan apa waktu konvergensi yang terkait.z

Seperti yang ditunjukkan oleh Douglas Zare dalam komentar, ini bisa sangat sulit dilakukan jika fungsinya mendekati nol atau sangat besar. Untungnya, fungsi yang perlu saya gunakan untuk ini dibatasi dari atas dan bawah, jadi anggaplah . Selain itu, kita juga dapat mengasumsikan bahwa adalah Lipschitz atau bahkan dapat dibedakan jika itu membantu.1f(x)2f

Robinson
sumber
1
Jika Anda tidak memiliki informasi lebih lanjut, Anda dapat memiliki perilaku yang sangat buruk. Bayangkan bahwa adalahf0 antara 1/3 dan 2/3, dan 01/3f(x) dx1/2. Sedikit perubahan ke f akan membuat median lompatan dari bawah 1/3 ke atas 2/3.
Douglas Zare
@robinson Bisakah Anda memberikan informasi lebih lanjut tentang f? Atau apakah Anda tertarik untuk memecahkan masalah untuk kepadatan apa punf?
@DouglasZare - Terima kasih atas komentarnya; lihat hasil edit saya.
robinson
@Prastrastator - Saya mengedit pertanyaan dengan sedikit informasi lebih lanjut.
robinson
3
(+1) Untuk pembaruan. Membagi sisi kiri dengan kanan, orang dapat melihat bahwa ini mengurangi untuk menemukan median dari distribusi probabilitas yang tidak diketahui didukung pada[0,1].
kardinal

Jawaban:

6

Seperti yang ditunjukkan kardinal dalam komentarnya, pertanyaan Anda dapat disajikan kembali sebagai berikut.

Dengan aljabar sederhana, persamaan integral dapat ditulis ulang sebagai

0zg(x)dx=12,
di mana g adalah fungsi kepadatan probabilitas didefinisikan sebagai
g(x)=f(x)01f(t)dt.

Membiarkan X menjadi variabel acak dengan kerapatan g. Menurut definisi,P{Xz}=0zg(x)dx, jadi persamaan integral Anda setara dengan

P{Xz}=12,
yang berarti bahwa masalah Anda dapat dinyatakan sebagai:

"Membiarkan X menjadi variabel acak dengan kerapatan g. Temukan medianX. "

Untuk memperkirakan median X, gunakan metode simulasi apa pun untuk menggambar sampel nilai X dan ambil sebagai estimasi Anda median sampel.

Salah satu kemungkinan adalah menggunakan algoritma Metropolis-Hastings untuk mendapatkan sampel poin dengan distribusi yang diinginkan. Karena ekspresi probabilitas penerimaan dalam algoritma Metropolis-Hastings, kita tidak perlu mengetahui nilai konstanta normalisasi01f(t)dt kepadatan g. Jadi, kita tidak perlu melakukan integrasi ini.

Kode di bawah ini menggunakan bentuk sederhana dari algoritma Metropolis-Hastings yang dikenal sebagai Indepence Sampler, yang menggunakan proposal yang distribusinya tidak bergantung pada nilai rantai saat ini. Saya telah menggunakan proposal seragam independen. Sebagai perbandingan, skrip menampilkan minimum Monte Carlo dan hasilnya ditemukan dengan optimasi standar. Titik sampel disimpan dalam vektor chain, tetapi kami membuang yang pertama10000 titik yang membentuk apa yang disebut periode "terbakar" dalam simulasi.

BURN_IN = 10000
DRAWS   = 100000

f = function(x) exp(sin(x))

chain = numeric(BURN_IN + DRAWS)

x = 1/2

for (i in 1:(BURN_IN + DRAWS)) {
    y = runif(1) # proposal
    if (runif(1) < min(1, f(y)/f(x))) x = y
    chain[i] = x
}

x_min = median(chain[BURN_IN : (BURN_IN + DRAWS)])

cat("Metropolis minimum found at", x_min, "\n\n")

# MONTE CARLO ENDS HERE. The integrations bellow are just to check the results.

A = integrate(f, 0, 1)$value

F = function(x) (abs(integrate(f, 0, x)$value - A/2))

cat("Optimize minimum found at", optimize(F, c(0, 1))$minimum, "\n")

Berikut ini beberapa hasilnya:

Metropolis minimum found at 0.6005409 
Optimize minimum found at 0.601365

Kode ini dimaksudkan hanya sebagai titik awal untuk apa yang benar-benar Anda butuhkan. Karenanya, gunakan dengan hati-hati.

Zen
sumber
Terima kasih atas jawaban anda. Saya tidak tahu R, jadi saya tidak yakin bagaimana menguraikan apa yang Anda lakukan. Bisakah Anda menyatakan dengan kata-kata / formula prosedur Anda? Terima kasih. Secara khusus, saya bertanya-tanya apakah Anda menghargai batasan bahwa satu-satunya hal yang dapat Anda lakukan adalah mengevaluasi f - Anda tidak diperbolehkan, misalnya, untuk mengintegrasikanf, (walaupun Anda tentu saja dapat membentuk perkiraan Monte-Carlo untuk integral berdasarkan evaluasi acak).
robinson
Ya, saya hanya mengevaluasi funtuk mendapatkan estimasi Monte Carlo.
Zen
Kode hanyalah sebuah contoh. Sintaks R mirip dengan bahasa lain. Adakah pernyataan khusus yang tidak Anda mengerti? Lihatlah halaman Wikipedia pada algoritma Metropolis-Hastings. Tentu saja, gagasan umum lebih penting. Anda dapat mencicipi dari menuf/fmenggunakan metode apa pun yang Anda miliki.
Zen
Apakah Anda mengambil kursus pengantar tentang proses stokastik, yang mencakup rantai Markov waktu diskrit?
Zen
1
BTW: Penunda dunia, bersatu! Tapi tidak hari ini ...
Zen
3

Kualitas perkiraan integral, setidaknya dalam kasus sesederhana 1D, diberikan oleh (Teorema 2.10 dalam Niederreiter (1992) ):

|1Nn=1Nf(xn)01f(u)du|ω(f;DN(x1,,xN))
dimana
ω(f;t)=sup{|f(u)f(v)|:u,v[0,1],|uv|t,t>0}
adalah modulus fungsi kontinuitas (terkait dengan variasi total, dan mudah diekspresikan untuk fungsi Lipshitz), dan
DN(x1,,xN)=supu|1Nn1{xn[0,u)}u|=12N+maxn|xn2n12N|
adalah perbedaan (ekstrim), atau perbedaan maksimum antara fraksi hit dengan urutan x1,,xN dari interval semi-terbuka [0,u) dan ukuran Lebesgue-nya u. Ekspresi pertama adalah definisi, dan ekspresi kedua adalah properti dari urutan 1D di[0,1] (Teorema 2.6 dalam buku yang sama).

Jadi jelas untuk meminimalkan kesalahan dalam perkiraan integral, setidaknya dalam RHS persamaan Anda, Anda perlu mengambil xn=(2n1)/2N. Sekrup evaluasi acak, mereka berisiko memiliki kesenjangan acak pada fitur penting dari fungsi.

Kerugian besar dari pendekatan ini adalah Anda harus berkomitmen pada suatu nilai Nuntuk menghasilkan urutan yang didistribusikan secara seragam ini. Jika Anda tidak puas dengan kualitas perkiraan yang diberikannya, yang dapat Anda lakukan hanyalah menggandakan nilaiN dan tekan semua titik tengah dari interval yang dibuat sebelumnya.

Jika Anda ingin memiliki solusi di mana Anda dapat meningkatkan jumlah poin secara lebih bertahap, Anda dapat terus membaca buku itu, dan belajar tentang urutan van der Corput dan invers radikal. Lihat urutan perbedaan rendah di Wikipedia, ini memberikan semua detail.

Perbarui: untuk dipecahkanz, tentukan jumlah parsial

Sk=1Nn=1kf(2n12N).
Temukan k seperti yang
Sk12SN<Sk+1,
dan interpolasi untuk menemukan
zN=2k12N+SN/2SkN(Sk+1Sk).
Interpolasi ini mengasumsikan itu f()kontinu. Jika tambahanf() dua kali dapat dibedakan, maka pendekatan ini dengan mengintegrasikan ekspansi orde kedua untuk dimasukkan Sk1 dan Sk+2, dan memecahkan persamaan kubik untuk z.
Tugas
sumber
1
Saya suka intinya. Saya pikir akan sangat membantu untuk membuat lebih eksplisit strategi yang Anda usulkan untuk menyelesaikan pertanyaan OP. Saat ini, jawabannya berbunyi (untuk saya) sebagian besar seolah-olah membahas bagaimana menghitung RHS dari persamaan dalam pertanyaan.
kardinal
1
(+1) Pembaruan yang bagus. SN hanya dapat dilihat sebagai pendekatan Riemann-sum ke integral di mana kita menggunakan nilai fdi titik tengah setiap interval yang ditentukan oleh partisi, daripada titik akhir kiri atau kanan. :-)
cardinal
Iya; namun menarik bahwa jumlah Riemann ini memiliki pembenaran optimalitas ini.
Tugas