如何從具有不可計算 CDF 的分佈中採樣？

半計算機科學模擬相關問題在這裡。

我有一個發行版

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

對於一些常數 b 和 n，並且 x 是一個整數，使得 $ 0\leq x \leq n $ .

現在，我需要從這個分佈中取樣。它有一個可逆的 CDF，所以理論上可以直接做到這一點。問題是涉及的數字很大。實際上如此之大，以至於它們都溢出了傳統格式的變量，並且至少需要幾分鐘（在某些時候我放棄了……）才能使用任意精度格式進行計算。基本上，逆 CDF 仍然涉及一個術語 $ e^{b(n+1)} $ ， 為了 $ 350 < n < 3500 $ . 儘管如此，輸出數字仍將在範圍內 $ 0-n $ ，所以似乎應該有一種方法可以做到這一點。

我正在尋找的是一種從可計算的分佈中近似採樣的方法。是否有其他抽樣方法？這些是什麼？

CDF 很容易反轉。 反演公式導致必須是最簡單和最方便的可能解決方案之一。

首先觀察結果的概率,, 與. 因此，如果我們生成一個統一的值之間和=，我們只需要找到最大的為此

簡單代數給出了解決方案

這是一個與所有其他隨機數生成器一樣構造的R實現：它的第一個參數指定要生成多少個iid值，其餘參數命名參數（作為b和作為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))
}

作為其使用示例，讓我們根據此分佈生成一百萬個隨機變量：

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

(需要幾秒鐘。）

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)

(被添加到每個值以創建更好的直方圖：R的hist過程有一個特質（=錯誤），其中當左端點設置為零時，第一個條形太高。）紅色曲線是此參考分佈模擬嘗試重現。讓我們用卡方檢驗評估擬合優度：

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

p 值為: 漂亮的合身。

引用自：https://stats.stackexchange.com/questions/55533