最近在做一个机器学习的项目,其中用到了泊松随机数。查维基百科 Poisson distribution 发现了一个算法,可以生成泊松随机数:3 ^3 ]8 G! [2 j0 F/ H W, X4 r8 Y
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.3 p. D2 Z1 P5 \2 `: f2 c g& j1 D
词条里面没有解释为什么这个算法可以生成泊松随机数,我在此给出证明。
" i. {; h8 y& {& @( f第一节:算法描述 上面的这个算法可以描述为: 第一步:给定一个参数 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> , 生成一系列的随机数,这些随机数服从 分布,也就是这些随机数在开区间 (0, 1) 之间均匀分布。 第二步:求这些随机数的乘积,当乘积小于或者等于 时,程序停止。记下此时参与求乘积的随机数的个数。 第三步:程序终止时参与乘积的随机数的个数减一服从参数为 的泊松分布。
$ S& P" x" E1 I2 ~1 D- {2 Y2 G第二节:算法的数学表达 为了证明这个算法确实可以生成泊松随机数,我们记
4 r+ j6 k/ @2 F7 j5 |- }: C) |5 _: h' r. L; P; l
% r$ E* \5 I7 w6 @6 I' O z' T: j这就等价于 所以变量 的概率密度为 已知 所以 这是一个指数分布。 因为随机变量 ,所以我们要计算一系列服从指数分布的随机变量的和。已知,对于独立随机变量 ,它们的和 的概率密度为 这是两个概率密度函数的卷积。做傅里叶变换,得到 的概率分布的特征函数是 两个随机变量的概率密度分布的特征函数的乘积。为了计算 的概率分布,我们先要计算指数分布的特征函数。根据特征函数的定义,我们有 所以变量 的概率密度的特征函数为 . 第三节:根据概率密度的特征函数计算所对应的概率密度 现在我们已经知道了概率密度的特征函数,接下来我们要根据这个特征函数计算所对应的概率密度。做傅里叶逆变换可以得到所对应的概率密度分布: 因为变量 是一系列负数的求和,所以上面的积分中, . 选择如下图所示的一个积分围道: ) f7 c. N: e% {' M( o
计算在这个围道上的积分: 这个积分可以分为两部分,第一部分是沿着横轴求积分,第二部分是沿着外面的大圆求积分。可以证明沿着大圆的积分为零,因为 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;"> 时, 以指数速度衰减到零。 所以我们就有 根据柯西积分定理,左边的积分为 上面式子最右边的积分为 所以围道积分为 最终我们得到随机变量 所服从的概率密度函数为 这个分布的名字叫做 分布。显然,根据上式,当 的时候,上面的分布退化为指数分布。
% Q# `7 Q* K8 D/ Z' @# |7 Y第四节:计算随机变量 的概率密度函数 已经知道了 服从 分布,那么计算 的分布也很简单了。已知 所以
/ m4 r: |5 d* S, _6 R& G第五节:计算 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;"> 因为 ,所以我们有
& R* U F4 |$ ]9 q6 K* H) ]第六节:根据对概率的两种等价解释得到泊松分布 现在我们已经算出来了当我们用 个 [0, 1] 均匀分布的随机数连乘时,所得到的乘积小于 0" style="max-width: 100%; vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;"> 的概率。这里,我们相当于是固定了 ,扫描不同的参数 ,得到了概率。我们可以换一个角度。这个概率也可以看作是我们固定了参数 ,计算需要多少个 [0, 1] 之间均匀分布的随机数连乘才能让最后的乘积小于 . 也就是, 根据第五节的结果,我们知道 所以,假设 个 [0, 1] 之间均匀分布的随机数连乘后刚好小于 , 那么 服从这样的概率分布: 这就是泊松分布。
% ^7 K" E3 E. F. z第七节:程序实现 我已经写了一个程序来实现这个算法,并且得到了测试结果。程序GitHub地址为 PrimerLi/Poisson % r/ \6 ^: R7 ?3 s5 r
第八节:实验结果
, @# {1 W; h2 q. k
. n* x' |/ c& Q8 u( y; U
图中显示了 所对应的泊松分布的概率曲线。横轴为 ,纵轴为 . 可以看出,对于不同的参数 ,理论计算出来的结果和用Monte Carlo模拟出来的结果相差不大。
; W3 C2 M7 o: Q |