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