泊松過程是一個常見的隨機過程，泊松隨機數也是一種常見的隨機數。有時候需要模擬生成泊松隨機數。Knuth曾經提出過這樣的一個演算法用來生成泊松隨機數

［1］

：

algorithm poisson random number （Knuth）：

init：

Let L ← e^{−λ}， k ← 0 and p ← 1。

do：

k ← k + 1。

Generate uniform random number u in ［0，1］ and let p ← p × u。

while p > L。

return k − 1。

維基百科的詞條裡面沒有解釋為什麼這個演算法可以生成泊松隨機數，我在此給出證明。這其實是一個停時問題。

第一節：演算法描述

上面的這個演算法可以描述為：

第一步：給定一個引數

\lambda > 0

， 生成一系列的隨機數，這些隨機數服從

\text{Uniform}(0, 1)

分佈，也就是這些隨機數在開區間 （0， 1） 之間均勻分佈。

第二步：求這些隨機數的乘積，當乘積小於或者等於

e^{-\lambda}

時，程式停止。記下此時參與求乘積的隨機數的個數。這是一個停時特徵。

第三步：程式終止時參與乘積的隨機數的個數減一服從引數為

\lambda

的泊松分佈。

第二節：演算法的數學表達

為了證明這個演算法確實可以生成泊松隨機數，我們記

P = \prod_{i = 1}^{n}X_{i}, X_i \sim \text{Uniform}(0, 1)

這就等價於

\log P = \sum_{i = 1}^{n} \log X_{i}

已知隨機變數

X

的機率密度為

f_X(x)

，令

Y = \log X

。

p(Y \le y) = \int_{-\infty}^{y} f_{Y}(y^{\prime})dy^{\prime} = p(\log X \le y) = p(X \le e^{y}) = \int_{-\infty}^{e^{y}} f_{X}(x)dx

所以變數

Y

的機率密度為

f_{Y}(y) = \frac{d}{dy}\int_{-\infty}^{e^{y}} f_{X}(x)dx = f_{X}(e^{y}) e^{y}

已知

f_{X}(x) =\begin{cases} 1 & 0 \leq x\leq 1 \\ 0 & \text{otherwise} \end{cases}

所以

f_{Y}(y) = \begin{cases} e^{y} & -\infty < y \leq 0 \\ 0 & \text{otherwise} \end{cases}

這是一個指數分佈。

因為隨機變數

\log P = \sum_{i = 1}^{n} \log X_{i} := \sum_{i = 1}^{n} Y_{i}

，所以我們要計算一系列服從指數分佈的隨機變數的和。已知，對於獨立隨機變數

X, Y

，它們的和

Z = X + Y

的機率密度為

f_{Z}(z) = \int_{-\infty}^{\infty} f_{X}(\xi) f_{Y}(z - \xi) d\xi

這是兩個機率密度函式的卷積。做傅立葉變換，得到

Z

的機率分佈的特徵函式是

X, Y

兩個隨機變數的機率密度分佈的特徵函式的乘積。為了計算

\log P

的機率分佈，我們先要計算指數分佈的特徵函式。根據特徵函式的定義，我們有

\hat{f}_{Y}(\eta) = \int_{-\infty}^{\infty} f_{Y}(y) e^{i\eta y} dy = \frac{1}{i\eta + 1}

所以變數

\log P

的機率密度的特徵函式為

\frac{1}{(i\eta + 1)^n}

。

第三節：根據機率密度的特徵函式計算所對應的機率密度

現在我們已經知道了機率密度的特徵函式，接下來我們要根據這個特徵函式計算所對應的機率密度。做傅立葉逆變換可以得到所對應的機率密度分佈：

I(y) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{1}{(i\eta + 1)^n} e^{-i \eta y } d\eta

因為變數

\log P

是一系列負數的求和，所以上面的積分中，

y < 0

。

選擇如下圖所示的一個積分圍道：

泊松隨機數的生成演算法：數學推導和程式實現

計算在這個圍道上的積分：

\frac{1}{2\pi} \oint_{\gamma_R} \frac{1}{(iz + 1)^n}e^{-izy} dz

這個積分可以分為兩部分，第一部分是沿著橫軸求積分，第二部分是沿著外面的大圓求積分。可以證明沿著大圓的積分為零，因為

\Bigg\vert\frac{1}{2\pi} \int_{z = R e^{i\theta}, \sin\theta > 0} \frac{1}{(iz + 1)^n}e^{-iyR(\cos\theta + i\sin\theta)} dz\Bigg\vert \leq \frac{1}{2\pi}\int_{z = R e^{i\theta}} \frac{1}{(R+1)^n} e^{yR\sin\theta} Rd\theta\rightarrow 0

當大圓半徑為無窮大的時候該積分趨近於零，這是因為當

y < 0, \sin\theta > 0

時，

e^{yR\sin\theta}

以指數速度衰減到零。

所以我們就有

\frac{1}{2\pi} \oint_{\gamma_R} \frac{1}{(iz + 1)^n}e^{-izy} dz = \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{1}{(i\eta + 1)^n} e^{-i\eta y} d\eta

根據柯西積分定理，左邊的積分為

\frac{1}{2\pi} \oint_{\gamma_R} \frac{1}{(iz + 1)^n}e^{-izy} dz = \frac{1}{2\pi} \oint_{z = i + \delta e^{i\theta}} \frac{1}{(i\delta e^{i\theta})^n}e^{-iy(i + \delta e^{i\theta})}\delta e^{i\theta}i d\theta = \frac{1}{2\pi} e^{y}\oint\frac{e^{-iy\delta e^{i\theta}}}{(i\delta e^{i\theta})^{n-1}}d\theta

上面式子最右邊的積分為

\oint\frac{e^{-iy\delta e^{i\theta}}}{(i\delta e^{i\theta})^{n-1}}d\theta = \oint \sum_{m = 0}^{\infty} \frac{(-y)^m (i\delta e^{i\theta})^{m-n+1}}{m!}d\theta = \sum_{m = 0}^{\infty}\frac{(-y)^m}{m!} 2\pi \delta_{m, n-1} = 2\pi \frac{(-y)^{n-1}}{(n-1)!}

所以圍道積分為

\frac{1}{2\pi} \oint_{\gamma_R} \frac{1}{(iz + 1)^n}e^{-izy} dz = e^{y}\frac{(-y)^{n-1}}{(n-1)!}, \text{ }y < 0

最終我們得到隨機變數

\log P

所服從的機率密度函式為

f_{\log P}(y) = \begin{cases} \frac{(-y)^{n-1}}{(n-1)!}e^{y} & -\infty < y \leq 0 \\ 0 & \text{otherwise} \end{cases}

這個分佈的名字叫做

\Gamma

分佈。顯然，根據上式，當

n = 1

的時候，上面的分佈退化為指數分佈。

第四節：計算隨機變數

#FormatImgID_42#

的機率密度函式

已經知道了

\log P

服從

\Gamma

分佈，那麼計算

P

的分佈也很簡單了。已知

p(\log P \le y) = \int_{-\infty}^{y} f_{\log P}(y^{\prime}) dy^{\prime} = p(P \le e^{y}) = \int_{-\infty}^{e^{y}} f_{P}(p)dp

所以

f_{P}(p) = p^{-1}f_{\log P}(\log p) = \begin{cases} \frac{(-\log p)^{n-1}}{(n-1)!} & 0 < p \le 1 \\ 0 & \text{otherwise} \end{cases}

第五節：計算

#FormatImgID_48#

的機率

我們已經知道了變數

P

的分佈函式，那麼就可以計算

P < e^{-\lambda}

的機率為

p(P < e^{-\lambda}) = \int_{0}^{e^{-\lambda}} \frac{(-\log p)^{n-1}}{(n-1)!} dp = \frac{1}{(n-1)!} \int_{\lambda}^{\infty} e^{-t} t^{n-1}dt

因為這個機率依賴於

n

，所以可以將這個機率重新寫作

p_n(P < e^{-\lambda} ) = \frac{1}{(n-1)!} \int_{\lambda}^{\infty} e^{-t} t^{n-1}dt

利用分部積分，可以得到機率的遞迴關係為

p_{n}(P < e^{-\lambda}) = \frac{\lambda^{n-1}}{(n-1)!} e^{-\lambda} + p_{n-1}( P < e^{-\lambda}), n > 1

因為

p_{1}(P < e^{-\lambda}) = e^{-\lambda}

，所以我們有

p_{n}(P < e^{-\lambda}) = \sum_{k = 0}^{n-1}\frac{\lambda^k}{k!}e^{-\lambda}

第六節：根據對機率的兩種等價解釋得到泊松分佈

現在我們已經算出來了當我們用

n

個均勻分佈在

(0, 1)

的隨機數

X

連乘時，所得到的乘積小於

e^{-\lambda}, \forall \lambda > 0

的機率。嚴格意義上，這個連乘的結果應該記做

P_{n} = \prod_{i = 1}^{n}X_{i}, X_i \sim \text{Uniform}(0, 1)

前面已經知道這個機率為

p(P_n < e^{-\lambda}) = \sum_{k = 0}^{n-1}\frac{\lambda^k}{k!}e^{-\lambda}, n \ge 1

這是

n

個

X\sim U(0, 1)

連乘後乘積小於

e^{-\lambda}

的機率。記

N

為可以讓

X

連乘後乘積小於

e^{-\lambda}

所需

X

的個數，則有

p(P_{n} < e^{-\lambda}) = p(N \le n)

但是我們要計算的是剛好

n

個

X\sim U(0, 1)

連乘後的乘積首次小於

e^{-\lambda}

的機率，這個機率為

\begin{align} p(N = n) & = p(N \le n) - p(N \le n-1)\\ &= p(P_{n} < e^{-\lambda}) - p(P_{n-1} < e^{-\lambda}) \\ &= \frac{\lambda^{n-1}}{(n-1)!} e^{-\lambda} \end{align}

這就是泊松分佈。

第七節：程式實現

我已經寫了一個程式來實現這個演算法，並且得到了測試結果。程式GitHub地址為

PrimerLi/Poisson

第八節：實驗結果

泊松隨機數的生成演算法：數學推導和程式實現

圖中顯示了

\lambda = 1, \lambda = 4, \lambda = 10

所對應的泊松分佈的機率曲線。橫軸為

k

，縱軸為

p(k) = \frac{\lambda^k}{k!} e^{-\lambda}

。 對於每一個

\lambda

，我都畫了兩條線，一條線是解析解，另一條是蒙特卡洛模擬的結果。可以看出，對於不同的引數

\lambda

，理論計算出來的結果和用Monte Carlo模擬出來的結果相差不大。

參考

^

Generating Poisson-distributed random variables

https：//en。wikipedia。org/wiki/Poisson_distribution