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