法不同于前面幾章所介紹的確定性數(shù)值方法

1、第八章 Monte Carlo 法§8.1 概述Monte Carlo 法不同于前面幾章所介紹的確定性數(shù)值方法，它是用來解決數(shù)學和物理問題的非確定性的（概率統(tǒng)計的或隨機的）數(shù)值方法。Monte Carlo 方法（MCM），也稱為統(tǒng)計試驗方法，是理論物理學兩大主要學科的合并：即隨機過程的概率統(tǒng)計理論（用于處理布朗運動或隨機游動實驗）和位勢理論，主要是研究均勻介質(zhì)的穩(wěn)定狀態(tài)1.R.Hersch and R.J.Griego,“Brownian motion and potential theory,”Sci.Amer.,Mar.1969,pp.67-74.。它是用一系列隨機數(shù)來近似解決問

2、題的一種方法，是通過尋找一個概率統(tǒng)計的相似體并用實驗取樣過程來獲得該相似體的近似解的處理數(shù)學問題的一種手段。運用該近似方法所獲得的問題的解in spirit更接近于物理實驗結(jié)果，而不是經(jīng)典數(shù)值計算結(jié)果。普遍認為我們當前所應用的MC技術(shù)，其發(fā)展約可追溯至1944年，盡管在早些時候仍有許多未解決的實例。MCM的發(fā)展歸功于核武器早期工作期間Los Alamos（美國國家實驗室中子散射研究中心）的一批科學家。Los Alamos 小組的基礎工作刺激了一次巨大的學科文化的迸發(fā)，并鼓勵了MCM在各種問題中的應用?！癕onte Carlo”的名稱取自于Monaco（摩納哥）內(nèi)以賭博娛樂而聞名的一座城市。Mo

3、nte Carlo 方法的應用有兩種途徑：仿真和取樣。仿真是指提供實際隨機現(xiàn)象的數(shù)學上的模仿的方法。一個典型的例子就是對中子進入反應堆屏障的運動進行仿真，用隨機游動來模仿中子的鋸齒形路徑。取樣是指通過研究少量的隨機的子集來演繹大量元素的特性的方法。例如，在上的平均值可以通過間歇性隨機選取的有限個數(shù)的點的平均值來進行估計。這就是數(shù)值積分的Monte Carlo 方法。MCM已被成功地用于求解微分方程和積分方程，求解本征值，矩陣轉(zhuǎn)置，以及尤其用于計算多重積分。任何本質(zhì)上屬隨機組員的過程或系統(tǒng)的仿真都需要一種產(chǎn)生或獲得隨機數(shù)的方法。這種仿真的例子在中子隨機碰撞，數(shù)值統(tǒng)計，隊列模型，戰(zhàn)略游戲，以及其它

4、競賽活動中都會出現(xiàn)。Monte Carlo 計算方法需要有可得的、服從特定概率分布的、隨機選取的數(shù)值序列。§8.2 隨機數(shù)和隨機變量的產(chǎn)生510全面的論述了產(chǎn)生隨機數(shù)的各類方法。其中較為普遍應用的產(chǎn)生隨機數(shù)的方法是選取一個函數(shù)，使其將整數(shù)變換為隨機數(shù)。以某種方法選取，并按照產(chǎn)生下一個隨機數(shù)。最一般的方程具有如下形式： （8.1）其中 初始值或種子（） 乘法器（）增值（） 模數(shù)對于數(shù)位的二進制整數(shù)，其模數(shù)通常為。例如，對于31位的計算機即可取。這里和都是整數(shù)，且具有相同的取值范圍。所需的隨機數(shù)序便可由下式得 （8.2）該序列稱為線性同余序列。例如，若且，則該序列為 7，6，9，0，7，

5、6，9，0 （8.3）可以證明，同余序列總會進入一個循環(huán)套；也就是說，最終總會出現(xiàn)一個無休止重復的數(shù)字的循環(huán)。（8.3）式中序列周期長度為4。當然，一個有用的序列必是具有相對較長周期的序列。許多作者都用術(shù)語乘同余法和混合同余法分別指代和時的線性同余法。選取和的法則可參見6,10。這里我們只關(guān)心在區(qū)間內(nèi)服從均勻分布的隨機數(shù)的產(chǎn)生。用字符來表示這些數(shù)字，則由式（8.2）可得 （8.4）這樣僅在數(shù)組中取值。（對于區(qū)間（0,1）內(nèi)的隨機數(shù)，一種快速檢測其隨機性的方法是看其均值是否為0.5。其它檢測方法可參見3,6。）產(chǎn)生區(qū)間內(nèi)均勻分布的隨機數(shù)，可用下式 （8.5） 用計算機編碼產(chǎn)生的隨機數(shù)（利用式（8

6、.2）和（8.4）并不是完全隨機的；事實上，給定序列種子，序列的所有數(shù)字都是完全可預測的。一些作者為強調(diào)這一點，將這種計算機產(chǎn)生的序列稱為偽隨機數(shù)。但如果適當選取和，序列的隨機性便足以通過一系列的統(tǒng)計檢測。它們相對于真隨機數(shù)具有可快速產(chǎn)生、需要時可再生的優(yōu)點，尤其對于程序調(diào)試。Monte Carlo 程序中通常需要產(chǎn)生服從給定概率分布的隨機變量。該步可用6，1315中的幾種方法加以實現(xiàn)，其中包括直接法和舍去法。直接法（也稱反演法或變換法），需要轉(zhuǎn)換與隨機變量相關(guān)的累積概率函數(shù)（即：為的概率）。顯然表明，通過產(chǎn)生(0,1)內(nèi)均勻分布隨機數(shù)，經(jīng)轉(zhuǎn)換我們可得服從分布的隨機樣本。為了得到這樣的具有概率

7、分布的隨機數(shù)，不妨設，即可得 （8.6）其中具有分布函數(shù)。例如，若是均值為呈指數(shù)分布的隨機變量，且 （8.7）在中解出可得 （8.8）由于本身就是區(qū)間(0,1)內(nèi)的隨機數(shù)，故可簡寫為 （8.9）有時(8.6)式所需的反函數(shù)不存在或很難獲得。這種情況可用舍去法來處理。令為隨機變量的概率密度函數(shù)。令時的，且上界為（即：），如圖8.1所示。我們產(chǎn)生區(qū)間(0,1)內(nèi)的兩個隨機數(shù)，則 （8.10）分別為在(a,b)和(0,M)內(nèi)均勻分布的隨機數(shù)。若 （8.11）則為的可選值，否則被舍去，然后再試新的一組。如此運用舍去法，所有位于以上的點都被舍去，而位于上或以下的點都由來產(chǎn)生。 圖8.1 舍去法產(chǎn)生概率密

8、度函數(shù)為的隨機變量 例8.1 設計一子程序使之產(chǎn)生0,1之間呈均勻分布的隨機數(shù)。用該程序產(chǎn)生隨機變，其概率分布由下式給定 解：生成的子程序如圖8.2所示。該子程序中，且。應用種子數(shù)（如1234），主程序中每調(diào)用一次子程序，就會生成一個隨機數(shù)。種子數(shù)可取1到間的任一整數(shù)。 0001 C*0002 C PROGRAM FOR GENERATING RANDOM VARIABLES WITH0003 C A GIVEN PROBABILITY DISTRIBUTION0004 C* 0005 0006 DOUBLE PRECISION ISEED 0007 0008 ISEED=1234.D0 00

9、09 DO 10 I=1,100 0010 CALL RANSOM(ISEED,R) 0011 THETA=ACOSD(1.0-2.0*R) 0012 WRITE(6,*)I,THETA 0013 10 CONTINUE 0014 STOP 0015 END 0001 C* 0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN 0003 C THE INTERVAL (0,1) 0004 C* 0005 0006 SUBROUTINE RANDOM (ISEED,R) 0007 DOUBLE PRECISION ISDDE,DEL,A 0008

10、DATA DEL,A/2147483647.D0,16807.D0/ 0009 0010 ISDDE=DMOD(A*ISDDE,DEL) 0011 R=ISDDE/DEL 0012 RETURN 0013 END 圖8.2 例8.1的隨機數(shù)生成器 圖8.2的子程序只是為了說明本章所介紹的一些概念。大多數(shù)計算機都有生成隨機數(shù)的子程序。 為了生成隨機變量，令 則有 據(jù)此，一系列具有給定分布的隨機變量便可由圖8.2所示主程序中生成。§8.3 誤差計算Monte Carlo 程序給出的解按大量的檢測統(tǒng)計都達到了平均值。因此，該解中包含了平均值附近的浮動量，而且不可能達到100的置信度。要計算

11、Monte Carlo 算法的統(tǒng)計偏差，就必須采用與統(tǒng)計變量相關(guān)的各種統(tǒng)計方法。我們只簡要介紹期望值和方差的概念，并利用中心極限定理來獲得誤差估計13,16。設是隨機變量。則的期望值或均值定義為 （8.12）這里是的概率密度分布函數(shù)。如果從中取些獨立的隨機樣本，那么的估計值就表現(xiàn)為個樣本值的均值。 （8.13）是的真正的平均值，而只是的有著準確期望值的無偏估計。雖然的期望值等于，但。因此，我們還需要的值在附近的分布測度。為了估計以及在附近的的值的分布，我們需要引入的方差，其定義為與差的平方的期望值，即 （8.14）由，故有 （8.15）或者 （8.16）方差的平方根稱為標準差，即 （8.17）

12、標準差給出了在均值附近的分布測度，并由此給出了誤差幅度的階數(shù)。的標準差與的標準差的關(guān)系表示為 （8.18）該式表明，如果用根據(jù)(8.13)式由個值構(gòu)成的來估計，那么結(jié)果中在附近的擴散范圍便與成比例，且隨著樣本數(shù)的增加而降低。為了估計的擴散范圍，我們定義樣本方差為 （8.19）由此式還可看出的期望值等于。因此樣本方差是的無偏估計。將（8.19）式中平方項乘出來，便可得樣本標準差為 （8.20）當較大時，系數(shù)可設為1。 作為獲得中心極限定理的一種方法，概率論的一個基本解可考慮二次項函數(shù) （8.21）該式表明次獨立隨機試驗中有次成功的概率。(8.21)中，是一次試驗中成功的概率，且。當和都較大時，便

13、可用Stirling 公式 （8.22）因此(8.21)式可近似為正態(tài)分布17 (8.23)其中且。因此當時，中心極限定理表明，描述由點Monte Carlo算法獲得的的分布的概率密度函數(shù)是（8.23）式所示的正態(tài)分布函數(shù)。也就是說，大量隨機變量的集合趨于呈正態(tài)分布。將（8.18）式代入（8.23）式可得 （8.24）正態(tài)（高斯）分布在工程，物理以及統(tǒng)計學的各類問題中都非常有用。高斯模型的顯著特性源于中心極限定理。因此，高斯模型經(jīng)常用于其影響程度取決于由許多不規(guī)則的、浮動的元素構(gòu)成的集合的情況。例8.2中我們給出了根據(jù)中心極限定理產(chǎn)生高斯隨機變量的算法。由于樣本數(shù)是有限的，所以Monte Ca

14、rlo 算法不可能達到絕對的確定性。故而在附近估計出某一范圍或區(qū)間以確保我們估計的落入該范圍內(nèi)。假設要得到位于和之間的概率。由定義 （8.25）令可得 （8.26a）或 （8.26b）其中是誤差函數(shù)，且是標準正態(tài)差的上分位點。對（8.26）式可做如下解釋：如果用來呈現(xiàn)獨立隨機觀測值并構(gòu)建相關(guān)隨機區(qū)間的Monte Carlo 程序以較大的值反復運行，則這些隨機區(qū)間的分位點將近似等于。隨機區(qū)間稱為置信區(qū)間，稱為置信度。大多數(shù)的Monte Carlo 算法都使用誤差，這表明是在實際均值的標準差范圍內(nèi)的。由（8.26）式可得，樣本均值位于區(qū)間 （8.27）其中是標準差的個數(shù)。 在（8.26）和（8.2

15、7）式中均假設總體標準差為已知。由于這種情況很少出現(xiàn)，故必須用由（8.20）式算得的樣本標準差來估計的值，從而使學生氏分布取代正態(tài)分布。我們知道當很大，比如時，分布近似趨于正態(tài)分布。（8.26）式等價于 （8.28）其中為自由度是的學生氏分布的上分位點；其值在任何標準統(tǒng)計學課本中都有表列可查。這樣置信區(qū)間的上、下限便可由下式給出 上限 （8.29） 下限 （8.30） Monte Carlo 算法中關(guān)于誤差估計的進一步討論，可參見18,19。例8.2 利用中心極限定理產(chǎn)生一高斯（正態(tài)）分布的隨機變量。根據(jù)中心極限定理，大量均值附近的獨立隨機變量的總和，無論其個體變量的分布如何，總趨向于一種高斯

16、分布。也就是說，對于任意隨機數(shù)，均值為，方差為， （8.31）漸漸與合并為零均值、單位標準差的正態(tài)分布。若是均勻分布的隨機變量（即），則，故而 （8.32）且變量 （8.33）近似等于均值為、方差為的正態(tài)變量。小于3時近似為我們所熟知的鐘形高斯分布。為簡化計算，通常實際中設，因為這樣可消除（8.32）式中的平方根項。然而的這種取值截掉了邊界處的分布，且無法產(chǎn)生超過的值。對于分布曲線的兩端比較重要的仿真，就必須用其它方案來產(chǎn)生高斯分布（參見2022）。 這樣，要產(chǎn)生一個均值為、標準差為的高斯隨機變量，就要遵循以下步驟： （1） 生成12個均勻分布的隨機數(shù) （2） 求得 （3） 令§8.

17、4 數(shù)值積分 對于一維積分，現(xiàn)已有一些求積公式，如3.10節(jié)中所述。而對多重積分的公式則相對較少。Monte Carlo 技術(shù)對此類多重積分的重要性至少存在兩方面的原因。多重積分的求積公式非常復雜，而多重積分的MCM幾乎保持不變。Monte Carlo 積分的收斂性與維數(shù)無關(guān)，但對求積公式并非如此。人們已經(jīng)發(fā)現(xiàn)，積分的統(tǒng)計方法是計算天線問題尤其是那些與較大結(jié)構(gòu)相關(guān)的問題中的二維或三維積分的很有效的方法23。這里將討論兩類Monte Carlo 積分方法，即簡易MCM 和具有對立變量的MCM。其它類型，如成功失敗法和控制變量法，可參見2426。本節(jié)還將簡要介紹MCM在不規(guī)范積分中的應用。

18、7;設要計算積分 （8.34）其中是維空間。令是在內(nèi)均勻分布的隨機變量。則也是隨機變量，其均值由下式給出27,28 (8.35)方差為 （8.36）其中 （8.37）如果取的個獨立樣本，即，它們與具有相同的分布，且構(gòu)成平均項 （8.38）我們便會期望該平均項接近于的均值。這樣，由（8.35）和（8.38）式，即可得 （8.39）Monte Carlo 公式可用于有限區(qū)域上的任何積分。為了舉例說明，現(xiàn)將（8.39）式應用于一維和二維積分中。 對于一維積分，設 （8.40）由（8.39）式可得 （8.41）其中是區(qū)間內(nèi)隨機變量，即 （8.42） 對于二維積分，有 （8.43）則相應的Monte C

19、arlo 公式為 （8.44）其中 （8.45） （8.39）式中無偏估計值收斂的很慢，這是由于估計值方差的級數(shù)是。一種改良的方法，即控制變量法，可通過減小估計值方差來提高其準確性和收斂性。§ 具有對立變量的Monte Carlo 積分方法 術(shù)語對立變量29,30用來描述任一系列估計值，它們能互相抵消掉彼此的方差。方便起見，我們假設積分范圍為區(qū)間(0,1)。設要得到下面單重積分的估計值 （8.46）我們期望表達式與相比具有更小的方差。如果太小，那反過來就很可能太大。因此，我們定義估計值 （8.47）其中是0和1之間的隨機數(shù)。此估計值方差的級數(shù)是，這是在（8.39）式基礎上的一個極大的

20、進步。對于兩維積分 （8.48）則相應的估計值為 （8.49）根據(jù)相似性，可將該方法延拓至更高階的積分。對于不是(0,1)的區(qū)間，可應用諸如（8.41）式到（8.45）式的轉(zhuǎn)換。例如， （8.50）其中，且。由（8.47）和（8.49）式可得，當積分維數(shù)增加時，每一維用來在簡易Monte Carlo 方法上提高效率的對立變量的最小值也會增加。這樣使得簡易Monte Carlo 方法在多維運算中更具優(yōu)越性。§不規(guī)范積分 積分式 （8.51）可用Monte Carlo 仿真法進行估計31。對于概率密度為的隨機變量，其中在區(qū)間(0,)上的積分等于1，則有 （8.52）因此，為計算（8.51

21、）式中的值，首先要得到個服從同一在區(qū)間(0,)上的積分等于1的概率密度函數(shù)的獨立隨機變量。其樣本均值 （8.53）便是的估計值。例8.3 用Monte Carlo 方法計算積分 解：該積分式表示振幅和相位分別服從某一分布的圓形孔徑天線的輻射狀況。之所以選擇該式，是因為它至少是輻射積分的一部分。且其解的閉合形式亦為可得，由此便可得Monte Carlo 解的精確性。其閉合解為 其中是一階Bessel 函數(shù)。 圖8.3給出由（8.44）和（8.45）式計算該積分的一個簡單的程序，其中以及。該程序在Vax 11/780 中調(diào)用子程序RANDU來產(chǎn)生隨機數(shù)和。對于不同的值，簡易和對立變量Monte C

22、arlo 方法都可用于計算輻射積分。表8.1中對時的結(jié)果進行了精確的比較。在應用（8.49）式時，用到了下面的對應式： 表8.1 例8.3輻射積分的MC方法積分結(jié)果N簡易MCM對立變量MCM500-0.2892-j0.0742-0.2887-j0.05851000 -0.5737+j0.0808 -0.4982-j0.00802000 -0.4922-j0.0040 -0.4682-j0.00824000 -0.3999-j0.0345-0.4216-j0.03236000-0.3608-j0.0270 -0.3787-j0.04408000-0.4327-j0.0378 -0.4139-j0

23、.024110,000 -0.4229-j0.0237 -0.4121-j0.0240精確解：-0.4116+j0§8.5 位勢問題的解位勢理論與布朗運動（或隨機游動）的關(guān)系是由Kakutani在1944年首次提出的32。自此，所謂的概率位勢理論便在諸如熱傳導3338、靜電學3946以及電力工程等許多學科領(lǐng)域得到廣泛應用。不同方程式的概率解或MC解的一個基本概念就是隨機游動。不同類型的隨機游動應用不同的Monte Carlo方法。最常見的類型是固定隨機游動和自動定位隨機游動。其它不常見類型有遷離法、縮減邊界法、內(nèi)切形法以及表面密度法。 0001 C* 0002 C INTEGRATI

24、ON USING CRUDE MONTE CARLO 0003 C AND ANTITHETIC METHODS 0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM 0005 C FOR ANY MULTIDIMENSIONAL INTEGRATION 0006 C* 0007 0008 DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/ 0009 DATA A,B,C/0.0,1.0,0.0/ 0010 0011 C SPECIFY THE INTEGRAND 0012 F(RHO,PHI)=RH

25、O*CEXP(J*ALPHA*RHO*COS(PHI) 0013 0014 J=(0.0,1.0) 0015 ALPHA=5.0 0016 PIE=3.1415927 0017 D=2.0*PIE 0018 DO 30 NRUN=500,10000,500 0019 SUM1=(0.0,0.0) 0020 SUM2=(0.0,0.0) 0021 DO 10 I=1,NRUN 0022 CALL RANDU(IS1,IS2,U1) 0023 CALL RANDU(IS3,IS4,U2) 0024 X1=A+(B-A)*U1 0025 X2=C+(D-C)*U2 0026 X3=(B-A)*(1.

26、0-U1) 0027 X4=(D-C)*(1.0-U2) 0028 SUM1=SUM1+F(X1,X2) 0029 SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4) 0030 10 CONTINUE 0031 AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN) 0032 AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN) 0033 PRINT *,NRUN,AREA1,AREA2 0034 WRITE(6,*) NRUN,AREA1,AREA2 0035 WRITE(6,20) NRUN,AREA1,AREA2 0036 20 FORMAT(2X,NRUN=,I5,3X,AREA1=,F12.6,3X,F12.6,AREA2=, 0037 1 F12.6,3X,F12.6,/) 0038 30 CONTINUE 0039 STOP 0040 END 圖8.3 例8.3中用Monte Carlo方法計算二維積分的程序§ 為具體起見，假設用固定隨機游動的MCM解拉普拉斯方程 （在區(qū)域內(nèi)） （8.54a）滿足Dirichlet邊界條件 （在邊界上） （8.54b）首先將分成網(wǎng)狀結(jié)構(gòu)，并用其有限差當量替代。（8.54a