MonteCarlo-method

“当时我立刻想到了一样东西：你听说过’蒙特卡洛法’吗？哦，那是一种计算不规则图形面积的计算机程序算法，具体做法是在软件中用大量的小球随机击打那块不规则图形，被击中的地方不再重复打击，这样，达到一定的数量后，图形的所有部分就会都被击中一次，这时统计图形区域内小球的数量，就得到了图形的面积，当然，球越小结果越精确。

这种方法虽然简单，却展示了数学中的一种用随机的蛮力对抗精确逻辑的思想方法，一种用数量得到质量的计算思想。这就是我解决三体问题的策略。”

《三体》中讲到了蒙特卡洛方法，挺有意思的。

随机方法

随机方法可以分为两类：

• 蒙特卡洛算法：采样越多，越近似最优解；

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

• 拉斯维加斯算法：采样越多，越有机会找到最优解；

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

蒙特卡洛算法基础

求圆周率

最简单的例子是求圆周率：

• 一个圆半径R，它有一个外切正方形边长2R。可知：

• 圆面积 = Pi*R^2，方形面积 2R x 2R=4R^2

• 从这个正方形内随机抽取一个点，对这个点的要求是在正方形内任意一点的概率平均分布。那么这个点在圆以内的概率大概就是pi*R^2/4R^2=pi/4

• 生成若干个这样的点，利用平面上两点间距离公式计算这个点到圆心的距离来判断是否在圆内。

• 当我们使用足够多的点来进行统计时，我们得到的概率值十分接近pi/4，这样就可以得到pi值

代码实现如下：

# -*- coding: utf-8 -*-
import random
import math
def main():
    print '请输入迭代的次数：'
    n = int(raw_input())   # n是随机的次数
    total = 0     # total是所有落入圆内的随机点
    for i in xrange(n):
        x = random.random()
        y = random.random()
        if math.sqrt(x**2 + y**2) < 1.0:  #判断是否落入圆内
            total += 1
        mypi = 4.0*total/n   #得到pi值
        print '迭代次数是：', n, 'Pi的值是：', mypi
        print '数学pi是：', math.pi
        print '误差是：', abs(math.pi-mypi)/math.pi

main()

蒙特卡洛算法应用

产品厚度

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

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

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

交通堵塞

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

• 当前速度是 v 。

• 如果前面没车，它在下一秒的速度会提高到 v + 1 ，直到达到规定的最高限速。

• 如果前面有车，距离为d，且 d < v，那么它在下一秒的速度会降低到 d - 1 。

• 此外，司机还会以概率 p 随机减速， 将下一秒的速度降低到 v - 1 。

上图中，横轴代表距离（从左到右），纵轴代表时间（从上到下），因此每一行就表示下一秒的道路情况。
可以看到，该模型会随机产生交通拥堵（图形上黑色聚集的部分）。这就证明了，单车道即使没有任何原因，也会产生交通堵塞。