Saya memiliki masalah pengambilan sampel sederhana, di mana loop batin saya terlihat seperti:
v = sample_gamma(k, a)
di mana sample_gamma
sampel dari distribusi Gamma membentuk sampel Dirichlet.
Ini bekerja dengan baik, tetapi untuk beberapa nilai k / a, beberapa proses hilir perhitungan underflow.
Saya mengadaptasinya untuk menggunakan variabel ruang log:
v = log(sample_gamma(k, a))
Setelah mengadaptasi semua sisa program, ia bekerja dengan benar (setidaknya itu memberi saya hasil yang persis sama pada kasus uji). Namun, ini lebih lambat dari sebelumnya.
Apakah ada cara untuk langsung mencicipi tanpa menggunakan fungsi lambat seperti ? Saya mencoba googling untuk ini, tetapi saya bahkan tidak tahu apakah distribusi ini memiliki nama yang sama (log-gamma?).
sampling
gamma-distribution
luispedro
sumber
sumber
Jawaban:
Pertimbangkan parameter bentuk kecil dekat 0, seperti . Dalam rentang antara 0 dan , kira-kira , jadi pdf Gamma kira-kira . Ini dapat diintegrasikan ke CDF perkiraan, . Menolaknya, kita melihat kekuatan : eksponen besar. Untuk ini menyebabkan kemungkinan underflow (nilai presisi ganda kurang dari , lebih atau kurang). Berikut adalah plot peluang untuk mendapatkan underflow sebagai fungsi dari basis-sepuluh logaritmaα α=1/100 α e−α 1 xα−1dx/Γ(α) Fα(x)=xααΓ(α) 1/α α=1/100 10−300 α :
Salah satu solusinya adalah dengan mengeksploitasi perkiraan ini untuk menghasilkan varian log (Gamma): pada dasarnya, cobalah untuk menghasilkan variate Gamma dan jika terlalu kecil, buat logaritma dari perkiraan distribusi daya ini (seperti yang ditunjukkan di bawah). (Lakukan ini berulang kali hingga log berada dalam kisaran underflow, sehingga ini merupakan pengganti yang valid untuk varian underflow yang asli.) Untuk perhitungan Dirichlet, kurangi maksimum semua logaritma dari masing-masing nilai log: ini secara implisit mengubah semua variasi Gamma sehingga tidak akan memengaruhi nilai Dirichlet. Perlakukan setiap log yang dihasilkan yang terlalu kecil (katakanlah, kurang dari -100) sebagai log dengan nol sebenarnya. Exponentiate log lainnya. Sekarang Anda dapat melanjutkan tanpa underflow.
Ini akan memakan waktu lebih lama dari sebelumnya, tetapi setidaknya itu akan berhasil!
Untuk menghasilkan perkiraan log Gamma variate dengan parameter bentuk , precompute . Ini mudah, karena ada algoritma untuk menghitung nilai log Gamma secara langsung . Buat float acak seragam antara 0 dan 1, ambil logaritma, bagi dengan , dan tambahkan ke sana.α C=log(Γ(α))+log(α) α C
Karena parameter skala hanya mengubah skala variasinya, tidak ada masalah mengakomodasinya dalam prosedur ini. Anda bahkan tidak memerlukannya jika semua parameter skala sama.
Edit
Di balasan lain OP menjelaskan metode di mana kekuatan dari varian seragam (a variate) dikalikan dengan . Ini berfungsi karena pdf dari distribusi bersama kedua varian ini sama dengan . Untuk menemukan pdf dari kami mengganti , dibagi dengan Jacobean , dan mengintegrasikan . Integral harus berkisar dari hingga karena , dari mana1/α B(α) Γ(α+1) (αxα−1)(yαe−ydy/Γ(α+1)) z=xy y→z/x x x z ∞ 0≤y≤1
yang merupakan pdf dari distribusi .Γ(α)
Intinya adalah bahwa ketika , nilai yang diambil dari tidak mungkin di-underflow dan dengan menjumlahkan log-nya dan kali log dari varian seragam independen kami akan memiliki log dari . Log cenderung sangat negatif, tetapi kita akan melewati konstruksi antilognya, yang akan melimpah dalam representasi titik mengambang.0<α<1 Γ(α+1) 1/α Γ(α)
sumber
Saya menjawab pertanyaan saya sendiri, tetapi saya menemukan solusi yang cukup bagus, bahkan jika saya tidak sepenuhnya memahaminya. Melihat kode dari Perpustakaan Ilmiah GNU, di sini adalah bagaimana sampel variabel gamma (α β
r
adalah nomor acak generator,a
adalah dan merupakan ):b
gsl_ran_gamma
adalah fungsi yang mengembalikan sampel gamma acak (jadi di atas adalah panggilan rekursif), sementaragsl_rng_uniform_pos
mengembalikan nomor yang terdistribusi secara seragam dalam ( untuk positif ketat karena dijamin tidak mengembalikan 0,0)._pos
Oleh karena itu, saya dapat mengambil log dari ekspresi terakhir dan menggunakannya
Untuk mendapatkan apa yang saya inginkan. Saya sekarang memiliki dua1/a 1/a
log()
panggilan (tetapi satu lebih sedikitpow()
), tetapi hasilnya mungkin lebih baik. Sebelumnya, seperti yang ditunjukkan oleh whuber, saya memiliki sesuatu yang dinaikkan menjadi kekuatan , berpotensi banyak. Sekarang, di logspace, saya mengalikan dengan . Jadi, kecil kemungkinannya untuk underflow.1 / asumber