[matlab]⼀种⽣成poisson随机数的算法
⼀种⽣成poisson随机数的算法
背景知识——naive monte carlo
我们知道，利⽤naive monte carlo 来求poisson 随机变量的期望可以表⽰为如下公式
其中 独⽴同分布，服从参数为 的poisson分布。
由⼤数定律可知，有
由中⼼极限定理知，上述逼近在弱收敛意义下的收敛速度是 。
因此，为了计算 ，⾃然⽽然的⼀个问题是如何⽣成 。
经典poisson随机变量⽣成
最经典，也是最常⽤的⽅法就是利⽤poisson分布函数的逆函数来产⽣服从poisson 分布的随机数，步骤如下：从（0，1）上的均匀分布采样出u
利⽤poisson 分布的逆函数得到poisson 变量
其中是poisson分布函数
但这种⽅法在很⼤时，计算量会⽐较⼤。
Proposed poisson variable generator
该新⽅法基于的理论⽀持是:poisson分布分布⼀致逼近正态分布 。参考⽂献()
该⽅法步骤如下：
第⼀步：⽣成服从分布的随机数 (参考另外⼀篇⽣成正态随机数的博客)
第⼆步：计算, 则服从(0,1)上的均匀分布
第三步：按如下步骤计算poisson分布函数的逆函数寻
如果 , 则将使直到，并记
如果 , 则将使直到，并记
代码⽰例
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 =
matlab生成随机数
1.9993
2.0010
>> hist(x,100)
注：
该⽅法更多细节参考⽂献() ————————————————————————————————————