蒙特·卡罗方法(Monte Carlo method)

蒙特·卡罗方法(Monte Carlo method),也称统计模拟方法,是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明,而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。

是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。与它对应的是确定性算法。

蒙特·卡罗方法在金融工程学,宏观经济学,计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。

核心原理

蒙特卡洛方法的理论基础是大数定律。

大数定律是描述相当多次数重复试验的结果的定律,根据这个定律知道样本数量越多,其平均就越趋近于真实值。

基本思想

当所求解问题是某种随机事件出现的概率,或者是某个随机变量的期望值时,通过某种“实验”的方法,以这种事件出现的频率估计这一随机事件的概率,或者得到这个随机变量的某些数字特征,并将其作为问题的解。

工作过程

蒙特卡罗方法的解题过程可以归结为三个主要步骤:构造或描述概率过程;实现从已知概率分布抽样;建立各种估计量。

蒙特卡罗方法解题过程的三个主要步骤:

构造或描述概率过程

对于本身就具有随机性质的问题,如粒子输运问题,主要是正确描述和模拟这个概率过 程,对于本来不是随机性质的确定性问题,比如计算定积分,就必须事先构造一个人为的概率过程,它的某些参量正好是所要求问题的解。

即要将不具有随机性质的问题转化为随机性质的问题。

实现从已知概率分布抽样

构造了概率模型以后,由于各种概率模型都可以看作是由各种各样的概率分布构成的,因此产生已知概率分布的随机变量(或随机向量),就成为实现蒙特卡罗方法模拟实验的基本手段,这也是蒙特卡罗方法被称为随机抽样的原因。

最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布(或称矩形分布)。随机数就是具有这种均匀分布的随机变量。

随机数序列就是具有这种分布的总体的一个简单子样,也就是一个具有这种分布的相互独立的随机变数序列。

产生随机数的问题,就是从这个分布的抽样问题。

在计算机上,可以用物理方法产生随机数,但价格昂贵,不能重复,使用不便。

另一种方法是用数学递推公式产生。这样产生的序列,与真正的随机数序列不同,所以称为伪随机数,或伪随机数序列。

不过,经过多种统计检验表明,它与真正的随机数,或随机数序列具有相近的性质,因此可把它作为真正的随机数来使用。

由已知分布随机抽样有各种方法,与从(0,1)上均匀分布抽样不同,这些方法都是借助于随机序列来实现的,也就是说,都是以产生随机数为前提的。

由此可见,随机数是我们实现蒙特卡罗模拟的基本工具。

建立各种估计量

一般说来,构造了概率模型并能从中抽样后,即实现模拟实验后,我们就要确定一个随机变量,作为所要求的问题的解,我们称它为无偏估计。建立各种估计量,相当于对模拟实验的结果进行考察和登记,从中得到问题的解。

数学应用:

通常蒙特·卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。

对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题,蒙特·卡罗方法是一种有效的求出数值解的方法。

一般蒙特·卡罗方法在数学中最常见的应用就是蒙特·卡罗积分。

拟蒙特卡洛算法

另外,拟蒙特卡洛算法在近几年也获得迅速发展。

这种方法是用确定性的超均匀分布代替蒙特卡洛算法中的随机数序列,对于某些特定问题计算速度比普通的蒙特卡洛算法高几百倍。

求解问题的关键

由此可见用拟蒙特卡罗方法求解问题的关键是如何找到一个均匀散布的点集。

目前常用的点集有GLP点集(好格子点集,good lattice point set)、GP点集(好点集,good point set)、Halton点集及其变体、Hammersley点集等。

求 pi 的值

1777年,法国数学家布丰(Georges Louis Leclere de Buffon,1707—1788)提出用投针实验的方法求圆周率π。

这被认为是蒙特卡罗方法的起源。

基本原理概述

假设有一个圆半径为1,所以四分之一圆面积就为PI,而包括此四分之一圆的正方形面积就为1。

如果随意的在正方形中投射飞标(点)好了,则这些飞标(点)有些会落于四分之一圆内,假设所投射的飞标(点)有 n点,在圆内的飞标(点)有 c 点

PI/4:1 = c:n
PI = 4*c/n

同理也可以计算椭圆的面积等等,换言之也可以计算积分。

java 代码实现

public void calcPiTest() {
    final int N = 100000;
    int count = 0;
    for (int i = 1; i < N; i++) {
        double x = Math.random();
        double y = Math.random();
        if ((x * x + y * y) < 1) {
            count++;
        }
    }
    double pi = (double) 4 * count / N;
    System.out.println("PI = " + pi);
}

输出结果:

PI = 3.14588

ps: 我尝试了提升测试的次数,其实结果并没有变得更加精确,猜测是随机数的问题。

蒙特卡罗积分

题目引入

接下来用蒙特卡洛积分求自然常数 e。

这是2015年阿里的一道笔试题。

解决方式

首先考虑方法如下:

image

接下来分别用蒙特卡洛积分和牛顿莱布尼兹公式计算,在蒙特卡洛方法中样本很多时,它们的值应该相等。

蒙特卡洛积分

利用蒙特卡洛方法,图像大致如下

image

上述积分的目的是求阴影部分的面积,所以先在所标矩形内取 n 对随机点

(x1,y1), (x2,y2), ..., (xn,yn)

对于每一对 (xi, yi),考察是否满足如下条件

1/xi ≥ yi

假设满足上述条件的点有 m 个,而全部的点有 n 个,所以得到近似公式为

image

而依据牛顿莱布尼兹公式可以得到

image

这两种方法结果应该是相等的,即有

image

java 编码

int max = 100000;
int count = 0;
for (int i = 1; i < max; i++) {
    double x = ThreadLocalRandom.current().nextDouble(1.0, 2.0);
    double y = ThreadLocalRandom.current().nextDouble(1.0);
    // 落在阴影之中
    if ((x * y) <= 1.0) {
        count++;
    }
}

double e = Math.pow(2.0, 1.0*max/count);
System.out.println(e);

其实就是随机 x,y 的值,让其满足在 (x,y) 的区间。

然后计算比例,根据公式推导 e 的值。

  • 输出
2.7199979218512764

拓展阅读

自然对数 e

《算法之美:隐匿在数据结构背后的原理》

参考资料

百度百科-蒙特卡罗方法

各种求圆周率π的算法(蒙特卡洛法的Java实现)

蒙特卡洛算法

蒲丰投针问题的推广

蒙特卡洛方法与定积分计算