《Introducing Monte Carlo Methods with R》第三章笔记

蒙特卡罗积分

引言

motivation

数值方法求积分在统计计算中具有十分重要的地位，很多需要求边缘概率密度分布的情况（尤其是算贝叶斯后验）需要用到数值积分。

R中普通一维数值积分方法

常用的数值积分方法在不管是本科还是研究生的数值分析中都进行了详细的介绍，包括基于梯形公式，辛普森公式，基于Romberg积分等方法，这篇文章中暂时不进行详细介绍。值得一提的是，R中实现了两种一维的数值积分法：area（from MASS）和integrate。

实验1 利用R内置函数计算一维数值积分

待求的积分如下：

$$ f(x)=\int\_{0}^{\infty}x^{\lambda-1}exp(-x)dx $$

采用integrate函数进行数值计算并和利用gamma函数求得的解析解进行对比，代码如下：

> ch=function(la){integrate(function(x){x^(la-1)*exp(-x)},0,Inf)$val}
> plot(lgamma(seq(0.01,10,le=50)),log(apply(as.matrix(seq(0.01,10,le=50)),1,ch)),xlab = "log(integrate(f))",ylab = expression(log(Gamma(lambda))))

得到的结果如下图所示，可以看出在$\lambda \in [0.01,10]$的范围内integrate函数求得数值积分和解析解并无明显差异。

第二个要进行数值积分的函数如下式所示，这个函数可以简单认为是当先验为均匀分布时10个柯西分布样本对模型参数的边缘后验（未进行归一化）。

$$ m(x)=\int\_{-\infty}^{\infty}\sum^{10}*{i=1} \frac{1}{\pi} \frac{1}{1+(x*{i}-\theta)} dx $$

设$\theta =350$，代码如下，从运行结论来看显然integrate函数算的有错，因为把积分上限从正无穷改为400后结果竟然变小了。

> cac=rcauchy(10)+350
> like=function(theta){
+     u=dcauchy(cac[1]-theta)
+     for(i in 2:10)
+         u=u*dcauchy(cac[i]-theta)
+     return(u)
+ }
> integrate(like,-Inf,Inf)
5.559304e-44 with absolute error < 1.1e-43
> integrate(like,-Inf,400)
5.882426e-31 with absolute error < 1.1e-30

重新生成10个样本，这次设置$\theta =0$，然后将integrate函数的结果和area函数的结果对比，代码如下。

> library(MASS)
> cac=rcauchy(10)
> nin=function(a){integrate(like,-a,a)$val}
> nan=function(a){area(like,-a,a)}
> x=seq(1,1000,le=10000)
> y=log(apply(as.matrix(x),1,nin))
> z=log(apply(as.matrix(x),1,nan))
> plot(x,y,ylim = range(cbind(y,z)),lwd=2)
> lines(x,z,col='sienna',lwd=6)

实验结果如图所示，图中黑色线为积分限-integrate函数积分值的曲线，很显然错误，黄色线为积分限-area函数积分值的曲线，其趋势基本正确。

经典蒙特卡罗积分

从上文的实验可以看出，当积分的形式变得比较复杂时R内建的一维数值积分算法可能会不稳定。这种常规数值积分方法在计算高维积分时也面临效率和精度问题。因此，面对统计推断中随处可见的高维复杂积分，蒙特卡罗方法被用于解决这个问题。
利用蒙特卡罗法进行数值积分的原理可以写成如下公式，其中$\mathbb{E}_{f}[h(X)]$表示蒙特卡罗生成的$X \sim f(x)$的大量样本的$h(x)$均值的理论表达式，将样本均值$\hbar_{n}$作为$\mathbb{E}_{f}[h(x)]$的估计量，就可以间接计算出在$ f(x)$支撑集上对$h(x)f(x)$进行积分的结果。

$$ \begin{aligned} \mathbb{E}_{f}[h(X)] &=\int_{\mathcal{X}}^{}h(x)f(x)dx \hbar_{n} &=\frac{1}{n}\sum^{n}_{j=1}h(x_{j}) \end{aligned} $$

在应用场景中，最简单直接的方法就是将被积函数作为$h(x)$，将$f(x)$设为积分区间上的均匀分布$U(a,b)$，生成大量的$X \sim U(a,b)$样本点并用$\hbar_{n}$作为积分结果的估计值。

实验2 利用经典蒙特卡罗方法计算数值积分

待求的积分如下：

$$ h(x)=[cos(50x)+sin(20x)]^2 \qquad x\in\ U(0,1) $$

R代码如下，从图中和integrate函数结果来看，蒙特卡罗积分的结果和integrate结果近似，在生成约4000个样本后可以认为$\hbar_{n}$收敛于$\mathbb{E}_{f}[h(X)]$即积分结果。

> h=function(x){(cos(50*x)+sin(20*x))^2}
> par(mar=c(2,2,2,1),mfrow=c(2,1))
> curve(h,xlab = "function",ylab = "",lwd=2)
> integrate(h,0,1)
0.9652009 with absolute error < 1.9e-10
> x=h(runif(10000))
> estint=cumsum(x)/(1:10000)
> esterr=sqrt(cumsum((x-estint)^2))/(1:10000)
> plot(estint,xlab = "mean and error range",lwd=2,ylim = mean(x)+20*c(-esterr[10000],esterr[10000]),ylab = "")
> lines(estint+2*esterr,col="gold",lwd=2)
> lines(estint-2*esterr,col="gold",lwd=2)

重要性采样（importance sampling）

这种方法是经典蒙特卡罗方法的一种改进，通过给抽样赋予不同权重的方法提高算法的性能和效率。

利用重要性采样求积分

重要性采样求积分的原理可以用下式描述：

$$ \begin{aligned} \mathbb{E}_{f}[h(X)] &=\int_{\mathcal{X}}^{}h(x)f(x)dx &=\int_{\mathcal{X}}^{}h(x)\frac{f(x)}{g(x)}g(x)dx &=\mathbb{E}_{g}[h(X)\frac{f(X)}{g(X)}] \end{aligned} $$

经过变换后，生成$X \sim f(x)$的样本点并计算$h(x)$的样本均值（经典蒙特卡罗）转变为生成$X \sim g(x)$的样本点并计算$h(x)\frac{f(x)}{g(x)}$的样本均值。当知道$f(x)$表达式且$X \sim f(x)$较难生成时可以采用重要性采样方法，用一款R内建分布作为$g(x)$，计算$h(x)\frac{f(x)}{g(x)}$的样本均值实现原积分的近似计算。

实验3 利用重要性采样求积分

对于$Z \sim N(0,1)$，求$P(Z>4.5)$。对于这种求长尾分布的尾部的概率密度函数积分问题，采用传统的蒙特卡罗方法效率极低，因为如果采用接受-拒绝采样的话接受率极低无比。因此，可以构造一个支撑集和这个长尾分布的尾部差不多的提议分布$g(x)$，然后通过采样提议分布并求权重和样本均值的方法进行数值逼近。这样做首先可以通过选择容易生成的提议分布的方法加速，其次每一个生成的样本都对结果有贡献，比带接受率的方法更高效，减少了算力浪费。
利用如下一个位移过的指数分布作为提议分布：

$$ g(x)=e^{-y}/\int^{\infty}_{4.5}e^{-x}dx=e^{-(y-4.5)} $$

这个重要性抽样中$h(X)=1$，权重$\frac{f(X)}{g(X)}$计算如下：

$$ \frac{1}{n}\sum_{i=1}^{n}\frac{f(Y_{i})}{g(Y_{i})}=\frac{1}{n}\sum_{i=1}^{n}\frac{e^{-Y_{i}^{2}+Y_{i}-4.5}}{\sqrt{2\pi}} $$

R代码如下：

> N=10000
> y=rexp(N)+4.5
> w=dnorm(y)/dexp(y-4.5)
> plot(cumsum(w)/1:N)
> abline(a=pnorm(-4.5),b=0,col='red',lwd=2)

从图中可以看出采样约4000个以后结果收敛到真实值（约3.5e-6），如果采用接受-拒绝采样，拒绝率会极高无比。

利用重要性采样模拟分布——采样重要性重采样（sampling importance resampling, SIR）

重要性采样技巧也可以用于生成给定分布的样本，这时候需要对生成的提议分布的样本点进行有放回的重采样，每个粒子的被抽中概率正比于重要性采样的权重。这种方法叫做sampling importance resampling，SIR。

“采样重要性重采样”这个翻译是我从西安交通大学出版社的中文版Givens《计算统计》中拿来的，感觉很奇怪，人民邮电出版社翻译的这本书也这么写，既然这样就直接采用英文名字和缩写方便一点。

在重要性采样选择提议分布$g(x)$时就要求其支撑集要包含待仿真分布$f(x)$的支撑集，这说明所有$X\sim f(x)$的样本都有可能出现在$Y\sim g(y)$中，只是Y的样本出现在X中的概率不同而已，这种概率的差异可以用权重$\frac{f(X)}{g(Y)}$表示。因此，利用生成的提议分布$g(x)$的样本重采样生成$f(x)$的样本是完全可行的。严格一些的论证方法如下：

$$ \begin{aligned} P(X^{*} \in A) &=\sum_{i=1}^{n}P(X^{} \in A \quad and \quad X^{*}=X_{i}) \\ &=\int_{A}^{}\frac{f(x)}{g(x)}g(x)dx\\ &=\int_{A}^{}f(x)dx\\ &=F(x) \end{aligned} $$

其中$X_{i}$为$g(x)$的样本点，$X^{*}$为重采样后的$f(x)$样本点，A代表重采样中被抽中的点，每个样本点的重采样抽中概率$P(X^{*} \in A| X^{*}=X_{i}) \propto \frac{f(X_{i})}{ng(X_{i})}$。

重要性函数选择

要使重要性采样求积分的的结果收敛（方差不能无限大），需要满足下式：

$$ \begin{aligned} \mathbb{E}_{g}[h^{2}(X)\frac{f^{2}(X)}{g^{2}(X)}] &=\mathbb{E}_{f}[h^{2}(X)\frac{f^{}(X)}{g^{}(X)}]\\ &=\int_{\mathcal{X}}^{}h^{2}(x)\frac{f^{2}(x)}{g(x)}dx<\infty \end{aligned} $$

如果上式不满足导致方差无限大，模拟中的积分估计结果不仅不会随着样本数的增大而收敛至真值，反而会出现多次的跳变。

实验4 方差无限大情况下的重要性采样

利用正态分布$N(0,1)$作为提议分布生成柯西分布$Cauchy(0,1)$并计算其在[2,6]区间上的积分，代码如下：

> x=rnorm(100000)
> wein=dcauchy(x)/dnorm(x)
> plot(cumsum(wein*(x>2)*(x<6))/cumsum(wein))
> abline(a=pcauchy(6)-pcauchy(2),b=0,col="sienna",lwd=2)

图中红色黑色线为计算出的积分结果，红色横线为真实值。虽然经过在生成了大量样本后重要性采样的结果接近真实值，但是其中清晰可见发生了一次剧烈的跳变，而且这种在理论上这种方差无限的计算结果不具有理论上的大样本收敛性，不能保证计算精度。

一种改进的重要性采样算法

一种改善收敛性并且可以用于提高模拟效率的方法（被称为“defensive sampling”）是用如下的混合提议分布$g_{2}(x)$取代原有的提议分布$g(x)$：

$$ g_{2}(x)=\rho g(x)+(1-\rho)l(x) $$

其中$l(x)$是比$g(x)$在收敛性方面表现更好的另一种提议分布。之后就从$g_{2}(2)$中抽样，并用公式$\frac{f(X)}{g_{2}(Y)}$计算权重。