[matlab]一种生成poisson随机数的算法

一种生成poisson随机数的算法

背景知识——naive monte carlo

我们知道，利用naive monte carlo 来求poisson 随机变量的期望可以表示为如下公式

$$\hat{X}_{n}=\frac{1}{n}\sum^{n}_{i=1}X_{i}$$
其中 $\{X_{i}\}_{i=1,\dots,n}$独立同分布，服从参数为 $\beta$的poisson分布。

由大数定律可知，有

$$\hat{X}_{n} \longrightarrow^{a.s.} \mathbb{E}[X_i]=\beta$$
由中心极限定理知，上述逼近在弱收敛意义下的收敛速度是 $\mathcal{O}(\frac{1}{\sqrt{n}})$。
因此，为了计算 $\hat{X}_n$，自然而然的一个问题是如何生成 $X_i$。

经典poisson随机变量生成

最经典，也是最常用的方法就是利用poisson分布函数的逆函数来产生服从poisson 分布的随机数，步骤如下：

• 从（0，1）上的均匀分布采样出u
• 利用poisson 分布的逆函数得到poisson 变量
  $$m^{*}:=\inf{\{m| F_{\beta}(m)>u\}}$$
  其中$F_{\beta}$是poisson分布函数
  $$F_{\beta}(m)= \sum^{m}_{k=0}\frac{\beta^{k}e^{-\beta}}{k!}\quad m\in \mathbb{Z}^{+}$$

但这种方法在$\beta$很大时，计算量会比较大。

Proposed poisson variable generator

该新方法基于的理论支持是:poisson分布$\mathcal{P}( \beta)$分布一致逼近正态分布$\mathcal{N}(\beta, \beta)$ 。参考文献(1)
该方法步骤如下：

• 第一步：生成服从$\mathcal{N}(0,1)$分布的随机数$z$ (参考另外一篇生成正态随机数的博客)
• 第二步：计算$u=\Phi(z)$, 则$u$服从(0,1)上的均匀分布
• 第三步：按如下步骤计算poisson分布函数的逆函数$F^{-1}_{\beta}$寻找$m^*$
  
  • $m_{0}=\max([\beta+z\sqrt{\beta}], 0)$
  • 如果 $F_{\beta}(m_{0})<u$, 则将使$m_{0}=1+m_{0}$直到$F_{\beta}(m_{0})>u$,并记$m^*=m_{0}$
  • 如果 $F_{\beta}(m_{0})>u$, 则将使$m_{0}=m_{0}-1$直到$F_{\beta}(m_{0})<u$,并记$m^*=m_{0}+1$

代码示例

function x=poisson(lamda, size)
%输入：poisson分布的强度lamda, 要产生的随机变量数目size
%输出：poisson 随机变量序列x
x=zeros(1,2*size);
for i=1:size
    u=unifrnd(0,1,1,2); %生成(0,1)上的均匀分布变量
   %% 生成独立的服从参数为（0，1）的正态分布变量(参考另外一篇博客中的介绍)
    X1=sqrt(-2*log(u(1)))*cos(2*pi*u(2));
    X2=sqrt(-2*log(u(1)))*sin(2*pi*u(2));
   %% 生成独立的服从参数为 lamda 的poisson分布变量
   z=normcdf([X1, X2]);
   x(2*i-1)=pois_rand(lamda, X1, z(1)); x(2*i)=pois_rand(lamda, X2, z(2));
end
end 

function m=pois_rand(lamda, x, z)
%搜索m*
    m0=max(floor(lamda+x*sqrt(lamda)), 0);
    if F(lamda, m0)<z
        m=m0+1;
        while F(lamda, m)<z
            m=m+1;
        end
    else 
        m=m0;
        while F(lamda, m)>=z
            m=m-1;
        end
        m=m+1;
    end
end

function F=F(lamda, m)
%poisson分布的分布函数
if m<0
    F=0;
else
     F=exp(-lamda);
     for i=1:m
         F=F+lamda^(i)*exp(-lamda)./factorial(i);
     end
end
end

运行结果：

>> x=poisson(2,1e5);[mean(x), var(x)]

ans =

    1.9993    2.0010
>> hist(x,100)

注：
该方法更多细节参考文献(2)
————————————————————————————————————