蒙特卡羅方法求解定積分

程錦華

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

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

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

1.前言

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

2.預(yù)備知識(shí)

2.1 均勻分布的隨機(jī)變量

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

P（x）=■， ?a

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

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

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

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

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

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

2.2 強(qiáng)大數(shù)定律

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

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

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

表格1：隨機(jī)變量X1的概率分布列

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

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

表格2：拋硬幣實(shí)驗(yàn)的結(jié)果

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

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

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

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

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

3.1求解定積分

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

+∞，于是應(yīng)用強(qiá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）。其次，我們考慮這種估計(jì)式的均方誤差，

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)點(diǎn)

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

（上接第142頁(yè)）

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

3.4求定積分的變式

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

例2：假設(shè)g是定義在[a，b]上的連續(xù)函數(shù)，計(jì)算■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ù)，計(jì)算■g（x）dx。

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

4.總結(jié)

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

參考文獻(xiàn)：

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

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