蒙特卡羅方法求解定積分

程錦華

【摘要】蒙特卡羅方法，又稱隨機模擬，是利用隨機數(shù)進(jìn)行統(tǒng)計試驗以確定隨機事件相應(yīng)的概率與數(shù)學(xué)期望的方法。本文首先介紹均勻分布和強大數(shù)定律的有關(guān)內(nèi)容，然后討論用蒙特卡羅方法求定積分的理論基礎(chǔ)，收斂速度以及在高維空間中的適用性。

【關(guān)鍵詞】均勻分布 ?強大數(shù)定律 ?蒙特卡羅方法

【中圖分類號】G633.6 【文獻(xiàn)標(biāo)識碼】A 【文章編號】2095-3089（2018）42-0142-02

1.前言

蒙特卡羅（Monte Carlo）是摩納哥國的世界著名賭城。第二次世界大戰(zhàn)期間，美國原子彈“曼哈頓”計劃的成員馮·諾依曼和烏拉姆對裂變物質(zhì)中子的隨機擴散進(jìn)行模擬，并以蒙特卡羅來命名這種方法。傳統(tǒng)的經(jīng)驗方法由于不能逼近真實的物理過程，所以很難得到令人滿意的結(jié)果，而蒙特卡羅方法能夠真實地模擬實際物理過程，故在解決實際問題時往往可以得到很圓滿的結(jié)果。比如18世紀(jì)法國科學(xué)家蒲豐就提出了著名的用投針試驗來估計圓周率的方法。本文也介紹了利用蒙特卡羅方法求解定積分的過程。

2.預(yù)備知識

2.1 均勻分布的隨機變量

若隨機變量的X的概率密度函數(shù)為：

P（x）=■， ?a

則稱X服從區(qū)間（a，b）上的均勻分布，記作X～U（a，b）。其分布函數(shù)為

F（x）=0， ? ? ? ? x

例1：假設(shè)我們從（0，1）區(qū)間上取點，取得的點在數(shù)軸上的坐標(biāo)記為隨機變量X，則X服從均勻分布。那么X落在[■，■]內(nèi)的概率是多少？

解：記隨機變量X的概率密度函數(shù)為p（x），則

P（■≤x≤■）=■p（x）dx=■1dx=■

實際上（0，1）上的均勻分布和高中所學(xué)的幾何概型息息相關(guān)，因為點落在[■，■]內(nèi)的概率就是此區(qū)間的長度除以（0，1）區(qū)間的長度。

2.2 強大數(shù)定律

定義2.1 假設(shè)我們有概率空間（Ω，F(xiàn)，P），定義在其上的隨機變量序列X1，X2，…，Xn…及X，滿足P■Xn=X=1，則稱序列{Xn}幾乎處處收斂到X，記為Xn→Xa.s.（n→+∞）。

定理2.1（強大數(shù)定律）令X1，X2，…是兩兩獨立同分布的隨機變量，且期望存在，即EX1<∞。令EX1=μ，Sn=X1+X2+…+Xn。則當(dāng)n→+∞時，■幾乎處處收斂到μ。

下面，我們通過一個實驗來解釋強大數(shù)定律的直觀含義?，F(xiàn)在我們要拋一枚質(zhì)地均勻的硬幣n次，若出現(xiàn)正面向上則記為1，否則記為0。令Xi（i=1，2，…，n）表示第i次實驗的結(jié)果。假設(shè)我們進(jìn)行n次重復(fù)的伯努利試驗，即可獲得一組樣本X1，X2，…，Xn。但在抽取樣本前無法預(yù)知它們的數(shù)值，因此樣本依然是隨機的，具有如下的概率分布列：

表格1：隨機變量X1的概率分布列

于是，我們可以證明X1的數(shù)學(xué)期望存在，且EX1=1/2，于是由強大數(shù)定律的結(jié)果可知■，即正面向上的次數(shù)所占的比例在n→+∞時的極限為1/2。

現(xiàn)在用C++語言進(jìn)行模擬試驗，得到拋硬幣的結(jié)果如下，

表格2：拋硬幣實驗的結(jié)果

結(jié)果表明我們的拋硬幣次數(shù)越多，正面向上的頻率就會越接近1/2，和強大數(shù)定律的理論結(jié)果相符。

3.蒙特卡羅方法在定積分中的應(yīng)用

有很多實際問題都需要計算定積分，比如在物理中討論一些規(guī)則物體的質(zhì)心、轉(zhuǎn)動慣量的計算問題等。根據(jù)微積分基本定理，對于任意一個在區(qū)間[a，b]上可積的函數(shù)f（x），如果能夠找到它的原函數(shù)F（x），則我們可以通過如下的牛頓-萊布尼茨公式求解定積分：

■f（x）dx=F（b）-F（a）

然而在實際使用這種求定積分方法的時候，往往會遇到很多困難，因為大量的函數(shù)，例如■，sinx2等找不到可以用初等函數(shù)表示的原函數(shù);因此我們有必要研究求解定積分的其他方法。

3.1求解定積分

令f是[0，1]上的一個連續(xù)函數(shù)，我們想要計算■f（x）dx，傳統(tǒng)的方法是將[0，1]區(qū)間劃分為一些小區(qū)間，然后用梯形公式或者矩形公式進(jìn)行近似求解。下面我們就介紹一種基于強大數(shù)定律的概率途徑。令{X1，X2，…，Xn…}表示一列獨立同分布服從[0，1]上均勻分布的隨機變量，則由f的有界性可知EF（X1）<

+∞，于是應(yīng)用強大數(shù)定律，可知■■■f（Xi）=Ef（X1）a.s.

而由數(shù)學(xué)期望的定義可知■f（x）dx=Ef（X1）

于是我們可以用來In（f）：=■■■■f（Xi）近似積分I（f）：=■f（x）dx

3.2收斂速度

首先，我們由期望的性質(zhì)可得EIn（f）=E■■■■f（Xi）=■■■■Ef（Xi）=I（f）。其次，我們考慮這種估計式的均方誤差，

E（In（f）-I（f））■=E■■f（Xi）-I（f）■

=■■E（f（Xi）-I（f））（f（Xj）-I（f））

=■E（f（X1）-I（f））■=■Var（f（X1））

于是我們可以得到In（f）-I（f）～■■

3.3蒙特卡羅方法的優(yōu)點

（下轉(zhuǎn)第145頁）

（上接第142頁）

當(dāng)我們處理高維空間情形的時候，求解定積分的復(fù)雜度上升。如果我們采用劃分小區(qū)間的方式，比如每個維度都分成n塊，然后分塊運用梯形公式或者矩形公式，那么在d維空間中就要有n■個小塊，運算復(fù)雜度呈指數(shù)增長。而當(dāng)我們采用蒙特卡羅方法求解定積分的時候，根據(jù)我們3.2節(jié)中的推導(dǎo)，我們發(fā)現(xiàn)收斂速度僅僅與模擬隨機點的個數(shù)有關(guān)系，而與空間維度無關(guān)，于是在高維的情形下，蒙特卡羅方法幾乎是唯一有用的途徑。

3.4求定積分的變式

根據(jù)3.1中的結(jié)果，我們已經(jīng)可以用蒙特卡羅方法求解在區(qū)間[0，1]上連續(xù)函數(shù)的積分，但是注意函數(shù)的連續(xù)性并不是必要的，我們只要求其可積性以及隨機變量f（X1）的期望存在即可，其中X1服從[0，1]上的均勻分布?，F(xiàn)在我們運用變量替換的方法推導(dǎo)其他形式積分的蒙特卡羅方法。

例2：假設(shè)g是定義在[a，b]上的連續(xù)函數(shù)，計算■g（x）dx。

解：作變量替換，令y=■，則dy=■■。我們有■g（x）dx=■g（a+（b-a）y）（b-a）dy=■h（y）dx，其中h（y）=g（a+（b-a）y）（b-a），那么新函數(shù)h是[0，1]上的連續(xù)函數(shù)，于是可以利用3.1中的結(jié)果進(jìn)行近似。

例3：假設(shè)g是定義在[0，+∞）上的連續(xù)函數(shù)，計算■g（x）dx。

解：作變量替換y=■，則dy=-■=-y■dx■， 我們有■g（x）dx=■h（y）dx， 其中h（y）=■。

4.總結(jié)

本文首先介紹了均勻分布和強大數(shù)定律，并給出了相應(yīng)的直觀解釋。然后在強大數(shù)定律的基礎(chǔ)上推導(dǎo)出了蒙特卡羅方法，并研究了它的收斂速度和在高維空間中的適用性。

參考文獻(xiàn)：

[1]茆詩松. 概率論與數(shù)理統(tǒng)計簡明教程[M].高等教育出版社，2012.

[2]Brian W K， Dennis M R. C程序設(shè)計語言[M].機械工業(yè)出版社，2004.