蒙特卡羅方法概述

1、第一章第一章 蒙特卡羅方法概述蒙特卡羅方法概述第二章第二章 隨機(jī)數(shù)隨機(jī)數(shù)第三章第三章 由已知分布的隨機(jī)抽樣由已知分布的隨機(jī)抽樣第四章第四章 蒙特卡羅方法解粒子輸運(yùn)問(wèn)題蒙特卡羅方法解粒子輸運(yùn)問(wèn)題l蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用許淑艷 編著原子能出版社l蒙特卡羅方法蒙特卡羅方法清華大學(xué)l蒙特卡羅方法及其在粒子輸運(yùn)問(wèn)題中的應(yīng)用蒙特卡羅方法及其在粒子輸運(yùn)問(wèn)題中的應(yīng)用裴鹿成 張孝澤 編著 科學(xué)出版社l蒙特卡羅方法蒙特卡羅方法徐鐘濟(jì) 編著上?？茖W(xué)技術(shù)出版社l電話電電子郵件電子郵件1. 蒙特卡羅方法的基本思想蒙特卡羅方法的基本思想2. 蒙特卡羅方法

2、的收斂性，誤差蒙特卡羅方法的收斂性，誤差3. 蒙特卡羅方法的特點(diǎn)蒙特卡羅方法的特點(diǎn)4. 蒙特卡羅方法的主要應(yīng)用范圍蒙特卡羅方法的主要應(yīng)用范圍作作 業(yè)業(yè)Monte Carlo方法：方法：蒙特卡羅方法，又稱(chēng)隨機(jī)抽樣或統(tǒng)計(jì)試驗(yàn)方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在本蒙特卡羅方法，又稱(chēng)隨機(jī)抽樣或統(tǒng)計(jì)試驗(yàn)方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在本世紀(jì)四十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來(lái)的。世紀(jì)四十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來(lái)的。亦稱(chēng)統(tǒng)計(jì)模擬方法，亦稱(chēng)統(tǒng)計(jì)模擬方法，statistical simulation method 利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法

3、Monte Carlo名字的由來(lái)：名字的由來(lái)：是由是由Metropolis在二次世界大戰(zhàn)期間提出的：在二次世界大戰(zhàn)期間提出的：Manhattan計(jì)劃，研究與原子彈有關(guān)計(jì)劃，研究與原子彈有關(guān)的中子輸運(yùn)過(guò)程；的中子輸運(yùn)過(guò)程；Monte Carlo是摩納哥（是摩納哥（monaco)的首都，該城以賭博聞名的首都，該城以賭博聞名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte Carlo模擬的應(yīng)用：模擬的應(yīng)用：自然現(xiàn)象的模擬：自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過(guò)程；宇宙射線在地球大氣中的傳輸過(guò)程；高能物理實(shí)驗(yàn)中的核相互作用過(guò)程；高能物

4、理實(shí)驗(yàn)中的核相互作用過(guò)程；實(shí)驗(yàn)探測(cè)器的模擬實(shí)驗(yàn)探測(cè)器的模擬數(shù)值分析：數(shù)值分析：利用利用Monte Carlo方法求積分方法求積分目前國(guó)外利用蒙特卡羅方法解決屏蔽問(wèn)題很普遍。蒙特卡羅方法由于直觀、逼真，應(yīng)用范圍越來(lái)越廣泛，在粒子輸運(yùn)問(wèn)題，統(tǒng)計(jì)物理，典型數(shù)學(xué)問(wèn)題，真空技術(shù)，激光技術(shù)以及醫(yī)學(xué)，生物，探礦等方面應(yīng)用很廣。在SCI 1998年論文檢索中，關(guān)于蒙特卡羅方法的論文有2892篇，EI 論文檢索中關(guān)于蒙特卡羅方法的論文也有2212篇，可見(jiàn)該方法在國(guó)內(nèi)、國(guó)際上應(yīng)用很廣。Monte CarloMonte Carlo方法簡(jiǎn)史方法簡(jiǎn)史簡(jiǎn)單地介紹一下簡(jiǎn)單地介紹一下Monte CarloMonte Carl

5、o方法的發(fā)展歷史方法的發(fā)展歷史1 1、BuffonBuffon投針實(shí)驗(yàn)：投針實(shí)驗(yàn)：17681768年，法國(guó)數(shù)學(xué)家年，法國(guó)數(shù)學(xué)家Comte de Buffon利用投針實(shí)驗(yàn)估計(jì)利用投針實(shí)驗(yàn)估計(jì) 的值的值dLp2dL2、1930年，年，Enrico Fermi利用利用Monte Carlo方法研究中子的擴(kuò)方法研究中子的擴(kuò)散，并設(shè)計(jì)了一個(gè)散，并設(shè)計(jì)了一個(gè)Monte Carlo機(jī)械裝置，機(jī)械裝置，F(xiàn)ermiac,用于計(jì)用于計(jì)算核反應(yīng)堆的臨界狀態(tài)算核反應(yīng)堆的臨界狀態(tài)3、Von Neumann是是Monte Carlo方法的正式奠基者方法的正式奠基者,他與他與Stanislaw Ulam合作建立了概率密度函

6、數(shù)、反累積分布函數(shù)合作建立了概率密度函數(shù)、反累積分布函數(shù)的數(shù)學(xué)基礎(chǔ)，以及偽隨機(jī)數(shù)產(chǎn)生器。在這些工作中，的數(shù)學(xué)基礎(chǔ)，以及偽隨機(jī)數(shù)產(chǎn)生器。在這些工作中， Stanislaw Ulam意識(shí)到了數(shù)字計(jì)算機(jī)的重要性意識(shí)到了數(shù)字計(jì)算機(jī)的重要性合作起源于合作起源于Manhattan工程：利用工程：利用ENIAC(Electronic Numerical Integrator and Computer)計(jì)算產(chǎn)額計(jì)算產(chǎn)額Monte Carlo模擬在物理研究中的作用模擬在物理研究中的作用 蒙特卡羅方法又稱(chēng)隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法。半個(gè)多世紀(jì)以來(lái)，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明 ，這種方法作為一種獨(dú)立的方法

7、被提出來(lái)，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。蒙特卡羅方法是一種計(jì)算方法，但與一般數(shù)值計(jì)算方法有很大區(qū)別。它是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。由于蒙特卡羅方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程，解決一些數(shù)值方法難以解決的問(wèn)題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。 二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨(dú)立的方法被提出來(lái)，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。但其基本思想并非新穎，人們?cè)谏a(chǎn)實(shí)踐和科學(xué)試驗(yàn)中就已發(fā)現(xiàn)，并加以利用。兩個(gè)例子 例1. 蒲豐氏問(wèn)題 例2. 射擊問(wèn)題（打靶游戲）基本思想計(jì)算機(jī)模擬試驗(yàn)過(guò)程 為了求得圓周率值，在十九世紀(jì)后期，

8、有很多人作了這樣的試驗(yàn)：將長(zhǎng)為2l的一根針任意投到地面上，用針與一組相間距離為2a（ la）的平行線相交的頻率代替概率P，再利用準(zhǔn)確的關(guān)系式： 求出值 其中為投計(jì)次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問(wèn)題。alP2)(22nNalaPl 一些人進(jìn)行了實(shí)驗(yàn)，其結(jié)果列于下表 ：實(shí)驗(yàn)者年份投計(jì)次數(shù)的實(shí)驗(yàn)值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929 設(shè)r表示射擊運(yùn)動(dòng)員的彈著點(diǎn)到靶心的距離，(r)表示擊中r處相應(yīng)的得分?jǐn)?shù)（環(huán)

9、數(shù)），f(r)為該運(yùn)動(dòng)員的彈著點(diǎn)的分布密度函數(shù)，它反映運(yùn)動(dòng)員的射擊水平。該運(yùn)動(dòng)員的射擊成績(jī)?yōu)?用概率語(yǔ)言來(lái)說(shuō)，是隨機(jī)變量(r)的數(shù)學(xué)期望，即 )(rgEg 0)()(drrfrgg 現(xiàn)假設(shè)該運(yùn)動(dòng)員進(jìn)行了次射擊，每次射擊的彈著點(diǎn)依次為r1，r2，rN，則次得分g(r1)，g(r2)，g(rN)的算術(shù)平均值 代表了該運(yùn)動(dòng)員的成績(jī)。換言之，為積分的估計(jì)值，或近似值。 在該例中，用次試驗(yàn)所得成績(jī)的算術(shù)平均值作為數(shù)學(xué)期望的估計(jì)值（積分近似值）。 NiiNrgNg1)(1 由以上兩個(gè)例子可以看出，當(dāng)所求問(wèn)題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的數(shù)學(xué)期望，或者是與概率、數(shù)學(xué)期望有關(guān)的量時(shí)，通過(guò)某種試驗(yàn)的

10、方法，得出該事件發(fā)生的頻率，或者該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值，通過(guò)它得到問(wèn)題的解。這就是蒙特卡羅方法的基本思想。 當(dāng)隨機(jī)變量的取值僅為1或0時(shí)，它的數(shù)學(xué)期望就是某個(gè)事件的概率?；蛘哒f(shuō)，某種事件的概率也是隨機(jī)變量（僅取值為1或0）的數(shù)學(xué)期望。 因此，可以通俗地說(shuō)，蒙特卡羅方法是用隨機(jī)試驗(yàn)的方法計(jì)算積分，即將所要計(jì)算的積分看作服從某種分布密度函數(shù)f(r)的隨機(jī)變量(r)的數(shù)學(xué)期望 通過(guò)某種試驗(yàn)，得到個(gè)觀察值r1，r2，rN（用概率語(yǔ)言來(lái)說(shuō)，從分布密度函數(shù)f(r)中抽取個(gè)子樣r1，r2，rN，），將相應(yīng)的個(gè)隨機(jī)變量的值g(r1)，g(r2)，g(rN)的算術(shù)平均值 作為積分的估計(jì)值（近似值

11、）。 NiiNrgNg1)(10)()(drrfrgg 為了得到具有一定精確度的近似解，所需試驗(yàn)的次數(shù)是很多的，通過(guò)人工方法作大量的試驗(yàn)相當(dāng)困難，甚至是不可能的。因此，蒙特卡羅方法的基本思想雖然早已被人們提出，卻很少被使用。本世紀(jì)四十年代以來(lái)，由于電子計(jì)算機(jī)的出現(xiàn)，使得人們可以通過(guò)電子計(jì)算機(jī)來(lái)模擬隨機(jī)試驗(yàn)過(guò)程，把巨大數(shù)目的隨機(jī)試驗(yàn)交由計(jì)算機(jī)完成，使得蒙特卡羅方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中發(fā)揮應(yīng)有的作用。 計(jì)算機(jī)模擬試驗(yàn)過(guò)程，就是將試驗(yàn)過(guò)程（如投針，射擊）化為數(shù)學(xué)問(wèn)題，在計(jì)算機(jī)上實(shí)現(xiàn)。以上述兩個(gè)問(wèn)題為例，分別加以說(shuō)明。 例1. 蒲豐氏問(wèn)題 例2. 射擊問(wèn)題（打靶游戲） 由上面兩個(gè)例題看

12、出，蒙特卡羅方法常以一個(gè)“概率模型”為基礎(chǔ)，按照它所描述的過(guò)程，使用由已知分布抽樣的方法，得到部分試驗(yàn)結(jié)果的觀察值，求得問(wèn)題的近似解。 設(shè)針投到地面上的位置可以用一組參數(shù)（x,）來(lái)描述，x為針中心的坐標(biāo)，為針與平行線的夾角，如圖所示。 任意投針，就是意味著x與都是任意取的，但x的范圍限于0，a，夾角的范圍限于0，。在此情況下，針與平行線相交的數(shù)學(xué)條件是針在平行線間的位置 sin lx 如何產(chǎn)生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均勻分布的，其分布密度函數(shù)為： 類(lèi)似地，的分布密度函數(shù)為： 因此，產(chǎn)生任意的（x,）的過(guò)程就變成了由f1(x)抽樣x及由f2()抽樣的過(guò)程了。由此得

13、到： 其中1，2均為（0,1）上均勻分布的隨機(jī)變量。 其他, 00,/1)(1axaxf其他, 00,/1)(2f21 ax 每次投針試驗(yàn)，實(shí)際上變成在計(jì)算機(jī)上從兩個(gè)均勻分布的隨機(jī)變量中抽樣得到（x,），然后定義描述針與平行線相交狀況的隨機(jī)變量s(x,)，為 如果投針次，則 是針與平行線相交概率的估計(jì)值。事實(shí)上， 于是有 其他當(dāng), 0sin, 1),(lxxsNiiiNxsNs1),(1aladxddxdfxfxsPl2)()(),(sin0021NsalaPl22產(chǎn)生任意的 過(guò)程就變成在f1中抽樣x，在f2中抽樣 的過(guò)程), x(21ax是隨機(jī)數(shù)，（0，1）間均勻分布的隨機(jī)變量投針N次，針與

14、平行線相交的隨機(jī)變量的平均只值為：N1iiiN),x( sN1sNsal 2aPl 2開(kāi)始n=n+1N=C21,21ax21sinlaN=N+1NnsNNsal 2 設(shè)射擊運(yùn)動(dòng)員的彈著點(diǎn)分布為 用計(jì)算機(jī)作隨機(jī)試驗(yàn)（射擊）的方法為，選取一個(gè)隨機(jī)數(shù)，按右邊所列方法判斷得到成績(jī)。 這樣，就進(jìn)行了一次隨機(jī)試驗(yàn)（射擊），得到了一次成績(jī) (r)，作次試驗(yàn)后，得到該運(yùn)動(dòng)員射擊成績(jī)的近似值 環(huán)數(shù) 78910概率 0.10.10.30.5環(huán)中命環(huán)命中環(huán)命中環(huán)命中1095 . 082 . 071 . 0NiiNrgNg1)(1 蒙特卡羅方法作為一種計(jì)算方法，其收斂性與誤差是普遍關(guān)心的一個(gè)重要問(wèn)題。收斂性誤差減小方

15、差的各種技巧 效率 由前面介紹可知，蒙特卡羅方法是由隨機(jī)變量X的簡(jiǎn)單子樣X(jué)1，X2，XN的算術(shù)平均值： 作為所求解的近似值。由大數(shù)定律可知， 如X1，X2，XN獨(dú)立同分布，且具有有限期望值（E(X)），則 即隨機(jī)變量X的簡(jiǎn)單子樣的算術(shù)平均值 ，當(dāng)子樣數(shù)充分大時(shí)，以概率1收斂于它的期望值E(X)。NiiNXNX111)(limXEXPNNNX 蒙特卡羅方法的近似值與真值的誤差問(wèn)題，概率論的中心極限定理給出了答案。該定理指出，如果隨機(jī)變量序列X1，X2，XN獨(dú)立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函數(shù)。則dtexXEXNPxxtNN2/221)(limdxxfXEx)()(

16、022 當(dāng)N充分大時(shí)，有如下的近似式 其中稱(chēng)為置信度，1稱(chēng)為置信水平。 這表明，不等式 近似地以概率 1成立，且誤差收斂速度的階為 。 通常，蒙特卡羅方法的誤差定義為 上式中 與置信度是一一對(duì)應(yīng)的，根據(jù)問(wèn)題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出 。122)(02/2dteNXEXPtNNXEXN)()(2/1NON 下面給出幾個(gè)常用的與的數(shù)值： 關(guān)于蒙特卡羅方法的誤差需說(shuō)明兩點(diǎn)：第一，蒙特卡羅方法的誤差為概率誤差，這與其他數(shù)值計(jì)算方法是有區(qū)別的。第二，誤差中的均方差是未知的，必須使用其估計(jì)值 來(lái)代替，在計(jì)算所求量的同時(shí)，可計(jì)算出 。 0.50.050.003 0.67451.9

17、632112)1(1NiiNiiXNXN 顯然，當(dāng)給定置信度后，誤差由和N決定。要減小，或者是增大N，或者是減小方差2。在固定的情況下，要把精度提高一個(gè)數(shù)量級(jí)，試驗(yàn)次數(shù)N需增加兩個(gè)數(shù)量級(jí)。因此，單純?cè)龃驨不是一個(gè)有效的辦法。 另一方面，如能減小估計(jì)的均方差，比如降低一半，那誤差就減小一半，這相當(dāng)于N增大四倍的效果。因此降低方差的各種技巧，引起了人們的普遍注意。后面課程將會(huì)介紹一些降低方差的技巧。 一般來(lái)說(shuō)，降低方差的技巧，往往會(huì)使觀察一個(gè)子樣的時(shí)間增加。在固定時(shí)間內(nèi)，使觀察的樣本數(shù)減少。所以，一種方法的優(yōu)劣，需要由方差和觀察一個(gè)子樣的費(fèi)用（使用計(jì)算機(jī)的時(shí)間）兩者來(lái)衡量。這就 是蒙特卡羅方法中效

18、率的概念。它定義為 ，其中c 是觀察一個(gè)子樣的平均費(fèi)用。顯然 越小，方法越有效。 c2c2優(yōu)點(diǎn)1) 能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程。2) 受幾何條件限制小。3) 收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)。4) 具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力。5) 誤差容易確定。6) 程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)。 缺點(diǎn)1) 收斂速度慢。2) 誤差具有概率性。3) 在粒子輸運(yùn)問(wèn)題中，計(jì)算結(jié)果與系統(tǒng)大小有關(guān)。 從這個(gè)意義上講，蒙特卡羅方法可以部分代替物理實(shí)驗(yàn)，甚至可以得到物理實(shí)驗(yàn)難以得到的結(jié)果。用蒙特卡羅方法解決實(shí)際問(wèn)題，可以直接從實(shí)際問(wèn)題本身出發(fā)，而不從方程或數(shù)學(xué)表達(dá)式出發(fā)。它有直觀、形象的特點(diǎn)。

19、在計(jì)算s維空間中的任一區(qū)域Ds上的積分 時(shí)，無(wú)論區(qū)域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產(chǎn)生N個(gè)點(diǎn) ，得到積分的近似值。 其中Ds為區(qū)域Ds的體積。這是數(shù)值方法難以作到的。 另外，在具有隨機(jī)性質(zhì)的問(wèn)題中，如考慮的系統(tǒng)形狀很復(fù)雜，難以用一般數(shù)值方法求解，而使用蒙特卡羅方法，不會(huì)有原則上的困難。 ssDdxdxdxxxxggs2121),( ),()()(2)(1isiixxxNiisiisNxxxgNDg1)()(2)(1),( 由誤差定義可知，在給定置信水平情況下，蒙特卡羅方法的收斂速度為，與問(wèn)題本身的維數(shù)無(wú)關(guān)。維數(shù)的變化，只引起抽樣時(shí)間及估計(jì)量計(jì)算時(shí)間的

20、變化，不影響誤差。也就是說(shuō)，使用蒙特卡羅方法時(shí)，抽取的子樣總數(shù)N與維數(shù)s無(wú)關(guān)。維數(shù)的增加，除了增加相應(yīng)的計(jì)算量外，不影響問(wèn)題的誤差。這一特點(diǎn)，決定了蒙特卡羅方法對(duì)多維問(wèn)題的適應(yīng)性。而一般數(shù)值方法，比如計(jì)算定積分時(shí)，計(jì)算時(shí)間隨維數(shù)的冪次方而增加，而且，由于分點(diǎn)數(shù)與維數(shù)的冪次方成正比，需占用相當(dāng)數(shù)量的計(jì)算機(jī)內(nèi)存，這些都是一般數(shù)值方法計(jì)算高維積分時(shí)難以克服的問(wèn)題。)(2/1NO 對(duì)于那些需要計(jì)算多個(gè)方案的問(wèn)題，使用蒙特卡羅方法有時(shí)不需要像常規(guī)方法那樣逐個(gè)計(jì)算，而可以同時(shí)計(jì)算所有的方案，其全部計(jì)算量幾乎與計(jì)算一個(gè)方案的計(jì)算量相當(dāng)。例如，對(duì)于屏蔽層為均勻介質(zhì)的平板幾何，要計(jì)算若干種厚度的穿透概率時(shí)，只需計(jì)算最厚的一種情況，其他厚度的穿透概率在計(jì)算最厚一種情況時(shí)稍加處理便可同時(shí)得到。 另外，使用蒙特卡羅方法還可以同時(shí)得到若干個(gè)所求量。例如，在模擬粒子過(guò)程中，可以同時(shí)得到不同區(qū)域的通量、能譜、角分布等，而不像常規(guī)方法那樣，需要逐一計(jì)算所求量。 對(duì)于一般計(jì)算方法，要給出計(jì)算結(jié)果與真值的誤差并不是一件容易的事情，而蒙特卡羅方法則不然。根據(jù)蒙特卡羅方法的誤差公式，可以在計(jì)算所求量的同時(shí)計(jì)算出誤差。對(duì)干很復(fù)雜的蒙特卡羅方法計(jì)算問(wèn)題，也是容易確定的。 一般計(jì)算方法常存在著有效位數(shù)損失問(wèn)題，而要解決這一問(wèn)題有時(shí)相當(dāng)困難，蒙特卡羅方法則不存在這一問(wèn)