Bagaimana cara mengambil sampel dari

Saya ingin sampel sesuai dengan kepadatan

f (a) \propto \frac{c^{a} d^{a - 1}}{Γ (a)} 1_{(1, \infty)} (a)

$f(a) \propto \frac{c^a d^{a-1}}{\Gamma(a)} 1_{(1,\infty)}(a)$ manadanbenar-benar positif. (Motivasi: Ini bisa berguna untuk pengambilan sampel Gibbs ketika parameter bentuk kepadatan Gamma memiliki seragam sebelumnya.)

c

$c$

d

$d$

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 , sampel dari seragam dalam kotak besar 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 besar $a^*$ $f$ $(a,u)$ $[0,10a^*]\times [0,f(a^*)]$ $u>f(a)$ $f(a^*)$ $c$ $d$ $c$ dan kira-kira pada .) $d$ $a=cd$

Terima kasih sebelumnya atas bantuannya!

+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)

$\Gamma(a)$

Γ

$\Gamma$

Γ

$\Gamma$

Ada perbedaan besar antara case

dan

. Apakah Anda perlu membahas kedua kasus ini?

c d < 2

$c d \lt 2$

c d \geq 2

$c d \ge 2$

Whuber

Itu benar - terima kasih. Kita dapat mengasumsikan bahwa

c d \geq 2

$cd\geq 2$

Jawaban:

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

Untuk sedikit menyederhanakan matematika, misalkan , tulis , dan catat itu $k = c d$ $x = a$

f (x) \propto \frac{k^{x}}{Γ (x)} d x

$f(x) \propto \frac{k^x}{\Gamma(x)} dx$

untuk . Pengaturan memberikan $x \ge 1$ $x = u^{3/2}$

f (u) \propto \frac{k^{u^{3 / 2}}}{Γ (u^{3 / 2})} u^{1 / 2} d u

$f(u) \propto \frac{k^{u^{3/2}}}{\Gamma(u^{3/2})} u^{1/2} du$

untuk . Ketika , distribusi ini sangat dekat dengan Normal (dan semakin dekat dengan semakin besar). Secara khusus, Anda bisa $u \ge 1$ $k \ge \exp(5)$ $k$

Temukan mode numerik (menggunakan, misalnya, Newton-Raphson). $f(u)$
Luaskan ke urutan kedua tentang modenya. $\log{f(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 , Anda mungkin perlu meningkatkan pdf Normal sedikit untuk memastikan dominasi.) $f(u)$ $k \lt \exp(5)$

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

Gambarkan nilai dari distribusi Normal yang mendominasi . $u$ $g(u)$
Jika atau jika baru seragam variate melebihi , kembali ke langkah 1. $u \lt 1$ $X$ $f(u)/(M g(u))$
Set . $x = u^{3/2}$

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

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 dimasukkan untuk memastikan logaritma positif di seluruh bagian utama distribusi; yaitu, untuk memastikan kecuali mungkin di daerah dengan probabilitas yang dapat diabaikan. Dengan membuat cukup besar Anda dapat menjamin bahwa $\log(\exp(0.004)g(u)/f(u))$ $M = \exp(0.004)$ $Mg(u) \ge f(u)$ $M$ $M \cdot g$ mendominasi dalam semua kecuali ekor yang paling ekstrim (yang praktis tidak memiliki peluang untuk dipilih dalam simulasi). Namun, semakin besar , semakin sering penolakan akan terjadi. Ketika tumbuh besar, dapat dipilih sangat dekat dengan , yang praktis tidak ada penalti. $f$ $M$ $k$ $M$ $1$

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

Plot for k=2

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

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 meningkat melebihi 5. $f(u)$ $k = c d$

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

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 = u^{3/2}$

Mengingat teori transformasi Box-Cox , adalah wajar untuk mencari beberapa transformasi kekuatan dari bentuk (untuk 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 ketiga $x = u^\alpha$ $\alpha$ $\alpha$ 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) = \frac{k^{u^{α}}}{Γ (u^{α})} u^{α - 1} .

$f(u) = \frac{k^{u^{\alpha}}}{\Gamma(u^{\alpha})} u^{\alpha-1}.$

Ambil logaritma dan menggunakan ekspansi asimtotik Stirling dari : $\log(\Gamma)$

\log (f (u)) \approx \log (k) u^{α} + (α - 1) \log (u) - α u^{α} \log (u) + u^{α} - \log (2 π u^{α}) / 2 + c u^{- α}

$\log(f(u)) \approx \log(k) u^\alpha + (\alpha - 1)\log(u) - \alpha u^\alpha \log(u) + u^\alpha - \log(2 \pi u^\alpha)/2 + c u^{-\alpha}$

(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$ $\alpha$

Hitung turunan ketiganya (yang bila dibagi , akan menjadi koefisien kekuatan ketiga 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 ) $3!$ $u$ $c$

- \frac{1}{2} u^{- (3 + α)} α (2 α (2 α - 3) u^{2 α} + (α^{2} - 5 α + 6) u^{α} + 12 c α) .

$-\frac{1}{2} u^{-(3+\alpha)} \alpha \left(2 \alpha(2 \alpha-3) u^{2 \alpha} + (\alpha^2 - 5\alpha +6)u^\alpha + 12 c \alpha \right).$

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

2 α - 3 = 0.

$2 \alpha-3 = 0.$

Itu sebabnya karya dengan baik: dengan pilihan ini, koefisien istilah kubik sekitar berperilaku puncak seperti , yang dekat dengan . Begitu melebihi 10 atau lebih, Anda bisa melupakannya, dan itu cukup kecil bahkan untuk turun ke 2. Kekuatan yang lebih tinggi, mulai dari yang keempat, memainkan peran yang semakin sedikit ketika bertambah besar, karena koefisien mereka tumbuh secara proporsional lebih kecil juga. Kebetulan, perhitungan yang sama (berdasarkan turunan kedua dari $\alpha = 3/2$ $u^{-3}$ $\exp(-2 k)$ $k$ $k$ $k$ pada puncaknya) menunjukkan standar deviasi dari perkiraan Normal ini sedikit kurang dari $log(f(u))$ , dengan error sebanding dengan. $\frac{2}{3}\exp(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

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 / (10 N)$

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

$1$

Γ (\exp (c d)) / Γ (x)

$\Gamma(\exp(c d))/\Gamma(x)$

x

$x$

\exp (k)

$\exp(k)$

Γ

$\Gamma$

1

$1$ dan

2

$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.

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

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

c d

$cd$ is large, it still breaks. This is basically because it has to calculate the equivalent of

(c d)^{x}

$(cd)^{x}$ for large

x

$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

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 dekat

f

$f$ daripada berkenaan dengan distribusi seragam. Mengenai perhitungan: hitung

a \log (c d) - \log (Γ (a))

$a \log(c d) - \log(\Gamma(a))$ (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

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

c d

$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