共计 5304 个字符，预计需要花费 14 分钟才能阅读完成。
蒙特卡罗积分
引言
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)}$ 计算权重。