一、蒙特卡羅方法簡介

蒙特卡羅（Monte Carlo）方法：簡單來說，蒙特卡洛的基本原理簡單描述是先大量模擬，然后計算一個事件發(fā)生的次數(shù)，再通過這個發(fā)生次數(shù)除以總模擬次數(shù)，得到想要的結(jié)果，精髓就是：用統(tǒng)計結(jié)果去計算頻率，從而得到真實值的近似值。蒙特卡洛方法可以應(yīng)用在很多場合，但求的是近似解，在模擬樣本數(shù)越大的情況下，越接近與真實值，但樣本數(shù)增加會帶來計算量的大幅上升。

二、實例

1.求圓周率pi的近似值：
（1）正方形內(nèi)部有一個相切的圓，它們的面積之比是π/4?，F(xiàn)在，在這個正方形內(nèi)部，隨機產(chǎn)生10000個點，計算它們與中心點的距離，從而判斷是否落在圓的內(nèi)部,若這些點均勻分布，則圓周率 pi=res/npi=res/n。
其中res:表示落到圓內(nèi)投點數(shù) n:表示總的投點數(shù)
（2）代碼

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# 投點次數(shù)
n = 10000
# 圓的半徑、圓心
r = 1.0
a,b = (0.,0.)
# 正方形區(qū)域
x_min, x_max = a-r, a+r
y_min, y_max = b-r, b+r
# 在正方形區(qū)域內(nèi)隨機投點
x = np.random.uniform(x_min, x_max, n) #均勻分布
y = np.random.uniform(y_min, y_max, n)
#計算點到圓心的距離
d = np.sqrt((x-a)**2 + (y-b)**2)
#統(tǒng)計落在圓內(nèi)點的數(shù)目
res = sum(np.where(d < r, 1, 0))
#計算pi的近似值（Monte Carlo:用統(tǒng)計值去近似真實值）
pi = 4 * res / n
print('pi: ',pi)
#畫個圖
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') #防止圖像變形
circle = Circle(xy=(a,b), radius=r ,alpha=0.5)
axes.add_patch(circle)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

（3）運行結(jié)果：

pi的值是：

可以看出存在一定誤差，模擬樣本越大，誤差也會隨之減小。
2.求定積分的近似值
（1）在一個1×1的正方形區(qū)域里，使用蒙特卡洛方法，隨機在這個正方形里面產(chǎn)生大量隨機點（數(shù)量為n），計算有多少點（數(shù)量為res）落在函數(shù)下方區(qū)域內(nèi)，res/n就是所要求的積分值，也即紅色區(qū)域的面積。
（2）代碼

n = 30000
#矩形邊界區(qū)域
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
#在矩形區(qū)域內(nèi)隨機投點x
x = np.random.uniform(x_min, x_max, n)#均勻分布
y = np.random.uniform(y_min, y_max, n)
#統(tǒng)計落在函數(shù)y=x^2下方的點
res = sum(np.where(y < f(x), 1 ,0))
#計算定積分的近似值
integral = res / n
print('integeral: ', integral)
# 畫個圖
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') # 防止圖像變形
axes.plot(np.linspace(x_min, x_max, 10), f(np.linspace(x_min, x_max, 10)), 'b-') # 函數(shù)圖像
#plt.xlim(x_min, x_max)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

（3）運行結(jié)果

積分值：