蒙特卡罗方法（Monte Carlo Method）_学习笔记

前言

蒙特卡罗方法（Monte Carlo Method） 是将确定性求解问题转换为随机求解问题，得到概率意义上的近似解，这为许多复杂难以求得确定性解释式的问题提供了一种近似求解的数值计算方法，在工程上的应用非常广泛。本文将简要介绍该方法原理，并使用其实现投针法求$\pi$的仿真，重点在于理解该方法的核心思想和关键要素。

背景

蒙特卡罗方法（Monte Carlo Method） 是由20世纪40年代由冯.诺伊曼( Von Neumann)和马尔钦.乌拉姆(Marcin Ulam)在从事研发美国核武器相关项目时提出的。出于对核武器项目保密的需要，他们的另一位同事康斯坦丁.梅特罗波利斯（Constantine Metropolis）建议使用“蒙特卡罗”作为他们工作的代号。这个名字源于摩纳哥的蒙特卡罗赌场，乌拉姆的叔叔会从亲戚那里借钱去那里赌博，而该方法的灵感来自于乌拉姆在养病期间玩一种名叫坎菲尔德（Canfield）单人纸牌游戏时希望计算获胜的概率。[1]值得一提的是，蒙特卡罗方法在此之前就已经存在，1777年，法国博物学家德.布丰（De Buffon）提出的布丰投针问题来求取圆周率π，被认为该方法的真正起源。

定义

目前关于蒙特卡罗方法这一术语的定义其实尚未达成一致意见，这里给出一种被普遍接受的表述。蒙特卡罗方法又被称为统计模拟，是一种随机模拟方法。以概率和统计理论方法为基础的一种数值计算方法，将求解问题与一定的概率模型相联系，多使用随机数（常见伪随机数）来实现随机抽样模拟，通过统计样本特性实现对问题的近似求解。

原理

蒙特卡罗方法的本质上是大数定律的应用，对问题的求解转换为概率意义上的求解。对于布丰投针实验，其原理即为伯努力大数定率（Bernoulli’s Law of Large Numbers）：
设${\mu }_n$是n次独立试验中事件A出现的次数，而$p(0<p<1)$是事件A在每次试验中出现的概率，则对任意$ϵ$有：
$$\begin{array}{r}\underset{n\to \infty }{\lim }P\left[\left|\frac{{\mu }_n}{n}-p\right|<ϵ\right]=1\end{array}\text{\hspace{1em}(1)}$$
此定律表明：当试验次数$n$很大时，事件A出现的频率$\frac{{\mu }_n}{n}$呈现出稳定性，逐渐稳定于事件A的概率p。此时可以用频率$\frac{{\mu }_n}{n}$代替概率p。更准确的表述是当$n$趋于无穷时，事件A出现的频率依概率收敛于事件A的概率。
下面简单分析下布丰投针实验中蒙特卡罗法原理，关于布丰投针实验详情请看视频资料[2]。
将每次投针可以看作1次试验，而每次投针的结果看作是（0-1）分布的随机变量，而针与平行线相交看作事件即结果为1。根据上述大数定律，当投针次数$n$足够多时，针与平行线相交的频率即为所求概率。而求解$\pi$的问题与所求概率建立联系的概率模型是通过几何概型设计实现的。

核心思想

核心思想：概率模型映射，不断抽样，不断统计，不断逼近，直到找到问题概率意义上的近似解。
直观理解：暴力！！！。

关键要素

蒙特卡罗方法实际上没有具体固定的形式，需要根据实际应用问题去设计，但会有如下关键要素：

1. 将求解问题映射为概率模型；
2. 确定输入可能值域范围；
3. 按照对象概率密度进行随机抽样；
4. 根据给定函数进行输出计算；
5. 对输出结果进行统计。

仿真

这里给出对布丰投针实验求$\pi$的改进版，即比较流行的正方形与圆形求$\pi$的版本。

原理阐述：如图1所示，以原点$(0,0)$为圆心分别作以半径为$r$的圆与以原点为左下角边长$r$的正方形。将$\frac{1}{4}$的原型面积记为$S_c$，正方形面积记为$S_s$，显然存在：
$$\begin{array}{r}\frac{S_c}{S_s}=\frac{1}{4}\frac{\pi R^2}{R^2}=\frac{1}{4}\pi \end{array}\text{\hspace{1em}(2)}$$
注意到$r$的取值对求解没有关系，不妨令$r=1$。改写式(2)为：
$$\begin{array}{r}\pi =4\frac{S_c}{S_s}\end{array}\text{\hspace{1em}(3)}$$
设随机变量序列 { Z 1 , Z 2 , . . . Z n } \left \{Z_{1},Z_{2},...Z_{n} \right \} {Z1​,Z2​,...Zn​}服从$(0-1)分布$表示每次投针的结果，事件A表示每次投针落入到$\frac{1}{4}$圆形区域内并取值为1，即对于第$i$次投针，随机变量$Z_i$表示为：
$$\begin{array}{rl}{Z}_i=\{\begin{array}{cc}1& d_i⩽1\\ 0& d_i>1\end{array},& i=1,\mathrm{2...}n\end{array}\text{\hspace{1em}(4)}$$
其中，$d$表示到该次投针点到原点的欧式距离。
由几何原理可知事件A的概率$P_A$为：
$$\begin{array}{r}{P}_A=\frac{S_c}{S_s}\end{array}\text{\hspace{1em}(5)}$$
根据式(1)伯努力大数定律，当$n$趋近于无穷时，事件A出现的频率收敛于其概率$P_A$，即:
$$\begin{array}{r}\underset{n\to \infty }{\lim }P\left[\left|\frac{{\sum }_{i=1}^n{Z}_i}{n}-P_A\right|<ϵ\right]=1\end{array}\text{\hspace{1em}(6)}$$
意味当投针次数$n$足够多时，事件A出现的频率近似为其概率$P_A$，即：
$$\begin{array}{r}{P}_A\approx \frac{{\sum }_{i=1}^n{Z}_i}{n}\end{array}\text{\hspace{1em}(7)}$$
由式(3)(5)(7)，求解$\pi$的问题就映射为求解概率$P_A$：
$$\begin{array}{r}\pi =4\frac{S_c}{S_s}=4P_A\approx 4\frac{{\sum }_{i=1}^n{Z}_i}{n}\end{array}\text{\hspace{1em}(8)}$$
对投针实验进行模拟，第$i$次投针在几何上存在一个落点$p_i$，为实现每次投针落点等概率，其由两个服从[0-1]均分分布且独立的随机变量$X_i$和$Y_i$组成，与随机变量$Z_i$由下式决定：
$$\begin{array}{rl}\{\begin{array}{rlr}{p}_i& =(X_i,Y_i)& \\ d_i& =\sqrt{X_i^2+Y_i^2}& \\ Z_i& =\{\begin{array}{cc}1& d_i⩽1\\ 0& d_i>1\end{array}& \end{array}& i=1,\mathrm{2...}n\end{array}\text{\hspace{1em}(9)}$$
Matlab代码如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 脚本说明：蒙特卡罗方法求Pi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bu_x = [0,1];bu_y = [1,1];   % 正方形上边界
br_x = [1,1];br_y = [1,0];   % 正方形右边界
bd_x = [1,0];bd_y = [0,0];   % 正方形下边界
bl_x = [0,0];bl_y = [0,1];   % 正方形左边界
c_X = 0:0.01:1.0;            % 圆弧上X坐标
c_Y = sqrt(1 - c_X.^2);      % 圆弧上Y坐标
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 绘制图形
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(bu_x,bu_y,'b-','LineWidth',2);
hold on;
plot(br_x,br_y,'b-','LineWidth',2);
plot(bd_x,bd_y,'m-','LineWidth',2);
plot(bl_x,bl_y,'m-','LineWidth',2);
plot(c_X, c_Y,'m-','LineWidth',2);
axis equal;
axis([0-0.005 1+0.005 0-0.005 1+0.005])
xlabel('X');
ylabel('Y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 主方法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p_i = [0,0]';
n = 5000;
z_sum = 0;
for i=1:n
  p_i(1) = rand;
  p_i(2) = rand;
  if sqrt(p_i(1)^2 + p_i(2)^2) <= 1.0
   z_sum = z_sum + 1;
   plot(p_i(1),p_i(2),'m.');
  else
   plot(p_i(1),p_i(2),'b.');
  end
  P = z_sum / i;
  r_pi = 4 * P;
  title('Monte Carlo Method Demo',['\pi = ',num2str(r_pi,'%.6f')]);
  pause(0.001);
end
hold off;

代码运行结果：

参考文献

[1] Monte Carlo Method
https://en.wikipedia.org/wiki/Monte_Carlo_method
[2] 令人惊叹的圆周率实验——蒲丰投针
https://www.bilibili.com/video/BV1uoDZYmEoL?t=3.5