[matlab]Variance reduction for poisson jump process (Antithetical method)

针对poisson 过程的一种Antithetical方差减小方法

背景介绍

在上一篇博客中，我们介绍了Naive monte carlo 方法生成Poisson random variables。但是我们知道Naive monte carlo 方法的收敛速度只有$\mathcal{O}(\frac{1}{\sqrt{n}})$, 可以说是相当的慢了。因此，为了提升收敛速度，variance reduction 方法应运而生。variance reduction 方法可以很大的降低采样的方差，从而仅通过少量样本便可以达到很高的精度，大大提升了mento carlo 方法的效率。有关variance reduction 的介绍和原理参考维基百科(1), 在此处，我们仅介绍一种Antithetical variance reduction 方法及其实现。

Antithetical variance reduction method

（1）我们知道poisson jump variable 的naive monte carlo采样步骤为：

• 生成（0，1）上均匀分布的随机变量 $$v \sim^{i.i.d} \text{Unif}(0,1)$$
• 利用poisson分布$F_{\beta}(m)$的逆函数$F^{-1}_{\beta}$ 产生poisson随机变量 $$m^{*}:=F^{-1}_{\beta} (v)$$

（2）利用Antithetical variance reduction 方法的话，采样步骤如下

• 生成（0，1）上均匀分布的随机变量 $$v \sim^{i.i.d} \text{Unif}(0,1)$$
• 利用poisson分布$F_{\beta}(m)$的逆函数$F^{-1}_{\beta}$ 产生随机变量
  $$m_{1}:=F^{-1}_{\beta} (u)$$ $$m_{2}:=F^{-1}_{\beta} (1-u)$$
• 产生最终的poisson 分布随机变量
  $$m^{*}=\frac{m_1+m_2}{2}$$

算法步骤梗概

1. 利用上一篇博客中的方法产生服从$\mathcal{N}(0,1)$分布的随机变量$z$
2. 使$u=\Phi(z)$则自然满足 $u\sim \text{Unif}(0,1)$
3. 利用上一篇博客中的方法（或者自己知道的高效算法）寻找 $$m_{1}:=F^{-1}_{\beta} (v)$$ $$m_{2}:=F^{-1}_{\beta} (1-v)$$
4. 得到poisson随机变量 $$m^{*}=\frac{m_1+m_2}{2}$$

实验示例

antithetical monte carlo 算法代码

function x=var_red_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]);
   %% variance reduction by antithetical mehtod 
   m0=max(floor(lamda+X1*sqrt(lamda)), 0); m1=max(floor(lamda-X1*sqrt(lamda)), 0);
   x(2*i-1)=(pois_rand(lamda, m0, z(1))+pois_rand(lamda, m1, 1-z(1)))./2;
   m0=max(floor(lamda+X2*sqrt(lamda)), 0); m1=max(floor(lamda-X2*sqrt(lamda)), 0);
   x(2*i)=(pois_rand(lamda, m0, z(2))+pois_rand(lamda, m1, 1-z(2)))./2;
end
end 

function m=pois_rand(lamda, m0, z)
%搜索m*
    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

在这里，我们比较naive monte carlo 和 antithetical monte carlo 方法产生poisson随机变量的效率

sizes=[1e1, 1e2, 1e3, 1e4, 1e5]; K=length(sizes);
means=zeros(2, K); vars=zeros(2, K);
lamda=2;
for i=1:K
    %naive monte carlo
    x=poisson(lamda, sizes(i));
    means(1,i)=mean(x); vars(1,i)=var(x);
    %antithetical monte carlo
    x=var_red_poisson(lamda, sizes(i));
    means(2,i)=mean(x); vars(2,i)=var(x);
end
means
vars

得到的结果如下

means =
    %naive monte carlo 算得的lamda 
    1.6500    1.9750    2.0220    2.0076    2.0038
    %antithetical monte carlo 算的的lamda
    2.1000    2.0075    1.9970    2.0023    2.0001
vars =
    %naive monte carlo 方法的方差
    1.5026    2.1451    2.0675    1.9923    2.0068
    %antithetical monte carlo 方法的方差
    0.2789    0.0866    0.1090    0.1160    0.1135

显然，我们可以看到antithetical monte carlo 收敛速度明显比naive monte carlo 快很多。
注： 更详细的细节和理论分析参考文献（2）

---