蒙特卡洛算法

蒙特卡洛方法，也就是统计模拟方法，以摩纳哥的著名赌城蒙特卡洛命名，该方法求解问题的基本步骤为：

• 构造概率模型，使得待求解问题恰好是概率模型的某个参数，比如概率模型的期望
• 依据构造的概率模型生成样本
• 由样本建立估计量，给出问题近似解

下面以圆周率和自然常数的求解为例，简单介绍模特卡洛方法

圆周率 $\pi$

• 构造概率模型

  • 给定一个边长为2的正方形和期内切圆，内切圆和正常性的面积比率$ \frac{内切圆面积}{正方形面积} = \frac{\pi * 1^2}{2 * 2} = \frac{\pi}{4}$
  • 在正方形内部随机生成一定数量的点，则点落入内切圆的概率等于内切圆与正方形面积之比
    • $X,Y\sim U(-1.0,1.0)$，则 $P(X^2+Y^2\le 1)=\frac{\pi }{4}$

• 抽样

  • $X\sim U(-1.0,1.0),Y\sim U(-1.0,1.0)$，构成样本$(x_1,y_1),...,(x_n,y_n)$,其中$n$为样本容量

• 建立估计量

  • 样本落入圆中的概率$P(X^2+Y^2<1)=\frac{\pi }{4}\approx \frac{\sum 1if{x}_i^2+y^2\le 1}{n}$
  • $\pi \approx 4\ast \frac{{\sum }_{i}1if{x}_i^2+y_i^2<=1}{n}=4\ast \frac{m}{n}$

import matplotlib.pyplot as plt
import numpy as np

n = 10000

x = np.random.uniform(-1.0, 1.0, n)
y = np.random.uniform(-1.0, 1.0, n)


m = 0
inner_x = []
inner_y = []
outer_x = []
outer_y = []
for i in range(n):
    if x[i]*x[i] + y[i]*y[i] <= 1:
        m += 1
        inner_x.append(x[i])
        inner_y.append(y[i])
    else:
        outer_x.append(x[i])
        outer_y.append(y[i])
        
plt.figure(figsize=(8, 8))
plt.scatter(inner_x, inner_y, color='red', s = 10)
plt.scatter(outer_x, outer_y, color = 'blue', s = 10)


pi = 4 * m / n
print('n=%d, pi = %f' % (n, pi))

自然常数 e

• 构造概率模型
  • ${\int }_1^2\frac{1}{x}dx=ln2-ln1=ln2$
  • ${\int }_1^2\frac{1}{x}dx$为上图中阴影部分的面积，采用蒙特卡洛方式计算其面积
    • 在正方形内生产均匀分布的随机点，则点落入阴影部分的概率，等于阴影面积与矩形形面积之比$p=\frac{A的面积}{矩形面积}=\frac{{\int }_1^2\frac{1}{x}dx}{(2-1)\ast 1}$
• 抽样
  • $X\sim U(1.0,2.0),Y\sim I(0.0,1.0)$，抽取样本$(x_1,y_1),...,(x_n,y_n)$, n为样本容量
• 建立估计量
  • $\hat{p}=\frac{{\sum }_{i}1if{y}_i\le \frac{1}{x_i}}{n}={\int }_1^2\frac{1}{x}dx=ln2=\frac{m}{n}$
  • $ln2=\frac{lo{g}_{2}2}{lo{g}_{2}e}=\frac{m}{n}$
  • $e=2^{\frac{n}{m}}$

import matplotlib.pyplot as plt
import numpy as np


n = 1000
x = np.random.uniform(1.0, 2.0, n)
y = np.random.uniform(0.0, 1.0, n)

m = 0
inner_x = []
inner_y = []
outer_x = []
outer_y = []
for i in range(n):
    if y[i] <= 1 / x[i]:
        m += 1
        inner_x.append(x[i])
        inner_y.append(y[i])
    else:
        outer_x.append(x[i])
        outer_y.append(y[i])
        
plt.figure(figsize=(10, 10))
plt.scatter(inner_x, inner_y, color='red', s = 10)
plt.scatter(outer_x, outer_y, color = 'blue', s = 10)

print("n=%d, e = %f" % (n, np.power(2, n / m)))