蒙特卡罗方法—举例说明（C++、python）

1.什么是蒙特卡洛方法(Monte Carlo method)

蒙特卡罗方法也称统计模拟方法，是1940年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。

20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法。因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法。
与它对应的是确定性算法。

其实蒙特卡罗算法并不是一种算法的名称，而是对一类随机算法的特性的概括。媒体说“蒙特卡罗算法打败武宫正树”，这个说法就好比说“我被一只脊椎动物咬了”，是比较火星的。实际上是ZEN的算法具有蒙特卡罗特性，或者说它的算法属于一种蒙特卡罗算法。

那么“蒙特卡罗”是一种什么特性呢？我们知道，既然是随机算法，在采样不全时，通常不能保证找到最优解，只能说是尽量找。那么根据怎么个“尽量”法儿，我们我们把随机算法分成两类：

蒙特卡罗算法：采样越多，越近似最优解；
拉斯维加斯算法：采样越多，越有机会找到最优解；

一些例子

1、挑苹果

假如筐里有100个苹果，让我每次闭眼拿1个，挑出最大的。于是我随机拿1个，再随机拿1个跟它比，留下大的，再随机拿1个……我每拿一次，留下的苹果都至少不比上次的小。拿的次数越多，挑出的苹果就越大，但我除非拿100次，否则无法肯定挑出了最大的。这个挑苹果的算法，就属于蒙特卡罗算法——尽量找好的，但不保证是最好的。

而拉斯维加斯算法，则是另一种情况。假如有一把锁，给我100把钥匙，只有1把是对的。于是我每次随机拿1把钥匙去试，打不开就再换1把。我试的次数越多，打开（最优解）的机会就越大，但在打开之前，那些错的钥匙都是没有用的。这个试钥匙的算法，就是拉斯维加斯的——尽量找最好的，但不保证能找到。

所以你看，这两个词并不深奥，它只是概括了随机算法的特性，算法本身可能复杂，也可能简单。这两个词本身是两座著名赌城，因为赌博中体现了许多随机算法，所以借过来命名。

2、机器下棋

对于机器围棋程序而言，因为每一步棋的运算时间、堆栈空间都是有限的，而且不要求最优解，所以ZEN涉及的随机算法，肯定是蒙特卡罗式的。

机器下棋的算法本质都是搜索树，围棋难在它的树宽可以达到好几百（国际象棋只有几十）。在有限时间内要遍历这么宽的树，就只能牺牲深度（俗称“往后看几步”），但围棋又是依赖远见的游戏，甚至不仅是看“几步”的问题。所以，要想保证搜索深度，就只能放弃遍历，改为随机采样——这就是为什么在没有MCTS（蒙特卡罗搜树）类的方法之前，机器围棋的水平几乎是笑话。而采用了MCTS方法后，搜索深度就大大增加了。比如，在ZEN与武宫正树九段的对局中，我们可以看这一步棋：
武宫正树九段（执白）第53步大飞，明显企图攻角，而ZEN（执黑）却直接不理，放弃整个右下角，转而把中腹走厚。这个交换究竟是否划算，就不在这里讨论了，但我们至少可以看出，ZEN敢于在此脱先，舍弃这么大的眼前利益，其搜索深度确实达到了人类专业棋手的水平。

这两类随机算法之间的选择，往往受到问题的局限。如果问题要求在有限采样内，必须给出一个解，但不要求是最优解，那就要用蒙特卡罗算法。反之，如果问题要求必须给出最优解，但对采样没有限制，那就要用拉斯维加斯算法。

3、π的计算

这个例子是，如何用蒙特卡罗方法计算圆周率π。

正方形内部有一个相切的圆，它们的面积之比是π/4。

现在，在这个正方形内部，随机产生10000个点（即10000个坐标对 (x, y)），计算它们与中心点的距离，从而判断是否落在圆的内部。

如果这些点均匀分布，那么圆内的点应该占到所有点的 π/4，因此将这个比值乘以4，就是π的值。

C++

#include <bits/stdc++.h>2 3 #define MAX_ITERS 10000004 5 using namespace std;6 7 double Rand(double L, double R)8 {9     return L + (R - L) * rand() * 1.0 / RAND_MAX;
10 }
11
12 double GetPi()
13 {14     srand(time(NULL));
15     int cnt = 0;
16     for(int i = 0; i < MAX_ITERS; i++)
17     {18         double x = Rand(-1, 1);
19         double y = Rand(-1, 1);
20         if(x * x + y * y <= 1)
21             cnt++;
22     }
23     return cnt * 4.0 / MAX_ITERS;
24 }
25
26 int main()
27 {28     for(int i = 0; i < 10; i++)
29         cout << GetPi() << endl;
30     return 0;
31 }

python

import randomdef calpai():n = 1000000r = 1.0a, b = (0.0, 0.0)x_neg, x_pos = a - r, a + ry_neg, y_pos = b - r, b + rcount = 0for i in range(0, n):x = random.uniform(x_neg, x_pos)y = random.uniform(y_neg, y_pos)if x*x + y*y <= 1.0:count += 1print (count / float(n)) * 4

4、积分的计算
上面的方法加以推广，就可以计算任意一个积分的值。

比如，计算函数 y = x2 在 [0, 1] 区间的积分，就是求出下图红色部分的面积。

这个函数在 (1,1) 点的取值为1，所以整个红色区域在一个面积为1的正方形里面。在该正方形内部，产生大量随机点，可以计算出有多少点落在红色区域（判断条件 y < x2）。这个比重就是所要求的积分值。

C++
求自然常数e，用S=∫121xdxS=\int_{1}^{2} \frac{1}{x} d xS=∫12x1dx积分

#include <bits/stdc++.h>2 3 #define MAX_ITERS 100000004 5 using namespace std;6 7 struct Point8 {9     double x, y;
10 };
11
12 double Rand(double L, double R)
13 {
14     return L + (R - L) * rand() * 1.0 / RAND_MAX;
15 }
16
17 Point getPoint()
18 {19     Point t;
20     t.x = Rand(1.0, 2.0);
21     t.y = Rand(0.0, 1.0);
22     return t;
23 }
24
25 double getResult()
26 {27     int m = 0;
28     int n = MAX_ITERS;
29     srand(time(NULL));
30     for(int i = 0; i < n; i++)
31     {32         Point t = getPoint();
33         double res = t.x * t.y;
34         if(res <= 1.0)
35             m++;
36     }
37     return pow(2.0, 1.0 * n / m);
38 }
39
40 int main()
41 {42     for(int i = 0; i < 20; i++)
43         cout << fixed << setprecision(6) << getResult() << endl;
44     return 0;
45 }

求∫01x2dx\int_{0}^{1} x^{2} d x∫01x2dx的值

def integral():n = 1000000x_min, x_max = 0.0, 1.0y_min, y_max = 0.0, 1.0count = 0for i in range(0, n):x = random.uniform(x_min, x_max)y = random.uniform(y_min, y_max)# x*x > y，表示该点位于曲线的下面。所求的积分值即为曲线下方的面积与正方形面积的比。if x*x > y:count += 1integral_value = count / float(n)print integral_value

5、交通堵塞
蒙特卡罗方法不仅可以用于计算，还可以用于模拟系统内部的随机运动。下面的例子模拟单车道的交通堵塞。

根据 Nagel-Schreckenberg 模型，车辆的运动满足以下规则。

当前速度是 v 。
如果前面没车，它在下一秒的速度会提高到 v + 1 ，直到达到规定的最高限速。
如果前面有车，距离为d，且 d < v，那么它在下一秒的速度会降低到 d - 1 。
此外，司机还会以概率 p 随机减速， 将下一秒的速度降低到 v - 1 。

在一条直线上，随机产生100个点，代表道路上的100辆车，另取概率 p 为 0.3 。

上图中，横轴代表距离（从左到右），纵轴代表时间（从上到下），因此每一行就表示下一秒的道路情况。

可以看到，该模型会随机产生交通拥堵（图形上黑色聚集的部分）。这就证明了，单车道即使没有任何原因，也会产生交通堵塞。

6、产品厚度
某产品由八个零件堆叠组成。也就是说，这八个零件的厚度总和，等于该产品的厚度。

已知该产品的厚度，必须控制在27mm以内，但是每个零件有一定的概率，厚度会超出误差。请问有多大的概率，产品的厚度会超出27mm？

取100000个随机样本，每个样本有8个值，对应8个零件各自的厚度。计算发现，产品的合格率为99.9979%，即百万分之21的概率，厚度会超出27mm。

7、证券市场
证券市场有时交易活跃，有时交易冷清。下面是你对市场的预测。

如果交易冷清，你会以平均价11元，卖出5万股。
如果交易活跃，你会以平均价8元，卖出10万股。
如果交易温和，你会以平均价10元，卖出7.5万股。

已知你的成本在每股5.5元到7.5元之间，平均是6.5元。请问接下来的交易，你的净利润会是多少？

取1000个随机样本，每个样本有两个数值：一个是证券的成本（5.5元到7.5元之间的均匀分布），另一个是当前市场状态（冷清、活跃、温和，各有三分之一可能）。

模拟计算得到，平均净利润为92, 427美元。

参考自：
阮一峰的网络日志