蒙特卡洛方法,也就是统计模拟方法,以摩纳哥的著名赌城蒙特卡洛命名,该方法求解问题的基本步骤为:
- 构造概率模型,使得待求解问题恰好是概率模型的某个参数,比如概率模型的期望
- 依据构造的概率模型生成样本
- 由样本建立估计量,给出问题近似解
下面以圆周率和自然常数的求解为例,简单介绍模特卡洛方法
圆周率 π \pi π
-
构造概率模型
- 给定一个边长为2的正方形和期内切圆,内切圆和正常性的面积比率$ \frac{内切圆面积}{正方形面积} = \frac{\pi * 1^2}{2 * 2} = \frac{\pi}{4}$
- 在正方形内部随机生成一定数量的点,则点落入内切圆的概率等于内切圆与正方形面积之比
- X , Y ∼ U ( − 1.0 , 1.0 ) X,Y \sim U(-1.0, 1.0) X,Y∼U(−1.0,1.0),则 P ( X 2 + Y 2 ≤ 1 ) = π 4 P(X^2 + Y^2 \le 1) = \frac{\pi}{4} P(X2+Y2≤1)=4π
-
抽样
- X ∼ U ( − 1.0 , 1.0 ) , Y ∼ U ( − 1.0 , 1.0 ) X \sim U(-1.0, 1.0), Y \sim U(-1.0, 1.0) X∼U(−1.0,1.0),Y∼U(−1.0,1.0),构成样本 ( x 1 , y 1 ) , . . . , ( x n , y n ) (x_1,y_1),...,(x_n,y_n) (x1,y1),...,(xn,yn),其中 n n n为样本容量
-
建立估计量
- 样本落入圆中的概率 P ( X 2 + Y 2 < 1 ) = π 4 ≈ ∑ 1 i f x i 2 + y 2 ≤ 1 n P(X^2 + Y^2 < 1) = \frac{\pi}{4}\approx \frac{\sum 1\quad if \quad x_i^2 + y^2 \le 1 }{n} P(X2+Y2<1)=4π≈n∑1ifxi2+y2≤1
- π ≈ 4 ∗ ∑ i 1 i f x i 2 + y i 2 < = 1 n = 4 ∗ m n \pi \approx 4 * \frac{ \sum_i1\quad if\quad x_i^2 + y_i^2 <= 1 }{n} = 4 * \frac{m}{n} π≈4∗n∑i1ifxi2+yi2<=1=4∗nm
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
- 构造概率模型
- ∫ 1 2 1 x d x = l n 2 − l n 1 = l n 2 \int_1^2 \frac{1}{x}dx = ln2 - ln1 = ln2 ∫12x1dx=ln2−ln1=ln2
-
∫
1
2
1
x
d
x
\int_1^2\frac{1}{x}dx
∫12x1dx为上图中阴影部分的面积,采用蒙特卡洛方式计算其面积
- 在正方形内生产均匀分布的随机点,则点落入阴影部分的概率,等于阴影面积与矩形形面积之比 p = A 的 面 积 矩 形 面 积 = ∫ 1 2 1 x d x ( 2 − 1 ) ∗ 1 p = \frac{A的面积}{矩形面积} = \frac{\int_1^2\frac{1}{x}dx}{(2-1)* 1} p=矩形面积A的面积=(2−1)∗1∫12x1dx
- 抽样
- X ∼ U ( 1.0 , 2.0 ) , Y ∼ I ( 0.0 , 1.0 ) X \sim U(1.0, 2.0), Y \sim I(0.0, 1.0) X∼U(1.0,2.0),Y∼I(0.0,1.0),抽取样本 ( x 1 , y 1 ) , . . . , ( x n , y n ) (x_1,y_1),...,(x_n,y_n) (x1,y1),...,(xn,yn), n为样本容量
- 建立估计量
- p ^ = ∑ i 1 i f y i ≤ 1 x i n = ∫ 1 2 1 x d x = l n 2 = m n \hat{p} = \frac{\sum_i 1\quad if\quad y_i \le \frac{1}{x_i}}{n} = \int_1^2\frac{1}{x}dx = ln2 = \frac{m}{n} p^=n∑i1ifyi≤xi1=∫12x1dx=ln2=nm
- l n 2 = l o g 2 2 l o g 2 e = m n ln2 = \frac{log_22}{log_2e} = \frac{m}{n} ln2=log2elog22=nm
- e = 2 n m e=2^{\frac{n}{m}} e=2mn
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)))