泊松随机数的生成算法：数学推导和程序实现

浅夏110

最近在做一个机器学习的项目，其中用到了泊松随机数。查维基百科 Poisson distribution 发现了一个算法，可以生成泊松随机数：
4 E- v$ N: R) X2 e7 b3 r) d; zalgorithm 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.
6 \5 D: y, v; ]7 ^

词条里面没有解释为什么这个算法可以生成泊松随机数，我在此给出证明。

7 @% ~. @( \- f' B7 a" M; t

第一节：算法描述

上面的这个算法可以描述为：

第一步：给定一个参数 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> , 生成一系列的随机数，这些随机数服从 分布，也就是这些随机数在开区间 (0, 1) 之间均匀分布。

第二步：求这些随机数的乘积，当乘积小于或者等于 时，程序停止。记下此时参与求乘积的随机数的个数。

第三步：程序终止时参与乘积的随机数的个数减一服从参数为 的泊松分布。

- \0 i, Q6 O* y5 I: ^

第二节：算法的数学表达

为了证明这个算法确实可以生成泊松随机数，我们记

+ j& m8 W5 ^5 d. W

这就等价于

已知随机变量 的概率密度为 ，令 .

所以变量 的概率密度为

已知

所以

这是一个指数分布。

因为随机变量 ，所以我们要计算一系列服从指数分布的随机变量的和。已知，对于独立随机变量 ，它们的和 的概率密度为

这是两个概率密度函数的卷积。做傅里叶变换，得到 的概率分布的特征函数是 两个随机变量的概率密度分布的特征函数的乘积。为了计算 的概率分布，我们先要计算指数分布的特征函数。根据特征函数的定义，我们有

所以变量 的概率密度的特征函数为 .

第三节：根据概率密度的特征函数计算所对应的概率密度

现在我们已经知道了概率密度的特征函数，接下来我们要根据这个特征函数计算所对应的概率密度。做傅里叶逆变换可以得到所对应的概率密度分布:

因为变量 是一系列负数的求和，所以上面的积分中， .

选择如下图所示的一个积分围道：

计算在这个围道上的积分：

这个积分可以分为两部分，第一部分是沿着横轴求积分，第二部分是沿着外面的大圆求积分。可以证明沿着大圆的积分为零，因为

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" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;">

当大圆半径为无穷大的时候该积分趋近于零，因为当 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> 时， 以指数速度衰减到零。

所以我们就有

根据柯西积分定理，左边的积分为

上面式子最右边的积分为

所以围道积分为

最终我们得到随机变量 所服从的概率密度函数为

这个分布的名字叫做 分布。显然，根据上式，当 的时候，上面的分布退化为指数分布。

8 v9 T% ~, f8 ?  i. ]+ N

第四节：计算随机变量 的概率密度函数

已经知道了 服从 分布，那么计算 的分布也很简单了。已知

所以

第五节：计算 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> 的概率

我们已经知道了变量 的分布函数，那么就可以计算 的概率为

因为这个概率依赖于 ，所以可以将这个概率重新写作

利用分部积分，可以得到概率的递归关系为

1" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;">

因为 ，所以我们有

' t# c7 ^/ o, H( x

第六节：根据对概率的两种等价解释得到泊松分布

现在我们已经算出来了当我们用 个 [0, 1] 均匀分布的随机数连乘时，所得到的乘积小于 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> 的概率。这里，我们相当于是固定了 ，扫描不同的参数 ，得到了概率。我们可以换一个角度。这个概率也可以看作是我们固定了参数 ，计算需要多少个 [0, 1] 之间均匀分布的随机数连乘才能让最后的乘积小于 . 也就是，

根据第五节的结果，我们知道

所以，假设 个 [0, 1] 之间均匀分布的随机数连乘后刚好小于 ， 那么 服从这样的概率分布：

这就是泊松分布。

, Z/ n. `; A& k" I

第七节：程序实现

我已经写了一个程序来实现这个算法，并且得到了测试结果。程序GitHub地址为

PrimerLi/Poisson

$ E8 h* ]% N8 `$ ^5 }

第八节：实验结果

' G9 T0 [2 c( {- |$ n1 n
; ?/ k( _+ A9 f2 p) A图中显示了 所对应的泊松分布的概率曲线。横轴为 ，纵轴为 . 可以看出，对于不同的参数 ，理论计算出来的结果和用Monte Carlo模拟出来的结果相差不大。

德古拉

非常好, 请问是什么语言的程序?
6 e( u2 K# U7 }  ]+ |