汪 楠

（江西醫(yī)學(xué)高等?？茖W(xué)校，江西 上饒 334000）

引 言

上世紀(jì)40年代美國在研制原子彈的“曼哈頓工程”時(shí)，由烏拉姆和馮?諾伊曼首先提出蒙特卡羅方法。數(shù)學(xué)家馮?諾伊曼用賭城Monte Carlo—來命名這種方法，為它蒙上了一層神秘色彩[1]。1777年，法國數(shù)學(xué)家布豐用投針實(shí)驗(yàn)的方法求圓周率π，這被認(rèn)為是蒙特卡羅方法的起源。蒙特?卡羅方法（Monte Carlo method），也稱統(tǒng)計(jì)模擬方法，使用隨機(jī)數(shù)（或偽隨機(jī)數(shù)）來解決計(jì)算問題的方法。

可以利用蒙特卡羅方法預(yù)測誤差統(tǒng)計(jì)量，解決擴(kuò)展卡爾曼濾波器中近似誤差協(xié)方差方程，可消除由于誤差協(xié)方差方程中過于簡化的閉合，而導(dǎo)致的在擴(kuò)展卡爾曼濾波器中發(fā)現(xiàn)的無窮大誤差增長[2]。在生物遺傳方面，在蒙特卡羅模擬中引入理想組織等效正比計(jì)數(shù)器，并結(jié)合微劑量動(dòng)力學(xué)模型精確計(jì)算重離子生物有效劑量[3]。在數(shù)學(xué)方面，根據(jù)蒙特卡羅方法可用來計(jì)算反常積分[4]。在與計(jì)算機(jī)結(jié)合方面，通過建立改進(jìn)提升決策樹的近似模型，運(yùn)用算法合成新樣本點(diǎn)信息，然后采用蒙特卡羅方法預(yù)測，進(jìn)行可靠性分析，可以得出利用蒙特卡羅方法改進(jìn)梯度提升決策樹進(jìn)行可靠性分析[5-7]；在生物學(xué)方面通過構(gòu)建薄膜的三維生長動(dòng)力學(xué)模型[8]，實(shí)現(xiàn)了薄膜生長的數(shù)值模擬,可以得到了沉積速率、基底溫度和原子覆蓋度等因素對(duì)薄膜表面粗糙度的影響。在天氣氣象方面利用蒙特卡羅方法推求降水分布[9]，采用K-S檢驗(yàn)對(duì)模型進(jìn)行顯著性檢驗(yàn)，采用基于歐氏距離的層次聚類方法進(jìn)行狀態(tài)劃分，可以確定了原州區(qū)的降水分布[10]，建立了適用于原州區(qū)的滑動(dòng)平均加權(quán)馬爾可夫預(yù)測模型。蒙特?卡羅方法在金融工程學(xué)[11-12]，宏觀經(jīng)濟(jì)學(xué)，計(jì)算物理學(xué)（如粒子輸運(yùn)計(jì)算、量子熱力學(xué)計(jì)算、空氣動(dòng)力學(xué)計(jì)算）等領(lǐng)域應(yīng)用廣泛。

1 蒙特卡羅方法方法原理

蒙特卡羅方法又稱隨機(jī)抽樣方法或者隨機(jī)數(shù)法，蒙特卡羅方法是一種非常重要的數(shù)值計(jì)算方法，它是建立在以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。蒙特卡羅方法通過計(jì)算機(jī)隨機(jī)數(shù)模擬仿真來解決問題，同時(shí)可以通過模擬來檢驗(yàn)?zāi)Ｐ偷恼_性。

1.1 基本原理

一個(gè)事件的概率可以通過大量實(shí)驗(yàn)中的事件頻率來估計(jì)，事件的概率可以通過事件的頻率來表示，這個(gè)已經(jīng)通過大數(shù)定理進(jìn)行了證明，蒙特卡羅分析就是基于該理論。對(duì)影響事件可靠性的隨機(jī)變量進(jìn)行大量隨機(jī)抽樣，然后將抽樣結(jié)果代入已知的功能函數(shù)式，分析結(jié)果是否有效，最后得到結(jié)果的有效概率。

有一組獨(dú)立的隨機(jī)變量為Qi(i= 1,2,3…)，而這組隨機(jī)變量對(duì)應(yīng)的概率密度函數(shù)分別為Yq1,Yq1,Yq1…，對(duì)應(yīng)的功能函數(shù)式為F=g(q1，q2，q3…)。根據(jù)隨機(jī)變量的分布，產(chǎn)生一組隨機(jī)數(shù)q1，q2，q3，q4…，每個(gè)隨機(jī)數(shù)對(duì)應(yīng)的功能函數(shù)的函數(shù)值為Fi=g(q1,q2，q3，q4…)(i= 1,2,3…)，并將Fi記錄下來產(chǎn)生一個(gè)新的數(shù)組，那么n趨于無窮大時(shí)時(shí)，根據(jù)伯努利大數(shù)定理及正態(tài)隨機(jī)變量的特性得出所需要的結(jié)果即結(jié)構(gòu)失效概率和可靠指標(biāo)。

1.2 蒙特卡羅方法的步驟

1) 構(gòu)造概率過程；一些事件的本身就是具有隨機(jī)性質(zhì)的問題，比如分子的運(yùn)動(dòng)、原子的碰撞和粒子的運(yùn)輸問題等，對(duì)這類微觀方面的問題，主要工作是描述和模擬分子運(yùn)輸和運(yùn)動(dòng)的概率過程。對(duì)于一些確定性非隨機(jī)性質(zhì)問題，例如數(shù)學(xué)中的定積分的計(jì)算，需要構(gòu)造一個(gè)的概率過程，它的某些參數(shù)就是所求問題的解。將無隨機(jī)性問題轉(zhuǎn)化為隨機(jī)性問題。

2) 將已知概率進(jìn)行分布抽樣；建立概率模型后開始試驗(yàn)，因?yàn)楦怕誓Ｐ褪强梢钥醋魇怯筛怕史植紭?gòu)成的，所以會(huì)產(chǎn)生的概率分布的隨機(jī)變量，這個(gè)隨機(jī)變量就是實(shí)現(xiàn)蒙特卡羅方法的過程。

3) 建立恰當(dāng)?shù)墓烙?jì)量。構(gòu)造概率模型，根據(jù)模型進(jìn)行抽樣，即進(jìn)行模擬實(shí)驗(yàn)，所求隨機(jī)變量的值是解的無偏估計(jì)。建立各種估計(jì)量，就是對(duì)模型模擬實(shí)驗(yàn)的結(jié)果和數(shù)據(jù)進(jìn)行觀察和記錄，從數(shù)據(jù)中找出所需要的結(jié)果。

2 蒙特卡羅方法在積分的應(yīng)用

蒙特卡羅解決微積分的核心問題就是求解這個(gè)隨機(jī)變量的期望值，無論這個(gè)積分函數(shù)是連續(xù)分布還是離散分布。按所需要的概率分布進(jìn)行采樣，采樣次數(shù)為N，當(dāng)N→+∞時(shí)，所求的這個(gè)期望值為E(X)，隨機(jī)變量連續(xù)性和離散型的期望值分別為

2.1 定積分的計(jì)算

假設(shè)f(x)是在(x1,x2)上連續(xù)有界可積的函數(shù)，且0≤f(x)≤p，求定積分

因?yàn)閒(x)是非初等函數(shù)，無法推導(dǎo)出其原函數(shù)，即無法用常規(guī)的方法求解。在此處，將其按麥克勞林公式展開：

可以求出：

用蒙特卡羅方法通過計(jì)算機(jī)進(jìn)行計(jì)算該積分。開始向該矩形區(qū)域隨機(jī)投點(diǎn)，投點(diǎn)的總數(shù)量為N，落在曲線下方X軸上方的點(diǎn)的數(shù)量為n。該矩形區(qū)域的面積為Q=3×1=3。很容易理解，隨機(jī)投擲的點(diǎn)落在函數(shù)下方發(fā)生的概率為定積分的值除以區(qū)域總面積，則S=nQ/N。

如圖1所示，當(dāng)投點(diǎn)數(shù)N= 100時(shí)，處于曲線下方的點(diǎn)的數(shù)量為n= 43，則，數(shù)據(jù)偏差為3.594%。

圖1 N = 100時(shí)試驗(yàn)圖

如圖2所示，當(dāng)投點(diǎn)數(shù)N= 1000時(shí)，處于曲線下方的點(diǎn)的數(shù)量為n= 438，則，，數(shù)據(jù)偏差為1.80%。

圖2 N = 1000時(shí)試驗(yàn)圖

如圖3所示，當(dāng)投點(diǎn)數(shù)N= 10 000時(shí)，處于曲線下方的點(diǎn)的數(shù)量為n= 4428，則1.3284，數(shù)據(jù)偏差為0.724%；從實(shí)驗(yàn)結(jié)果可以看出，投的點(diǎn)數(shù)N越大，得到的數(shù)據(jù)更接近真實(shí)值，那么可以得出結(jié)論：當(dāng)N→∞時(shí)，實(shí)驗(yàn)得出的結(jié)果就等于真實(shí)值，數(shù)據(jù)的偏差趨于0。

圖3 N = 10 000時(shí)試驗(yàn)圖

在（3，5）之間，取n個(gè)隨機(jī)數(shù)，計(jì)算出這n個(gè)數(shù)的平均值。平均值可以近似等于g(x)在(3，5)之間的數(shù)學(xué)期望E[g(x)]。即有

編寫程序代碼，由計(jì)算機(jī)隨機(jī)生成n個(gè)隨機(jī)數(shù)，當(dāng)n取1000時(shí)，多次運(yùn)行程序得到運(yùn)行結(jié)果近似等于40.5。

2.2 計(jì)算二重積分

在(1,4)范圍內(nèi)取n個(gè)x的隨機(jī)數(shù)，在(2,5)范圍內(nèi)取n個(gè)y的隨機(jī)數(shù)，g(x,y) 為(x,y)所對(duì)應(yīng)的函數(shù)值，求出所有g(shù)(x,y)的平均值，可以近似于g(x,y)在x∈(1,4),y∈(2,5)之間的數(shù)學(xué)期望E[g(xy)]，代入公式，可得：

采用蒙特卡羅法編程計(jì)算，當(dāng)n取1000時(shí)，該二重積分的值約為217.6。

3 結(jié) 論

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，蒙特卡羅方法被廣泛應(yīng)用于各個(gè)行業(yè)中。由于現(xiàn)代計(jì)算機(jī)的軟硬件計(jì)算成本降低，在使用蒙特卡羅方法時(shí)，可以增大實(shí)驗(yàn)?zāi)M次數(shù)，從而減少計(jì)算誤差，得到更精確的值。