Bagaimana cara mengambil sampel dari

19

Saya ingin sampel sesuai dengan kepadatan

f(a)cada1Γ(a)1(1,)(a)
manadanbenar-benar positif. (Motivasi: Ini bisa berguna untuk pengambilan sampel Gibbs ketika parameter bentuk kepadatan Gamma memiliki seragam sebelumnya.)cd

Adakah yang tahu cara mengambil sampel dari kepadatan ini dengan mudah? Mungkin itu standar dan hanya sesuatu yang tidak saya ketahui?

Saya bisa memikirkan algoritma penolakan penolakan bodoh yang akan lebih atau kurang bekerja (menemukan mode dari f , sampel ( a , u ) dari seragam dalam kotak besar [ 0 , 10 a ] × [ 0 , f ( a ) ] dan menolak jika ), tetapi (i) sama sekali tidak efisien dan (ii) akan terlalu besar untuk ditangani oleh komputer dengan mudah bahkan untuk dan . (Perhatikan bahwa mode untuk besaraf(a,u)[0,10a]×[0,f(a)]u>f(a)f(a)cdcdan kira-kira pada .)da=cd

Terima kasih sebelumnya atas bantuannya!

NF
sumber
+1 pertanyaan bagus. Saya tidak yakin apakah ada pendekatan standar.
suncoolsu
Sudahkah Anda memeriksa (untuk gagasan) di tempat "jelas", seperti, misalnya, teks Devroye ?
kardinal
Ya, saya sudah mencoba sejumlah ide dari teks Devroye. The telah membuatnya menjadi sulit bagi saya untuk mendapatkan di mana saja dengan sebagian besar dari mereka, meskipun ... paling pendekatan tampaknya memerlukan integrasi baik (untuk menemukan cdf yang), dekomposisi ke dalam fungsi sederhana, atau melompat-lompat dengan fungsi sederhana ... tetapi fungsi Γ membuat semua ini sulit. Jika ada yang punya ide tentang di mana mencari pendekatan untuk subproblem ini - misalnya, di mana lagi fungsi Γ muncul dalam cara "penting" seperti di sini (bukan hanya sebagai konstanta normalisasi) dalam statistik - yang bisa sangat membantu saya ! Γ(a)ΓΓ
NF
Ada perbedaan besar antara case dan c d 2 . Apakah Anda perlu membahas kedua kasus ini? cd<2cd2
Whuber
1
Itu benar - terima kasih. Kita dapat mengasumsikan bahwa . cd2
NF

Jawaban:

21

Sampel penolakan akan bekerja dengan sangat baik ketika dan masuk akal untuk c d exp ( 2 ) .cdexp(5)cdexp(2)

Untuk sedikit menyederhanakan matematika, misalkan , tulis x = a , dan catat ituk=cdx=a

f(x)kxΓ(x)dx

untuk . Pengaturan x = u 3 / 2 memberikanx1x=u3/2

f(u)ku3/2Γ(u3/2)u1/2du

untuk . Ketika k exp ( 5 ) , distribusi ini sangat dekat dengan Normal (dan semakin dekat dengan k semakin besar). Secara khusus, Anda bisau1kexp(5)k

  1. Temukan mode numerik (menggunakan, misalnya, Newton-Raphson).f(u)

  2. Luaskan ke urutan kedua tentang modenya.logf(u)

Ini menghasilkan parameter dari distribusi Normal yang mendekati perkiraan. Untuk akurasi yang tinggi, Normal yang mendekati ini mendominasi kecuali pada ekor yang ekstrem. (Ketika k < exp ( 5 ) , Anda mungkin perlu meningkatkan pdf Normal sedikit untuk memastikan dominasi.)f(u)k<exp(5)

Setelah melakukan pekerjaan pendahuluan ini untuk setiap nilai diberikan , dan setelah memperkirakan konstanta M > 1 (seperti yang dijelaskan di bawah), memperoleh varian acak adalah masalah:kM>1

  1. Gambarkan nilai dari distribusi Normal yang mendominasi g ( u ) .ug(u)

  2. Jika atau jika baru seragam variate X melebihi f ( u ) / ( M g ( u ) ) , kembali ke langkah 1.u<1Xf(u)/(Mg(u))

  3. Set .x=u3/2

Jumlah yang diharapkan dari evaluasi karena perbedaan antara g dan f hanya sedikit lebih besar dari 1. (Beberapa evaluasi tambahan akan terjadi karena penolakan varian kurang dari 1 , tetapi bahkan ketika k serendah 2 frekuensi seperti kejadiannya kecil.)fgf1k2

Plot of f and g for k=5

Plot ini menunjukkan dengan logaritma dari g dan f sebagai fungsi u untuk . Karena grafiknya sangat dekat, kita perlu memeriksa rasio mereka untuk melihat apa yang terjadi:k=exp(5)

plot of log ratio

Ini menampilkan log rasio ; faktor M = exp ( 0,004 ) dimasukkan untuk memastikan logaritma positif di seluruh bagian utama distribusi; yaitu, untuk memastikan M g ( u ) f ( u ) kecuali mungkin di daerah dengan probabilitas yang dapat diabaikan. Dengan membuat M cukup besar Anda dapat menjamin bahwa M glog(exp(0.004)g(u)/f(u))M=exp(0.004)Mg(u)f(u)M.Mgmendominasi dalam semua kecuali ekor yang paling ekstrim (yang praktis tidak memiliki peluang untuk dipilih dalam simulasi). Namun, semakin besar M , semakin sering penolakan akan terjadi. Ketika k tumbuh besar, M dapat dipilih sangat dekat dengan 1 , yang praktis tidak ada penalti.fM.kM.1

Pendekatan serupa bekerja bahkan untuk , tetapi nilai-nilai M yang cukup besar mungkin diperlukan ketika exp ( 2 ) < k < exp ( 5 ) , karena f ( u ) terasa asimetris. Misalnya, dengan k = exp ( 2 ) , untuk mendapatkan g yang cukup akurat kita perlu menetapkan M = 1 :k>exp(2)M.exp(2)<k<exp(5)f(u)k=exp(2)gM=1

Plot for k=2

Kurva merah atas adalah grafik sedangkan kurva biru bawah adalah grafik log ( f ( u ) ) . Sampel penolakan f relatif terhadap exp ( 1 ) g akan menyebabkan sekitar 2/3 dari semua penarikan uji coba ditolak, tiga kali lipat upaya: masih tidak buruk. Ekor kanan ( u > 10 atau x > 10 3 / 2 ~ 30log(exp(1)g(u))log(f(u))fexp(1)gu>10x>103/230) akan kurang terwakili dalam sampel penolakan (karena tidak lagi mendominasi f di sana), tetapi ekor itu kurang dari exp ( - 20 ) 10 - 9 dari total probabilitas.exp(1)gfexp(20)109

Untuk meringkas, setelah upaya awal untuk menghitung mode dan mengevaluasi istilah kuadrat dari rangkaian daya sekitar mode - upaya yang paling banyak membutuhkan beberapa evaluasi fungsi - Anda dapat menggunakan sampel penolakan di biaya yang diharapkan antara 1 dan 3 (atau lebih) evaluasi per varian. Pengganda biaya dengan cepat turun menjadi 1 karena k = c d meningkat melebihi 5.f(u)k=cd

Bahkan ketika hanya satu gambar dari diperlukan, metode ini masuk akal. Ia muncul dengan sendirinya ketika banyak undian independen diperlukan untuk nilai k yang sama , karena kemudian perhitungan awal diamortisasi selama banyak undian.fk


Tambahan

@ Cardinal telah meminta, cukup masuk akal, untuk mendukung beberapa analisis melambaikan tangan di masa depan. Secara khusus, mengapa harus transformasi make distribusi sekitar Normal?x=u3/2

Mengingat teori transformasi Box-Cox , adalah wajar untuk mencari beberapa transformasi kekuatan dari bentuk (untuk α yang konstan , semoga tidak terlalu berbeda dari persatuan) yang akan membuat distribusi "lebih" normal. Ingatlah bahwa semua distribusi Normal hanya ditandai: logaritma pdf mereka murni kuadratik, dengan nol istilah linear dan tanpa syarat pemesanan lebih tinggi. Oleh karena itu kita dapat mengambil pdf apa pun dan membandingkannya dengan distribusi normal dengan memperluas logaritma sebagai rangkaian daya di sekitar puncaknya (tertinggi). Kami mencari nilai α yang membuat (setidaknya) yang ketigax=uαααkekuatan lenyap, setidaknya kira-kira: itulah yang paling bisa kita harapkan bahwa satu koefisien bebas akan tercapai. Seringkali ini bekerja dengan baik.

Tetapi bagaimana cara menangani distribusi khusus ini? Setelah melakukan transformasi kekuatan, pdf-nya adalah

f(u)=kuαΓ(uα)uα1.

Ambil logaritma dan menggunakan ekspansi asimtotik Stirling dari :log(Γ)

log(f(u))log(k)uα+(α1)log(u)αuαlog(u)+uαlog(2πuα)/2+cuα

(untuk nilai , yang tidak konstan). Ini berfungsi asalkan α positif, yang akan kita anggap sebagai kasus (karena kalau tidak kita tidak dapat mengabaikan sisa ekspansi).cα

Hitung turunan ketiganya (yang bila dibagi , akan menjadi koefisien kekuatan ketiga u dalam seri daya) dan manfaatkan fakta bahwa pada puncaknya, turunan pertama harus nol. Ini sangat menyederhanakan turunan ketiga, memberi (kira-kira, karena kita mengabaikan turunan c )3!uc

12u(3+α)α(2α(2α3)u2α+(α25α+6)uα+12cα).

Ketika tidak terlalu kecil, kamu memang akan besar di puncak. Karena α positif, istilah dominan dalam ungkapan ini adalah kekuatan 2 α , yang dapat kita atur menjadi nol dengan membuat koefisiennya menghilang:kuα2α

2α3=0.

Itu sebabnya karya dengan baik: dengan pilihan ini, koefisien istilah kubik sekitar berperilaku puncak seperti u - 3 , yang dekat dengan exp ( - 2 k ) . Begitu k melebihi 10 atau lebih, Anda bisa melupakannya, dan itu cukup kecil bahkan untuk k turun ke 2. Kekuatan yang lebih tinggi, mulai dari yang keempat, memainkan peran yang semakin sedikit ketika k bertambah besar, karena koefisien mereka tumbuh secara proporsional lebih kecil juga. Kebetulan, perhitungan yang sama (berdasarkan turunan kedua dari l o g ( fα=3/2u3exp(2k)kkk pada puncaknya) menunjukkan standar deviasi dari perkiraan Normal ini sedikit kurang dari 2log(f(u)), dengan error sebanding denganexp(-k/2).23exp(k/6)exp(k/2)

whuber
sumber
(+1) Jawaban bagus. Mungkin Anda dapat memperluas secara singkat pada motivasi untuk pilihan Anda dari variabel transformasi.
kardinal
Tambahan yang bagus. Ini membuat jawaban yang sangat, sangat lengkap!
kardinal
11

Saya sangat suka jawaban @ whuber; sepertinya sangat efisien dan memiliki analisis yang indah. Tapi itu membutuhkan beberapa wawasan yang mendalam sehubungan dengan distribusi khusus ini. Untuk situasi di mana Anda tidak memiliki wawasan tersebut (jadi untuk distribusi yang berbeda), saya juga menyukai pendekatan berikut yang berfungsi untuk semua distribusi di mana PDF dua kali dapat dibedakan dan turunan kedua memiliki banyak akar. Ini membutuhkan sedikit kerja untuk menyiapkan, tetapi kemudian setelah itu Anda memiliki mesin yang berfungsi untuk sebagian besar distribusi yang dapat Anda gunakan.

Pada dasarnya, idenya adalah menggunakan batas atas linear piecewise untuk PDF yang Anda adaptasi saat Anda melakukan rejection sampling. Pada saat yang sama Anda memiliki linear piecewise lebih rendahterikat untuk PDF yang mencegah Anda dari harus mengevaluasi PDF terlalu sering. Batas atas dan bawah diberikan oleh akor dan garis singgung ke grafik PDF. Pembagian awal ke dalam interval adalah sedemikian rupa sehingga pada setiap interval, PDF adalah semua cekung atau semua cembung; setiap kali Anda harus menolak suatu titik (x, y) Anda membagi interval itu di x. (Anda juga dapat melakukan subdivisi tambahan di x jika Anda harus menghitung PDF karena batas bawahnya benar-benar buruk.) perkiraan PDF Anda pada dasarnya gratis. Rinciannya sedikit rumit untuk mendapatkan hak, tapi saya sudah mencoba untuk menjelaskan sebagian besar dari mereka dalam seri dari blog posts - terutamayang terakhir .

Posting-posting itu tidak membahas apa yang harus dilakukan jika PDF tidak dibatasi baik dalam domain maupun dalam nilai; Saya akan merekomendasikan solusi yang agak jelas baik melakukan transformasi yang membuat mereka terbatas (yang akan sulit untuk diotomatisasi) atau menggunakan cutoff. Saya akan memilih cutoff tergantung pada jumlah total poin yang Anda harapkan untuk menghasilkan, katakan N , dan pilih cutoff sehingga bagian yang dihapus memiliki probabilitas kurang dari . (Ini cukup mudah jika Anda memiliki formulir tertutup untuk CDF; jika tidak, mungkin juga rumit.)1/(10N)

Metode ini diimplementasikan dalam Maple sebagai metode default untuk distribusi kontinu yang ditentukan pengguna. (Pengungkapan penuh - Saya bekerja untuk Maplesoft.)


Saya melakukan contoh menjalankan, menghasilkan 10 ^ 4 poin untuk c = 2, d = 3, menentukan [1, 100] sebagai rentang awal untuk nilai-nilai:

grafik

Ada 23 penolakan (merah), 51 poin "dalam masa percobaan" yang pada saat itu berada di antara batas bawah dan PDF yang sebenarnya, dan 9949 poin yang diterima setelah memeriksa hanya ketidaksetaraan linear. Itu 74 evaluasi dari total PDF, atau sekitar satu evaluasi PDF per 135 poin. Rasio akan menjadi lebih baik karena Anda menghasilkan lebih banyak poin, karena perkiraannya menjadi lebih baik dan lebih baik (dan sebaliknya, jika Anda hanya menghasilkan beberapa poin, rasionya lebih buruk).

Erik P.
sumber
Dan omong-omong - jika Anda perlu mengevaluasi PDF hanya sangat jarang karena Anda memiliki batas bawah yang baik untuk itu, Anda dapat mengambil lebih lama untuk itu, jadi Anda bisa menggunakan perpustakaan bignum (mungkin bahkan MPFR?) Dan mengevaluasi fungsi Gamma di dalamnya tanpa terlalu banyak rasa takut meluap.
Erik P.
(+1) Ini adalah pendekatan yang bagus. Terima kasih telah membagikannya.
whuber
1Γ(exp(cd))/Γ(x)xexp(k)Γ1 dan 2 ditambah sejumlah kecil logaritma: tidak ada overflow di sana.
whuber
@whuber re: Gammas: Ah ya - saya melihat bahwa Anda telah menyarankan ini juga. Terima kasih!
Erik P.
3

You could do it by numerically executing the inversion method, which says that if you plug uniform(0,1) random variables in the inverse CDF, you get a draw from the distribution. I've included some R code below that does this, and from the few checks I've done, it is working well, but it is a bit sloppy and I'm sure you could optimize it.

If you're not familiar with R, lgamma() is the log of the gamma function; integrate() calculates a definite 1-D integral; uniroot() calculates a root of a function using 1-D bisection.

# density. using the log-gamma gives a more numerically stable return for 
# the subsequent numerical integration (will not work without this trick)
f = function(x,c,d) exp( x*log(c) + (x-1)*log(d) - lgamma(x) )

# brute force calculation of the CDF, calculating the normalizing constant numerically
F = function(x,c,d) 
{
   g = function(x) f(x,c,d)
   return( integrate(g,1,x)$val/integrate(g,1,Inf)$val )
}

# Using bisection to find where the CDF equals p, to give the inverse CDF. This works 
# since the density given in the problem corresponds to a continuous CDF. 
F_1 = function(p,c,d) 
{
   Q = function(x) F(x,c,d)-p
   return( uniroot(Q, c(1+1e-10, 1e4))$root )
}

# plug uniform(0,1)'s into the inverse CDF. Testing for c=3, d=4. 
G = function(x) F_1(x,3,4)
z = sapply(runif(1000),G)

# simulated mean
mean(z)
[1] 13.10915

# exact mean
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*x/nc
integrate(h,1,Inf)$val
[1] 13.00002 

# simulated second moment
mean(z^2)
[1] 183.0266

# exact second moment
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*(x^2)/nc
integrate(h,1,Inf)$val
[1] 181.0003

# estimated density from the sample
plot(density(z))

# true density 
s = seq(1,25,length=1000)
plot(s, f(s,3,4), type="l", lwd=3)

The main arbitrary thing I do here is assuming that (1,10000) is a sufficient bracket for the bisection - I was lazy about this and there might be a more efficient way to choose this bracket. For very large values, the numerical calculation of the CDF (say, >100000) fails, so the bracket must be below this. The CDF is effectively equal to 1 at those points (unless c,d are very large), so something could probably be included that would prevent miscalculation of the CDF for very large input values.

Edit: When cd is very large, a numerical problem occurs with this method. As whuber points out in the comments, once this has occurred, the distribution is essentially degenerate at it's mode, making it a trivial sampling problem.

Macro
sumber
1
The method is correct, but awfully painful! How many function evaluations do you suppose are needed for a single random variate? Thousands? Tens of thousands?
whuber
There is a lot of computing, but it doesn't actually take very long - certainly much faster than rejection sampling. The simulation I showed above took less than a minute. The problem is that when cd is large, it still breaks. This is basically because it has to calculate the equivalent of (cd)x for large x. Any solution proposed will have that problem though - I'm trying to figure out if there's a way to do this on the log scale and transforming back.
Macro
1
Satu menit untuk 1.000 variasi tidak terlalu bagus: Anda akan menunggu berjam-jam untuk satu simulasi Monte-Carlo yang bagus. Anda dapat melakukan empat kali lipat lebih cepat menggunakan sampel penolakan. Caranya adalah dengan menolak dengan perkiraan dekatfdaripada berkenaan dengan distribusi seragam. Mengenai perhitungan: hitungSebuahcatatan(cd)-catatan(Γ(Sebuah))(dengan menghitung log Gamma secara langsung, tentu saja), kemudian secara eksponensial. Itu menghindari meluap.
whuber
That is what I do for the computation - it still doesn't avoid overflow. You can't exponentiate a number greater than around 500 on a computer. That quantity gets much larger than that. I mean "pretty good" comparing it with the rejection sampling the OP mentioned.
Makro
1
I did notice that the "standard deviation rule" that normals follow (68% within 1, 95% within 2, 99.7% within 3) did apply. So basically for large cd it's a point mass at the mode. From what you say, the threshold where this occurs before the numerical problems, so this still works. Thanks for the insight
Macro