蒙特卡洛方法以及python實現
- 1.什么是蒙特卡洛方法(Monte Carlo method)
- 2.蒙特卡洛方法的基本思想
- 3.應用： 蒙特卡洛求定積分常見方法
- 3.1投點法：
- 3.2期望法：
- 3.3蒙特卡洛求定積分
- 4.蒙特卡洛方法python實例
1.什么是蒙特卡洛方法(Monte Carlo method)
蒙特卡羅方法也稱統計模擬方法，是1940年代中期由于科學技術的發展和電子計算機的發明，而提出的一種以概率統計理論為指導的數值計算方法。是指使用隨機數（或更常見的偽隨機數）來解決很多計算問題的方法。
20世紀40年代，在馮·諾伊曼，斯塔尼斯拉夫·烏拉姆和尼古拉斯·梅特羅波利斯在洛斯阿拉莫斯國家實驗室為核武器計劃工作時，發明了蒙特卡羅方法。因為烏拉姆的叔叔經常在摩納哥的蒙特卡洛賭場輸錢得名，而蒙特卡羅方法正是以概率為基礎的方法。
2.蒙特卡洛方法的基本思想
百度百科:
蒙特卡羅法也稱統計模擬法、統計試驗法。是把概率現象作為研究對象的數值模擬方法。是按抽樣調查法求取統計值來推定未知特性量的計算方法。蒙特卡羅是摩納哥的著名賭城，該法為表明其隨機抽樣的本質而命名。故適用于對離散系統進行計算仿真試驗。在計算仿真中，通過構造一個和系統性能相近似的概率模型，并在數字計算機上進行隨機試驗，可以模擬系統的隨機特性。
通常蒙特卡羅方法可以粗略地分成兩類：一類是所求解的問題本身具有內在的隨機性，借助計算機的運算能力可以直接模擬這種隨機的過程。例如在核物理研究中，分析中子在反應堆中的傳輸過程。中子與原子核作用受到量子力學規律的制約，人們只能知道它們相互作用發生的概率，卻無法準確獲得中子與原子核作用時的位置以及裂變產生的新中子的行進速率和方向。科學家依據其概率進行隨機抽樣得到裂變位置、速度和方向，這樣模擬大量中子的行為后，經過統計就能獲得中子傳輸的范圍，作為反應堆設計的依據。
另一種類型是所求解問題可以轉化為某種隨機分布的特征數，比如隨機事件出現的概率，或者隨機變量的期望值。通過隨機抽樣的方法，以隨機事件出現的頻率估計其概率，或者以抽樣的數字特征估算隨機變量的數字特征，并將其作為問題的解。這種方法多用于求解復雜的多維積分問題
3.應用： 蒙特卡洛求定積分常見方法
3.1投點法：
這個方法也常常被用來求
π π
π
值。現在我們用它來求函數的定積分。如下圖所示，有一個函數
f ( x ) {f(x)}
f
(
x
)
，若要求它從a到b的定積分，其實就是求曲線下方的面積。
這時我們可以用一個比較容易算得面積的矩型罩在函數的積分區間上（假設其面積為
A r e a Area
A
r
e
a
）。然后隨機地向這個矩形框里面投點，其中落在函數
f ( x ) {f(x)}
f
(
x
)
下方的點為綠色，其它點為紅色。然后統計綠色點的數量占所有點（紅色+綠色）數量的比例為
r r
r
，那么就可以據此估算出函數
f ( x ) {f(x)}
f
(
x
)
從aa到bb的定積分為
A r e a × r Area×r
A
r
e
a
×
r
。
注意由蒙特卡洛法得出的值并不是一個精確值，而是一個近似值。而且當投點的數量越來越大時，這個近似值也越接近真實值。
3.2期望法：
下面我們來重點介紹一下利用蒙特卡洛法求定積分的第二種方法——期望法，有時也稱為平均值法。詳情移步該博主文章：https://blog.csdn.net/baimafujinji/article/details/53869358
3.3蒙特卡洛求定積分
蒙特卡洛方法的一個重要應用就是求定積分。來看下面的一個例子(參考文獻1)
當我們在[a,b]之間隨機取一點x時，它對應的函數值就是f(x)。接下來我們就可以用f(x) * (b - a)來粗略估計曲線下方的面積，也就是我們需要求的積分值，當然這種估計（或近似）是非常粗略的,如下圖所示：
進行四次隨機采樣如下：
進行四次進行四次隨機采樣如下：
在此圖中，做了四次隨機采樣，得到了四個隨機樣本
x 1 , x 2 , x 3 , x 4 {x_1},{x_2},{x_3},{x_4}
x
1
?
,
x
2
?
,
x
3
?
,
x
4
?
，并且得到了這四個樣本的f(xi)f(xi)的值分別為
f ( x 1 ) , f ( x 2 ) , f ( x 3 ) , f ( x 4 ) f(x_1), f(x_2), f(x_3), f(x_4)
f
(
x
1
?
)
,
f
(
x
2
?
)
,
f
(
x
3
?
)
,
f
(
x
4
?
)
對于這四個樣本，每個樣本能求一個近似的面積值，大小為
f ( x i ) ? ( b ? a ) f({x_i})*(b - a)
f
(
x
i
?
)
?
(
b
?
a
)
。
為什么能這么干么？對照圖下面那部分很容易理解，每個樣本都是對原函數
f ( x ) f({x})
f
(
x
)
的近似，所以我們認為
f ( x ) f({x})
f
(
x
)
的值一直都等于
f ( x i ) f({x_i})
f
(
x
i
?
)
。
按照圖中的提示，求出上述面積的數學期望，就完成了蒙特卡洛積分。
如果用數學公式表達上述過程：
S = 1 4 ( f ( x 1 ) ( b ? a ) + f ( x 2 ) ( b ? a ) + f ( x 3 ) ( b ? a ) + f ( x 4 ) ( b ? a ) ) = 1 4 ( b ? a ) ( f ( x 1 ) + f ( x 2 ) + f ( x 3 ) + f ( x 4 ) ) = 1 4 ( b ? a ) ∑ i = 1 4 f ( x i ) \begin{aligned} S & = \frac{1}{4}(f(x_1)(b-a) + f(x_2)(b-a) + f(x_3)(b-a) + f(x_4)(b-a)) \\ & = \frac{1}{4}(b-a)(f(x_1) + f(x_2) + f(x_3) + f(x_4)) \\ & = \frac{1}{4}(b-a) \sum_{i=1}^4 f(x_i) \end{aligned}
S
?
=
4
1
?
(
f
(
x
1
?
)
(
b
?
a
)
+
f
(
x
2
?
)
(
b
?
a
)
+
f
(
x
3
?
)
(
b
?
a
)
+
f
(
x
4
?
)
(
b
?
a
)
)
=
4
1
?
(
b
?
a
)
(
f
(
x
1
?
)
+
f
(
x
2
?
)
+
f
(
x
3
?
)
+
f
(
x
4
?
)
)
=
4
1
?
(
b
?
a
)
i
=
1
∑
4
?
f
(
x
i
?
)
?
對于更一般的情況，假設要計算的積分如下：
I = ∫ a b f ( x ) d x I = \int _a ^ b f(x) dx
I
=
∫
a
b
?
f
(
x
)
d
x
取N趨于無窮大，則有：
I ￣ = b ? a N ∑ i = 1 N f ( X i ) \overline I = \frac{b-a}{N}\sum_{i=1}^Nf(X_i) I = N b ? a ? i = 1 ∑ N ? f ( X i ? )
此時可以認為： I ￣ ≈ I \overline I \approx I I ≈ I
如果從直觀上理解這個式子也非常簡潔明了：
在[a,b]區間上按均勻分布取N個隨機樣本
x i {x_i}
x
i
?
,計算
f ( x 1 ) f(x_1)
f
(
x
1
?
)
并取均值，得到的相當于Y坐標值，然后乘以
( b ? a ) (b-a)
(
b
?
a
)
為X坐標長度，得到的即為對應矩形的面積，即積分值。
4.蒙特卡洛方法python實例
首先看一個經典的用蒙特卡洛方法求π值。
簡單介紹下思路：
正方形內部有一個相切的圓，它們的面積之比是π/4。現在，在這個正方形內部，隨機產生n個點，計算它們與中心點的距離，并且判斷是否落在圓的內部。若這些點均勻分布，則圓周率 pi=4 * count/n, 其中count表示落到圓內投點數 n:表示總的投點數。
//求π值
import
random
n
=
1000000
r
=
1.0
a
,
b
=
(
0.0
,
0.0
)
x_neg
,
x_pos
=
a
-
r
,
a
+
r
y_neg
,
y_pos
=
b
-
r
,
b
+
r
count
=
0.0
for
i
in
range
(
0
,
n
)
:
x
=
random
.
uniform
(
x_neg
,
x_pos
)
y
=
random
.
uniform
(
y_neg
,
y_pos
)
if
x
*
x
+
y
*
y
<=
1.0
:
count
+=
1
result
=
(
count
/
float
(
n
)
)
*
4
print
(
"result is "
,
result
)
運行結果：
result is 3.14036
再舉個栗子
看一個求定積分的例子。
假設我們想求
∫ 0 1 x 2 d x \int_0^1 x^2 dx
∫
0
1
?
x
2
d
x
的值。具體代碼如下:
（x*x > y，表示該點位于曲線的下面。所求的積分值即為曲線下方的面積與正方形面積的比。）
import
random
def
zero2onetox2
(
)
:
n
=
1000000
x_min
,
x_max
=
0.0
,
1.0
y_min
,
y_max
=
0.0
,
1.0
count
=
0
for
i
in
range
(
0
,
n
)
:
x
=
random
.
uniform
(
x_min
,
x_max
)
y
=
random
.
uniform
(
y_min
,
y_max
)
if
x
*
x
>=
y
:
count
+=
1
result
=
count
/
float
(
n
)
print
(
"result is "
,
result
)
然后運行：
zero2onetox2
(
)
運行結果：
result is 0.333084
由此可見，基于蒙特卡洛原理的求解值與解析值的結果誤差還是比較小的。
參考文獻：
此部分代碼已上傳github：https://github.com/ShuaiWang-Code/Monte_Carlo
1.https://blog.csdn.net/baimafujinji/article/details/53869358
2.http://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration
3.https://blog.csdn.net/bitcarmanlee/article/details/82716641 參考1中的期望法
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
