使用蒙特卡洛模拟计算定积分python 用蒙特卡洛求积分

最近在看有关蒙特卡洛积分的内容，发现网上很多博主写的证明过程跳步较为严重，而且过程晦涩，不太容易理解。我在自己阅读国外相关教材附录后发现证明蒙特卡洛积分方法并不难，利用的仅是概率论的基本知识，现整理下来与大家分享。

那么什么是蒙特卡洛积分？简而言之就是，在求积分时，如果找不到被积函数的原函数，那么利用经典积分方法是得不到积分结果的，但是蒙特卡洛积分方法告诉我们，利用一个随机变量对被积函数进行采样，并将采样值进行一定的处理，那么当采样数量很高时，得到的结果可以很好的近似原积分的结果。这样一来，我们就不用去求原函数的形式，就能求得积分的近似结果。

一、前提知识：

由概率论基本知识，假设一连续型随机变量$X$的样本空间为$D$，其概率密度分布函数为$p(x)$，则其数学期望为

$\begin{equation} E\left[X\right]=\int_D {xp(x){\rm{d}}x} \end{equation}$若另一连续随机变量Y满足$Y=f(X)$，则$Y$的数学期望$E[Y]$可由下式给出
$\begin{equation} E\left[Y\right]=\int_D {f(x)p(x){\rm{d}}x} \end{equation}$二、蒙特卡洛积分与重要性采样

根据以上叙述，假设这里我们要计算一个一维积分式
$\begin{equation} A=\int_a^b {f(x){\rm{d}}x} \end{equation}$根据经典的方法，我们需要求得$f(x)$的原函数$F(x)$，才能解出这个积分结果，但如果$f(x)$的原函数形式复杂，或者根本求不出来，总之在不知道$F(x)$的具体形式的情况下，如果我们还想计算这个积分，怎么办？这时候我们就需要借助蒙特卡洛积分(Monte Carlo Integration)方法。蒙特卡洛积分方法告诉我们，为求得积分结果，可以构造
$\begin{equation} \label{eq:FN} F_N = \frac{{b - a}}{N}\sum\limits_{i = 1}^N {f(X_i )} \end{equation}$其中的每一个$X_i(i=1,2,3,...,N)$为$[a,b]$之间的均匀连续随机变量，所有的$X_i$组成一个随机变量集合。下面证明，$F_N$的数学期望即为我们要求的结果$A$。
$\begin{equation} \begin{split} \label{eq:proof1} E\left[ {F_N } \right] &= E\left[ {\frac{{b - a}}{N}\sum\limits_{i = 1}^N {f(X_i )} } \right] \\ &=\frac{{b - a}}{N}\sum\limits_{i = 1}^N {E\left[ {f(X_i )} \right]} \\ &= \frac{{b - a}}{N}\sum\limits_{i = 1}^N {\int_a^b {f(x)p(x){\rm{d}}x} } \\ &= \frac{{b - a}}{N}\sum\limits_{i = 1}^N {\int_a^b {f(x)\frac{1}{{b - a}}{\rm{d}}x} } \\ &= \frac{{b - a}}{N}\frac{1}{{b - a}}\sum\limits_{i = 1}^N {\int_a^b {f(x){\rm{d}}x} } \\ &= \frac{1}{N}\sum\limits_{i = 1}^N {\int_a^b {f(x){\rm{d}}x} } \\ &= \int_a^b {f(x){\rm{d}}x} \end{split} \end{equation}$以上证明过程表明，若我们根据式(\ref{eq:FN})来构造一个新的随机变量$F_N$，则$F_N$的期望就是积分$$\int_a^b{f(x){\rm{d}}x}$$的结果，随着$N$的增加，$F_N$就越逼近理论上$A$的值，即$F_N$是$A$的一个无偏估计。这样我们实际上就避开了求$f(x)$的原函数$F(x)$的过程。整个求积分近似值的过程则可以用文字描述如下：首先从区间$[a,b]$上对均匀分布的随机变量$X$连续取样$N$次，得到$N$个取样值$\{x_1,x_2,x_3,...,x_N\}$，对每个取样值$x_i(i=1,2,3,...,N)$计算$f(x_i)$得到$\{f(x_1),f(x_2),f(x_3),...,f(x_N)\}$，再计算它们的和$$\sum\limits_{i=1}^{N}{f(x_i)}$$，最后乘系数$\frac{b-a}{N}$即可得到对理论值的一个估计。

进一步对上述过程分析，我们发现这里的$X$被规定为与原积分区间相同的均匀分布随机变量。那么对于与原积分区间相同，但却不是均匀分布的一般随机变量，上述结论是否仍成立？结论是，如果这里的随机变量$X$的概率密度分布函数已知，那么上述结论还是成立的。假设其概率密度分布函数为$p(x)$，证明如下

类似式(\ref{eq:FN})，这次我们构造
$\begin{equation} \label{eq:FN2} F_N = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{f(X_i )}}{{p(X_i )}}} \end{equation}$再构造随机变量
$$Y=\frac{f(X)}{p(X)}$$
式(\ref{eq:FN2})的所有量都是已知的。则$F_N$的数学期望$E\left[F_N\right]$即为

$\begin{equation} \begin{split} \label{eq:proof2} E\left[ {F_N } \right] &= E\left[ {\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{f(X_i )}}{{p(X_i )}}} } \right] \\ &= \frac{1}{N}\sum\limits_{i = 1}^N {E\left[ {\frac{{f(X_i )}}{{p(X_i )}}} \right]} \\ &= \frac{1}{N}\sum\limits_{i = 1}^N {E\left[ {Y_i } \right]} \\ &= \frac{1}{N}\sum\limits_{i = 1}^N {\int_a^b {\left[ {\frac{{f(x)}}{{p(x)}}} \right]p(x){\rm{d}}x} } \\ &= \frac{1}{N}\sum\limits_{i = 1}^N {\int_a^b {f(x){\rm{d}}x} } \\ &= \int_a^b {f(x){\rm{d}}x} \end{split} \end{equation}$我们发现式(5)其实是式(7)的特殊形式。与之前的要求不同，这里仅要求随机变量$X$的概率密度分布函数$p(x)$已知且在$X$的样本空间内$p(x)\ne0$。综合上述叙述，我们可以得到蒙特卡洛积分方法如下
$\begin{equation} \int_D {f(x){\rm{d}}x} = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{f(X_i )}}{{p(X_i )}}} \end{equation}$我们可以进一步地将求解一维积分的方法扩展到求解高维积分中去。考虑求解如下积分

$\begin{equation} \int_{z_1 }^{z_2 } {\int_{y_1 }^{y_2 } {\int_{x_1 }^{x_2 } {f(x,y,z){\rm{d}}x{\rm{d}}y{\rm{d}}z} } } \end{equation}$只需要在积分域$$\int_{z_1 }^{z_2 } {\int_{y_1 }^{y_2 } {\int_{x_1 }^{x_2 } } } $$定义的盒型区域内选取概率密度分布为
$$$ p(x) = \frac{1}{{\left( {x_2 - x_1 } \right)\left( {y_2 - y_1 } \right)\left( {z_2 - z_1 } \right)}} $$$的均匀随机变量$X$，则积分结果为
$$$ \frac{{\left( {x_2 - x_1 } \right)\left( {y_2 - y_1 } \right)\left( {z_2 - z_1 } \right)}}{N}\sum\limits_{i = 1}^N {f\left( {X_i } \right)} $$$现在我们分析一下随着样本数量$N$的增加，估计值$F_N$的方差$$\sigma ^2 \left[ {F_N } \right]$$的变化情况，以便得出蒙特卡洛积分方法的收敛速度特性。在式(7)的基础上，我们继续计算$$\sigma ^2 \left[ {F_N } \right]$$如下

$\begin{equation} \begin{split} \label{eq:proof3} \sigma ^2 \left[ {F_N } \right] &= \sigma ^2 \left[ {\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{f(X_i )}}{{p(X_i )}}} } \right] \\ &= \frac{1}{{N^2 }}\sum\limits_{i = 1}^N {\sigma ^2 \left[ {\frac{{f(X_i )}}{{p(X_i )}}} \right]} \\ &= \frac{1}{{N^2 }}\sum\limits_{i = 1}^N {\sigma ^2 \left[ {Y_i } \right]} \\ &= \frac{1}{{N^2 }}\left( {N\sigma ^2 \left[ Y \right]} \right) \\ &= \frac{1}{N}\sigma ^2 \left[ Y \right] \end{split} \end{equation}$
所以
$\begin{equation} \sigma \left[ {F_N } \right] = \frac{1}{{\sqrt N }}\sigma \left[ Y \right] \end{equation}$上式告诉我们，估计值的不稳定来源于随机变量$Y$的取值不稳定。换句话说，如果$$\frac{{f(X_i )}}{{p(X_i )}}$$因不同$X_i$的取值变化地越剧烈，就会造成$Y$的方差较大，则会造成估计值的收敛速度越慢。这启示我们，若$p(x)$的形状越接近$f(x)$，则有益于最终结果的收敛。上述思想即为“重要性采样”方法，即对积分值有重要贡献($f(x)$较大)的被积函数区间，我们以较大概率生成处于这个区间附近的随机变量，用于快速逼近理论值。这也就是为什么我们要引入任意分布随机变量的蒙特卡洛积分方法，而不满足于利用均匀分布随机变量来求蒙特卡洛积分的原因。

利用蒙特卡洛方法，我们可以得到任意一个积分的结果，但是代价就是我们得不到其理论值，我们得到的只是一个对理论值的估计，估计值与理论值之间的误差可以通过增加样本数来减小，但收敛速率仅为$$O\left( {\sqrt N } \right)$$，也就是说，若想将误差降为现在的一半，我们需要再多计算$4$倍的计算量才可以达到。即便如此，原始的蒙特卡洛积分方法也不失为是一种经典有效的方法。