蒙特卡罗方法（Monte Carlo Method，MC），是一类基于概率与统计理论的数值计算方法，常用于常规方法难以解决的数学与物理问题。广泛应用于计算物理学、人工智能、宏观经济学等领域。

历史背景

蒙特卡罗方法的历史最早可以追溯到法国数学家布丰（Buffon）于1777年提出的布丰投针试验（Buffon’s Needle）^[1]。如图1所示，平面上有一组间距为$d$的平行线，试验者将一根长为$l (l<d)$的针任意投向这个平面。容易证明，在理想情况下，针与平行线相交的概率为$\frac{2l}{\pi d}$。因此，若用$p_n$表示$n$次试验中相交的频率，则可得到$\pi$的值为：
$$
\pi = \mathop {\lim }\limits_{n \to \infty } \frac{2l}{p_n d}
$$
当试验次数无限增大时该极限将收敛到$\pi$。

图1 布丰投针试验

直到上世纪40-50年代电子计算机的出现，蒙特卡罗方法才真正在科学计算领域得以应用。第二次世界大战期间，在研制原子弹的曼哈顿计划中，为了计算中子在核燃料中的随机扩散问题，以及求解临界质量，Stanislaw Ulam和John von Neumann等人提出了一种基于统计抽样的数值计算方法。作为本方法的重要推广者，Nick Metropolis用摩纳哥的著名赌城Monte Carlo为其命名，暗示了它的概率论与统计学基础。

基本原理

很多实际科学计算问题均可归结于计算如下的积分：
$$
I = \int_D {g({\bf{x}}){\rm{d}}{\bf{x}}}
$$
其中积分区域$D$通常为一高维空间$\mathbb{R}^n$的子集，$ g({\bf{x}})$为我们感兴趣的目标函数。如果我们从$D$中均匀的抽取$m$个样本$\left\{ {\bf{x}}^{(1)},{\bf{x}}^{(2)}, \cdots , {\bf{x}}^{(m)} \right\}$，设$I$的一个估计值为：
$$
{\hat I_m} = \frac{\xi}{m}\left[ g \left({\bf{x}}^{(1)}\right) + g \left({\bf{x}}^{(2)}\right) + \cdots + g \left({\bf{x}}^{(m)}\right) \right]
$$
式中：
$$
\xi = \int_D {\rm{d} \bf{x}}
$$
若$n$维随机变量$\bf{X}$服从$D$上的独立均匀分布，则$\bf{X}$的概率密度函数$f(\bf{x})$为：
$$
f(\bf{x}) = \left\{ \begin{array}{*{20}{c}}
{\xi ^{ - 1}}&{ {\bf{x}} \in D}\\
0&{ {\bf{x}} \notin D}
\end{array} \right.
$$
设随机变量$Y=g({\bf{X}})$，根据定义，$Y$的数学期望为：
$$
\begin{aligned}
E(Y) &= E\left[ {g({\bf{X}})} \right] \\
&= \int_D {f({\bf{x}})g({\bf{x}}){\rm{d}}{\bf{x}}} \\
&= \frac{1}{\xi }\int_D {g({\bf{x}}){\rm{d}}{\bf{x}}} = \frac{I}{\xi }
\end{aligned}
$$
根据大数定理^[2]，对于$\forall \varepsilon > 0$，有：
$$
\mathop {\lim }\limits_{m \to \infty } P\left( {\left| {\frac{g \left({\bf{x}}^{(1)}\right) + g \left({\bf{x}}^{(2)}\right) + \cdots + g \left({\bf{x}}^{(m)}\right)}{m} - E(Y)} \right| < \varepsilon } \right) = 1
$$
即：
$$
\mathop {\lim }\limits_{m \to \infty } P\left( \left| {\frac{ {\hat I}_m - I}{\xi}} \right| < \varepsilon \right) = 1
$$

$$
\mathop {\lim }\limits_{m \to \infty } P\left( {\left| { {\hat I}_m - I} \right| < \varepsilon ^ {\prime}} \right) = 1 , \varepsilon ^ {\prime} = \varepsilon \left| \xi \right|
$$
因此当$m$无限增大时，估计值$\hat I_m$依概率收敛于$I$。根据中心极限定理^[2]，当$m \to \infty$时：
$$
\frac{ {\hat I}_m - I}{\xi } \to N\left( {0,\frac{\sigma ^ 2}{m}} \right)
$$
式中：$\sigma = D(g({\bf{x}}))$。因此，该蒙特卡罗方法的误差为$O\left(m^{ - 1/2}\right)$。

虽然针对不同问题的蒙特卡罗方法在具体实现上有所差异，但其基本过程无外乎如下几步：

先构造随机过程，定义所有可能出现结果的范围（即样本空间）；
在已经定义的样本空间中按已知概率分布抽样；
对抽样结果进行确定性计算，以得到各种所需的统计量；
根据得出的统计量来估计概率模型的特征值，以得到问题的解。

下面将通过几个计算实例来讨论蒙特卡罗方法的具体实现过程。

计算实例

计算圆周率$\pi$

如图2所示，正方形边长为1，显然有正方形的面积$S=1$，四分之一圆的面积$S_o = \frac{\pi}{4}$。若在正方形中随机取一点$\left( x_0, y_0 \right)$，设事件$A$为点$\left( x_0, y_0 \right)$落在四分之一圆内，则由几何概型可知事件$A$发生的概率为：
$$
P(A) = \frac{S_o}{S} = \frac{\pi }{4}
$$
我们可以进行$n$次重复试验，并统计事件$A$发生的次数$m$，得到事件$A$发生的频率$p_n=\frac{m}{n}$。根据大数定理^[2]，当$n$足够大时，可以用事件$A$发生的频率来近似$P(A)$。故可得到$\pi$的估计值：
$$
{\hat \pi} = 4p_n
$$

图2 计算圆周率$\pi$

在此随机过程中样本空间为：
$$
\Omega = \left\{ (x,y)|0 \le x \le 1,0 \le y \le 1 \right\}
$$
表示事件$A$的区域：
$$
A = \left\{ (x,y)|{x^2} + {y^2} \le 1 \right\}
$$
在$\Omega$中均匀抽样，统计样本点落在$A$内的频率，进而得到$\pi$的估计值。可将此算法编写为如下C语言程序：

程序清单：mc_pi.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double rand01()
/* Generate random number between 0 - 1 */
{
    return ((double)(rand() % 1001)) / 1000.0;
}

int main()
{
    int n;
    int i;
    int m = 0;
    double x, y;
    double pi;
    
    printf("Input NR test: ");
    scanf("%d", &n);
    
    for(i = 0; i < n; i++)
    {
        x = rand01();
        y = rand01();
        if(x * x + y * y <= 1.0)
            m++;
    }
    
    pi = 4.0 * (double)m / (double)n;
    printf("pi = %f\n", pi);
    
    return 0;
}

计算结果如下：

试验次数n	$100$	$1000$	$10^4$	$10^5$	$10^6$	$10^7$
估计值${\hat \pi}$	3.120	3.100	3.108	3.135	3.137	3.139

计算定积分

投点法

若当$a \le x \le b$时，$0 \le g(x) \le h$，定积分$I = \int_a^b {g(x){\rm{d}}x}$的值可表示为图3中阴影部分曲边梯形的面积。

图3 投点法计算定积分

设随机变量$(X,Y)$服从矩形区域$\Omega = \left\{ {(x,y)|a \le x \le b,0 \le y \le h} \right\}$上的均匀分布，且$X$与$Y$相互独立，设事件$A$为：
$$
A = \left\{ {Y \le g\left( X \right)} \right\}
$$
则事件$A$发生的概率为：
$$
P(A) = \frac{I}{h(b - a)}
$$
参照实例1，可以用$n$次重复试验中事件$A$出现的频率$p_n=\frac{m}{n}$来作为$P(A)$的近似值，故可得到定积分的估计值：
$$
\hat I = \frac{mh(b - a)}{n}
$$
若给定$g(x) = \sin x$，$a = 0$，$b = \pi $，取$h = 1$，则可用如下程序求得定积分$\int_0^\pi {\sin x{\rm{d}}x} $的近似值。

程序清单：mc_int.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double rand01()
/* Generate random number between 0 - 1 */
{
    return ((double)(rand() % 1001)) / 1000.0;
}

int main()
{
    int n;
    int i;
    int m = 0;
    double x, y;
    double s;
    
    printf("Input NR test: ");
    scanf("%d", &n);
    
    for(i = 0; i < n; i++)
    {
        x = M_PI * rand01();
        y = rand01();
        if(y <= sin(x))
            m++;
    }
    
    s = M_PI * (double)m / (double)n;
    printf("I = %f\n", s);
    
    return 0;
}

计算结果如下：

试验次数n	$100$	$1000$	$10^4$	$10^5$	$10^6$	$10^7$
估计值${\hat I}$	2.136	1.945	1.962	1.988	1.997	1.998

平均值法

设随机变量$X$服从$[a,b]$上的均匀分布，则随机变量$X$的函数$Y = g(X)$的数学期望为：
$$
E(Y) = E\left[ {g(X)} \right] = \frac{1}{b - a}\int_a^b {g(x){\rm{d}}x}
$$
在$[a,b]$上均匀抽样，根据大数定理^[2]，当样本点足够多时，可以用样本平均值$\frac{1}{n}\sum\limits_{i = 1}^n {g\left( {x_i} \right)}$来近似数学期望$E(Y)$，进而得到定积分$I = \int_a^b {g(x){\rm{d}}x}$的估计值为：
$$
\hat I = \frac{b - a}{n}\sum\limits_{i = 1}^n {g\left( {x_i} \right)}
$$
注意到，如果取${x_i} = a + \frac{i}{n}(b - a)$，则上式将变成Riemann积分^[3]的形式。
而实际上，平均值法亦有着较为直观的几何意义。在图4中，阴影部分面积即为定积分$\int_a^b {g(x){\rm{d}}x}$的值。设$\bar y$为$g(x)$在闭区间$[a,b]$上的平均值，则阴影部分面积可表示为：
$$
S = \int_a^b {g(x){\rm{d}}x} = \bar y(b - a)
$$

图4 平均值法的几何意义

仍以定积分$\int_0^\pi {\sin x{\rm{d}}x}$为例，平均值法的C语言实现为：

程序清单：mc_inta.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double rand01()
/* Generate random number between 0 - 1 */
{
    return ((double)(rand() % 1001)) / 1000.0;
}

int main()
{
    int n;
    int i;
    double x, y;
    double s = 0;
    
    printf("Input NR test: ");
    scanf("%d", &n);
    
    for(i = 0; i < n; i++)
    {
        x = M_PI * rand01();
        y = sin(x);
        s += y;
    }
    
    s = M_PI * s / (double)n;
    printf("I = %f\n", s);
    
    return 0;
}

计算结果如下：

试验次数n	$100$	$1000$	$10^4$	$10^5$	$10^6$	$10^7$
估计值${\hat I}$	2.096	2.032	1.991	1.997	1.996	1.998

计算二重积分

我们知道二重积分的几何意义可表示为曲顶柱体的体积。若当$(x,y) \in D$时，$0 \le g(x,y) \le h$，则二重积分$I = \iint\limits_D {g(x,y){\rm{d}}x{\rm{d}}y}$可表示为图5中的曲顶柱体体积$V$。

图5 二重积分的几何意义

将示例2中的投点法由二维拓展到三维，若有界区域
$$
D \subset \left\{ {(x,y)|a \le x \le b,c \le y \le d} \right\}
$$
设随机变量$(X,Y,Z)$服从区域
$$
\Omega = \left\{ {(x,y,z)|a \le x \le b,c \le y \le d,0 \le z \le h} \right\}
$$
上的均匀分布，且$X$，$Y$，$Z$两两相互独立，设事件$A$为：
$$
A = \left\{ {(X,Y) \in D \wedge Z \le f(X,Y)} \right\}
$$
则事件$A$发生的概率为：
$$
P(A) = \frac{I}{h(b-a)(d-c)}
$$
采用$n$次重复试验中$A$出现的频率${p_n} = \frac{m}{n}$来近似$P(A)$，则可得二重积分的估计值为：
$$
\hat I = \frac{mh(b - a)(d - c)}{n}
$$
下面以二重积分$\iint\limits_D { {\rm{e}}^{x+y}{\rm{d}}x{\rm{d}}y}$为例，积分区域为：
$$
D = \left\{ {(x,y)|\left| x \right| + \left| y \right| \le 1} \right\}
$$
积分区域$D$如图6中阴影区域所示：

图6 积分区域$D$

该积分的解析解为：
$$
\begin{aligned}
I &= \iint\limits_D { {\rm{e}}^{x+y}{\rm{d}}x{\rm{d}}y} \\
&= \int_{ - 1}^0 { {\rm{e}}^x{\rm{d}}x} \int_{ - x - 1}^{x + 1} { {\rm{e}}^y{\rm{d}}y} + \int_0^1 { {\rm{e}}^x{\rm{d}}x\int_{x - 1}^{ - x + 1} { {\rm{e}}^y{\rm{d}}y} } \\
&= {\rm{e}} - \frac{1}{\rm{e}} = 2.3504 \cdots
\end{aligned}
$$

易知当$(x,y) \in D$时，有$x + y \le 1$，由于指数函数在定义域$\mathbb{R}$上为严格单调增函数，故有$0 < { {\rm{e}}^{x + y}} \le \rm{e}$。当采用上述蒙特卡罗方法时，可设$a = - 1$，$b = 1$，$c = - 1$，$d = 1$，$h = \rm{e}$即：
$$
D \subset \left\{ {(x,y)| - 1 \le x \le 1, - 1 \le y \le 1} \right\}
$$
样本空间为：
$$
\Omega = \left\{ {(x,y,z)| - 1 \le x \le 1, - 1 \le y \le 1,0 \le z \le {\rm{e}}} \right\}
$$
将上述算法编写为C语言程序，如下：

程序清单：mc_int2.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double rand01()
/* Generate random number between 0 - 1 */
{
    return ((double)(rand() % 1001)) / 1000.0;
}

int main()
{
    int n;
    int i;
    int m = 0;
    double x, y, z;
    double s;
    
    printf("Input NR test: ");
    scanf("%d", &n);
    
    for(i = 0; i < n; i++)
    {
        x = 2.0 * rand01() - 1.0;
        y = 2.0 * rand01() - 1.0;
        z = M_E * rand01();
        if(fabs(x) + fabs(y) <= 1 && z <= exp(x + y))
            m++;
    }
    
    s = 4.0 * M_E * (double)m / (double)n;
    printf("S = %f\n", s);
    
    return 0;
}

计算结果如下：

试验次数n	$100$	$1000$	$10^4$	$10^5$	$10^6$	$10^7$
估计值${\hat I}$	2.936	2.349	2.314	2.340	2.341	2.350

本方法亦可拓展到三维或更高维空间，用于计算三重积分或多重积分。事实上，蒙特卡罗方法是求解高维问题的一种较为有效的方法。这是因为当采用确定性数值积分算法时，随着维度的增加，计算量将呈指数增长，这将导致高维灾难（Curse of Dimensionality）^[4]的发生；且多重积分的边界条件往往非常复杂，很难转化为累次定积分。而蒙特卡罗方法可以有效的避免上述问题的产生。

总结

蒙特卡罗方法具有如下优点：

思路直观清晰，无需复杂数学推导，程序简洁明了；
不受问题多维度与复杂边界条件的影响，是解决高维问题的有效方法；
具有模拟性质，便于对真实物理过程进行仿真分析。

蒙特卡罗方法的主要缺点是收敛速度较慢且以概率收敛（即增加抽样数量，误差并非确定性减小），若单纯以增加抽样数量来减小误差，就要增加较大的计算量。但不论如何，蒙特卡罗方法仍然为解决一些复杂计算问题提供了新的思路。

参考文献

[1] Wikipedia contributors, “Buffon’s needle,” Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/w/index.php?title=Buffon%27s_needle&oldid=671661324
[2] 盛骤，谢式千，潘承毅. 概率论与数理统计. 北京：高等教育出版社，2008
[3] 伍胜健. 数学分析. 北京：北京大学出版社，2009
[4] Wikipedia contributors, “Curse of dimensionality,” Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/w/index.php?title=Curse_of_dimensionality&oldid=667550420

本文作者：赵宝华
本文链接：https://blog.zbhitech.com/introduction-to-mc-method/
版权声明：本文采用CC BY-NC-ND 4.0协议授权，转载请注明出处与原作者信息，请勿用于商业用途或修改原文内容。