基于Monte Carlo方法計算積分的研究

汪 楠

（江西醫學高等?？茖W校，江西 上饒 334000）

引 言

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

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

1 蒙特卡羅方法方法原理

蒙特卡羅方法又稱隨機抽樣方法或者隨機數法，蒙特卡羅方法是一種非常重要的數值計算方法，它是建立在以概率統計理論為基礎的一種方法。蒙特卡羅方法通過計算機隨機數模擬仿真來解決問題，同時可以通過模擬來檢驗模型的正確性。

1.1 基本原理

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

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

1.2 蒙特卡羅方法的步驟

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

2) 將已知概率進行分布抽樣；建立概率模型后開始試驗，因為概率模型是可以看作是由概率分布構成的，所以會產生的概率分布的隨機變量，這個隨機變量就是實現蒙特卡羅方法的過程。

3) 建立恰當的估計量。構造概率模型，根據模型進行抽樣，即進行模擬實驗，所求隨機變量的值是解的無偏估計。建立各種估計量，就是對模型模擬實驗的結果和數據進行觀察和記錄，從數據中找出所需要的結果。

2 蒙特卡羅方法在積分的應用

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

2.1 定積分的計算

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

因為f(x)是非初等函數，無法推導出其原函數，即無法用常規的方法求解。在此處，將其按麥克勞林公式展開：

可以求出：

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

如圖1所示，當投點數N= 100時，處于曲線下方的點的數量為n= 43，則，數據偏差為3.594%。

圖1 N = 100時試驗圖

如圖2所示，當投點數N= 1000時，處于曲線下方的點的數量為n= 438，則，，數據偏差為1.80%。

圖2 N = 1000時試驗圖

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

圖3 N = 10 000時試驗圖

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

編寫程序代碼，由計算機隨機生成n個隨機數，當n取1000時，多次運行程序得到運行結果近似等于40.5。

2.2 計算二重積分

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

采用蒙特卡羅法編程計算，當n取1000時，該二重積分的值約為217.6。

3 結 論

隨著計算機技術的快速發展，蒙特卡羅方法被廣泛應用于各個行業中。由于現代計算機的軟硬件計算成本降低，在使用蒙特卡羅方法時，可以增大實驗模擬次數，從而減少計算誤差，得到更精確的值。