最近在看有關蒙特卡洛積分的內容，發現網上很多博主寫的證明過程跳步較為嚴重，而且過程晦澀，不太容易理解。我在自己閱讀國外相關教材附錄后發現證明蒙特卡洛積分方法並不難，利用的僅是概率論的基本知識，現整理下來與大家分享。
那么什么是蒙特卡洛積分？簡而言之就是，在求積分時，如果找不到被積函數的原函數，那么利用經典積分方法是得不到積分結果的，但是蒙特卡洛積分方法告訴我們，利用一個隨機變量對被積函數進行采樣，並將采樣值進行一定的處理，那么當采樣數量很高時，得到的結果可以很好的近似原積分的結果。這樣一來，我們就不用去求原函數的形式，就能求得積分的近似結果。
一、前提知識：
由概率論基本知識，假設一連續型隨機變量$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$倍的計算量才可以達到。即便如此，原始的蒙特卡洛積分方法也不失為是一種經典有效的方法。