Bagaimana saya bisa mengambil sampel dari distribusi dengan CDF yang tidak dapat dihitung?

Simulasi ilmu semi-komputer terkait masalah di sini.

Saya punya distribusi di mana

P (x) =$\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$

untuk beberapa konstanta b dan n, dan x adalah bilangan bulat sehingga .$0\leq x \leq n$

Sekarang, saya perlu sampel dari distribusi ini. Ini memiliki CDF yang dapat dibalik, jadi dimungkinkan untuk melakukan ini secara langsung dalam teori. Masalahnya adalah bahwa angka-angka yang terlibat adalah BESAR. Begitu besar pada kenyataannya, bahwa mereka berdua meluap variabel diformat secara konvensional, dan mengambil setidaknya beberapa menit (pada titik tertentu saya menyerah ...) untuk menghitung menggunakan format presisi sewenang-wenang. Pada dasarnya, CDF terbalik masih melibatkan istilah , untuk . Meskipun demikian, angka-angka output masih akan berada dalam kisaran , jadi sepertinya harus ada cara untuk melakukan ini.$e^{b(n+1)}$$350 < n < 3500$$0-n$

Apa yang saya cari adalah cara kira-kira pengambilan sampel dari distribusi ini yang dapat dihitung. Apakah ada metode pengambilan sampel alternatif? Apakah mereka?

CDF siap dibalik. Formula untuk inversi mengarah pada apa yang harus menjadi salah satu solusi yang paling sederhana dan paling bijaksana.

Mulailah dengan mengamati bahwa kemungkinan hasilnya $k$, $0 \le k \le n$, sebanding dengan $e^{-b k}$. Jadi, jika kita menghasilkan nilai yang seragam$q$ antara $0$ dan $q_{\max}=\sum_{k=0}^{n} e^{-b k}$= , kita hanya perlu menemukan terbesar yang$(1 - e^{-b(n+1)})/(1 - e^{-b})$$k$

Aljabar sederhana memberikan solusinya

Berikut ini adalah Rimplementasi yang dibangun seperti semua generator bilangan acak lainnya: argumen pertamanya menentukan berapa banyak nilai iid yang akan dihasilkan dan argumen lainnya menamai parameter ( as dan as ):$b$b$n$n.max

rgeom.truncated <- function(n=1, b, n.max) {
  a <- 1 - exp(-b)
  q.max <- (1 - exp(-b*(n.max+1))) / a
  q <- runif(n, 0, q.max)
  return(-ceiling(log(1 - q*a) / b))
}

Sebagai contoh penggunaannya, mari kita hasilkan satu juta variasi acak menurut distribusi ini:

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))

( detik diperlukan.)$0.10$

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)

( ditambahkan ke setiap nilai dalam rangka menciptakan histogram yang lebih baik: 's prosedur memiliki keistimewaan (= bug) di mana bar pertama adalah terlalu tinggi ketika endpoint kiri ditetapkan pada nol.) Kurva merah adalah distribusi referensi bahwa simulasi ini mencoba mereproduksi. Mari kita evaluasi kebaikan fit dengan uji chi-square:$1$Rhist

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)

Nilai p adalah : sangat pas.$0.84$

Anda sedang berhadapan dengan distribusi geometris terpotong dengan . Ada berbagai cara untuk mendekati ini.$p = 1-e^{-b}$

Saya akan menyarankan opsi yang berbeda dalam situasi yang berbeda; beberapa opsi akan melibatkan simulasi dari geometrik dan regenerasi ketika berada di luar rentang, mengambil bagian integer dari eksponensial terpotong yang sesuai ( seperti di sini ), atau menggunakan salah satu dari beberapa teknik cepat yang dirancang untuk distribusi diskrit pada rentang yang terbatas. Mengingat bahwa besar, mengambil dasar dari eksponensial terpotong cenderung relatif cepat, tetapi apakah itu pilihan terbaik juga tergantung pada .$n$$b$

Ini pertanyaan terkait tentang math.SE

Sebelum saya mencoba saran spesifik, apa kisaran nilai khas untuk ?$b$

Pertama, perhatikan bahwa yang, jika kontinu, akan terkait dengan distribusi eksponensial. Kemudian, apa yang dapat Anda lakukan adalah mensimulasikan dari distribusi eksponensial terpotong dan mengambil (bagian integer) pengamatan.$P(x)\propto e^{-bx}$$x$floor()

Cdf dari eksponensial terpotong adalah

Kemudian, jika kita membuat , kita memperoleh . Jika besar, maka yang menyarankan untuk mendekati .$F(x;n,b)=u$$x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$$bn$$e^{-bn}\approx 0$$x\approx -\dfrac{1}{b}\log[1-u]$

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

Cara untuk mengambil sampel dari target distribusi adalah dengan$p(k)\propto \exp\{-bk\}$

1. menjalankan eksperimen Metropolis-Hastings untuk menentukan (menarik) dukungan dari distribusi, yaitu di mana subset dari terkonsentrasi;$\{0,1,\ldots,n\}$

   metro=function(N,b,n){
   x=sample(0:n,N,rep=TRUE)
   for (t in 2:N){
     x[t]=prop=x[t-1]+sample(c(-1,1),1)
   
     if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
         x[t]=x[t-1]
     }
   return(x)
   }
   

2. Gunakan dukungan yang ditentukan, katakan, untuk menghitung probabilitas tepat sebagai untuk menghindari luberan.$\{k_0,\ldots,k_1\}$$p(k)\propto \exp\{-bk+bk_0\}$

Pembaruan: Ketika memikirkan lebih lanjut tentang itu, karena berkurang dalam k, dukungan distribusi yang efektif akan selalu dimulai pada . Jika cukup besar, dukungan ini akan berakhir dengan sangat cepat, dalam hal ini tidak menjadi masalah karena nilai tidak akan pernah dikunjungi. Jika sangat kecil, pdfnya hampir datar, yang berarti bahwa seseorang dapat menggunakan distribusi seragam pada sebagai proposal accept-reject. Dan gunakan log pada langkah penerimaan untuk menghindari luapan.$p(\cdot)$$k_0=0$$b$$n$$k$$b$$\{0,1,\ldots,n\}$