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 1000000
4
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 random
def calpai():
n = 1000000
r = 1.0
a, b = (0.0, 0.0)
x_neg, x_pos = a - r, a + r
y_neg, y_pos = b - r, b + r
count = 0
for 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 += 1
print (count / float(n)) * 4
4、积分的计算
上面的方法加以推广,就可以计算任意一个积分的值。
比如,计算函数 y = x2 在 [0, 1] 区间的积分,就是求出下图红色部分的面积。
这个函数在 (1,1) 点的取值为1,所以整个红色区域在一个面积为1的正方形里面。在该正方形内部,产生大量随机点,可以计算出有多少点落在红色区域(判断条件 y < x2)。这个比重就是所要求的积分值。
C++
求自然常数e,用
积分
#include <bits/stdc++.h>
2
3 #define MAX_ITERS 10000000
4
5 using namespace std;
6
7 struct Point
8 {
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 }
求 的值
def integral():
n = 1000000
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
count = 0
for 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 += 1
integral_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美元。
参考自:
阮一峰的网络日志