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

浅夏110

最近在做一个机器学习的项目，其中用到了泊松随机数。查维基百科 Poisson distribution 发现了一个算法，可以生成泊松随机数：
* ^6 j1 y5 k+ t5 halgorithm 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.

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

" I7 O. I& L6 B8 y" M/ X. f

第一节：算法描述

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

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

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

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

" r0 U/ P+ s: [. f+ A

第二节：算法的数学表达

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

% u0 R& ^7 y* f# k  o. V) e4 v# r

这就等价于

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

所以变量 的概率密度为

已知

所以

这是一个指数分布。

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

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

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

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

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

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

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

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

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

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;"> 时， 以指数速度衰减到零。

所以我们就有

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

上面式子最右边的积分为

所以围道积分为

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

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

0 v( m/ U/ N& H  E

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

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

所以

' q2 ~. {( G( B8 W9 _4 L2 t

第五节：计算 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;">

因为 ，所以我们有

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

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

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

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

这就是泊松分布。

) J5 F2 y5 I2 |  W+ e0 U

第七节：程序实现

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

PrimerLi/Poisson

第八节：实验结果

图中显示了 所对应的泊松分布的概率曲线。横轴为 ，纵轴为 . 可以看出，对于不同的参数 ，理论计算出来的结果和用Monte Carlo模拟出来的结果相差不大。
! ?: h0 |7 P+ D

德古拉

非常好, 请问是什么语言的程序?
: _$ x$ |# E  B9 @! x