蒙特卡羅(Monte Carlo)模擬

1、17. 蒙特卡羅(Monte Carlo)模擬蒙特卡羅方法，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。該方法源于美國在WWI后研制原子彈的“曼哈頓計(jì)劃”。 S. M. Ulam (1908-1984)和該計(jì)劃的主持人之一、數(shù)學(xué)家馮諾伊曼（J.von Neumann ）用馳名世界的賭城摩納哥的Monte Carlo來命名這種方法。其基本思想很早以前就被人們所發(fā)現(xiàn)和利用。17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來決定事件的“概率”。19世紀(jì)人們用投針試驗(yàn)的方法來決定。高速計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量模擬這樣的試驗(yàn)成為可能?？萍加?jì)算及社會中的問題比計(jì)算要復(fù)雜得多。比如金融

2、衍生產(chǎn)品（期權(quán)、期貨、掉期等）的定價(jià)及交易風(fēng)險(xiǎn)估算，問題的維數(shù)（即變量的個(gè)數(shù)）可能高達(dá)數(shù)百甚至數(shù)千。對這類問題，難度隨維數(shù)的增加呈指數(shù)增長，這就是所謂的“維數(shù)的災(zāi)難”(Curse of Dimensionality)，傳統(tǒng)的數(shù)值方法難以對付（即使使用速度最快的計(jì)算機(jī)）。Monte Carlo方法能很好地用來對付維數(shù)的災(zāi)難，因?yàn)樵摲椒ǖ挠?jì)算復(fù)雜性不再依賴于維數(shù)。以前那些本來是無法計(jì)算的問題現(xiàn)在也能夠計(jì)算量。為提高方法的效率，科學(xué)家們提出了許多所謂的“方差縮減”技巧。另一類形式與Monte Carlo方法相似，但理論基礎(chǔ)不同的方法“擬蒙特卡羅方法”(Quasi-Monte Carlo方法)近年來也

3、獲得迅速發(fā)展。我國數(shù)學(xué)家華羅庚、王元提出的“華王”方法即是其中的一例。這種方法的基本思想是“用確定性的超均勻分布序列(數(shù)學(xué)上稱為Low Discrepancy Sequences)代替Monte Carlo方法中的隨機(jī)數(shù)序列。對某些問題該方法的實(shí)際速度一般可比Monte Carlo方法提出高數(shù)百倍，并可計(jì)算精確度。http:/ call random_seed call random_number(a)Matlab: x=rand(N) 產(chǎn)生元素在(0, 1)間隨機(jī)分布的N*N矩陣 s= rand(state,0) 重設(shè)該生成函數(shù)到初始狀態(tài)計(jì)算程序產(chǎn)生的隨機(jī)數(shù)不是真正的隨機(jī)數(shù)，它們是確定的，但

4、看上去是隨機(jī)的，且能通過一些隨機(jī)性的檢驗(yàn)，故常稱為偽隨機(jī)數(shù)。如對32位字長的計(jì)算機(jī)real procedure random(xi)integer array (li)nreal array (xi)nl0 any integer that 1 l0 231-1for i=1 to n doli =(231-1)除以16807 li-1的余數(shù)xi li/ (231-1)endfor其它算法1. 取初值x0 (0,1)，let xi be the fractional part of (+ xi-1)5 2. Let u0 =1. For i=1,2,3, let ui be the remai

5、nder of (8t-3) ui-1 divided by 28, and xi = ui/2i. Here t can be any large integer注意：上述隨機(jī)數(shù)序列均具周期性，如上頁random子程序的周期約230。Maple has a collection of random number generators.如：With(stats)x :=uniform(0.1):seq(x(),i=1.10);Matlab: x=rand(10,1)產(chǎn)生10個(gè)隨機(jī)數(shù)。如 (b-a)r+a x (a,b) integer(n+1)r) x 0,1,2,n試產(chǎn)生1000個(gè)在橢圓x2

6、+4y2=4內(nèi)均勻分布的隨機(jī)點(diǎn)：方法：在 中均勻產(chǎn)生足夠 多隨機(jī)點(diǎn)，當(dāng)位于橢圓內(nèi)的點(diǎn)數(shù)為1000時(shí)停止?？捎呻S機(jī)數(shù)序列 r (0,1)構(gòu)造一系列隨機(jī)數(shù)序列1122yx，xy有時(shí)隨機(jī)數(shù)產(chǎn)生函數(shù)通不過嚴(yán)格的隨機(jī)性測試。而在一些計(jì)算（如MC積分）中隨機(jī)性非常重要。故使用更長字節(jié)的計(jì)算機(jī)更好。應(yīng)用之一估計(jì)面積和體積Numerical integration選取(0, 1)中隨機(jī)數(shù)序列x1， x2， x3， xn。則niixfnxxf110)(1d)(誤差約 ，它并不能和一些高級的數(shù)值積分算法比擬，但對多維情況，MC方法卻很有吸引力。n1 niiiizyxfnzyxzyxf1101010),(1ddd)

7、,(我們可產(chǎn)生一系列隨機(jī)數(shù)可簡單取3個(gè)隨機(jī)數(shù)構(gòu)成一個(gè)隨機(jī)點(diǎn)，即,.,654321),.,(),(654321相應(yīng)地，niibaxfnxxfab1)(1d)(1 niiidcbayxfnyxyxfabcd1),(1dd),()(1一般地，A) in points random nover of (average) of measure(d)(fAxfA例如： 求其中25. 0)5 . 0()5 . 0( : ),(22yxyxyxyxdd) 1ln(sin在該區(qū)域中取5000個(gè)點(diǎn)，可得該積分約0.57.計(jì)算體積Pseudo-code1sin10 10 102yezxzyxzyxProgram v

8、olume_regionInteger parameter n 5000, iprt 1000Real array integer i, mReal vol,x,y,zCall random(rij)For I=1 to n dox ri1y ri2z ri3If then m m+1If mod(I,iprt)=0 thenVol real(m)/real(I)Output I, volEndifEndforEnd program3)(nijr1 sin2yezxzyx計(jì)算體積之作業(yè)：冰淇淋的體積2222221) 1(yxzzyxyxz應(yīng)用之二MC模擬生日問題生日問題假設(shè)有假設(shè)有n個(gè)人在一起

9、，各自的生日為個(gè)人在一起，各自的生日為365天之一，根據(jù)概率理論，與很多人的直天之一，根據(jù)概率理論，與很多人的直覺相反，只需覺相反，只需23個(gè)人便有大于個(gè)人便有大于50的幾的幾率人群中至少有率人群中至少有2個(gè)人生日相同。個(gè)人生日相同。365)1(365.3653633653643653651nn 理論幾率 模擬幾率100.117 0.110200.411 0.412230.507 0.520300.706 0.692450.941 0.93650 0.986 0.987Real function prob(n,npts)Integer i, nptsLogical birthReal sumS

10、um=0For I=1 to npts doIf birth(n) then sum sum+1EndforProb sum/real(npts)End function蒲豐的投針問題蒲豐的投針問題Buffons Needle problem蒲豐證明在相距蒲豐證明在相距d 的一系列平行線上的一系列平行線上投長為投長為 l 針。針與線的交點(diǎn)數(shù)除以投針。針與線的交點(diǎn)數(shù)除以投的次數(shù)等于的次數(shù)等于2 l / (d)。Real function prob(n,npts)Integer i, nptsLogical birthReal sumSum=0For I=1 to npts doIf birth(

11、n) then sum sum+1EndforProb sum/real(npts)End function不妨取不妨取d l。產(chǎn)生隨機(jī)數(shù)。產(chǎn)生隨機(jī)數(shù)u和和 ，其中，其中u為針中心距最近的線的距離，故為針中心距最近的線的距離，故u小小于等于于等于0.5， 為夾角，根據(jù)對稱性，為夾角，根據(jù)對稱性，可取可取0到到/2。usin21即為相交usin21算時(shí)需要用上是否不自洽（姜冰）作業(yè)中子屏蔽問題中子屏蔽問題Neutron Shielding problem假設(shè)假設(shè)鉛墻長為5，中子在鉛中的平均自由程為1，中子與鉛原子碰撞，中子與鉛原子碰撞后各向同性散射。令碰撞后各向同性散射。令碰撞8次后中次后中子能

12、量耗盡，試求穿透鉛墻的中子能量耗盡，試求穿透鉛墻的中子的比例。子的比例。暫不考慮垂直紙面的運(yùn)動，則中暫不考慮垂直紙面的運(yùn)動，則中子的水平位移是子的水平位移是 。1721cos.coscos1入口鉛墻（長為5）221點(diǎn)問題點(diǎn)問題Blackjack 21 problem參數(shù)擬合問題參數(shù)擬合問題Parameter fitting problem輻射轉(zhuǎn)移問題輻射轉(zhuǎn)移問題Radiative transfer problem1 Ahn,S.-H., Lee, H.-W. & Lee, H.-M. 2000, JKAS, 33, p29.2 Zheng Zheng & Jordi Miral

13、da-Escude. 2002, ApJ, 578, p33.3Watson, A.M. & Henney, W.J. 2001, RMxAA, 37, p221.4Lorne W. Avery & Lewis L. H. , 1968, ApJ, 152, p493.5Lawrence J. C., Peter D. N. & Jeffery D.S. 1972, ApJ, 176, 439.6David A. Neufeld, 1990, ApJ, 350, p216.輻射轉(zhuǎn)移問題對三維輻射轉(zhuǎn)移問題,蒙特卡羅模擬幾乎是必須的計(jì)算方法1. 將空間劃分為Nx*Ny*Nz個(gè)網(wǎng)格2. 選取空間中的隨機(jī)點(diǎn)產(chǎn)生光子包(photon package), 并此向隨機(jī)方向傳播3. 考慮光子與粒子的隨機(jī)碰撞產(chǎn)生具有某種分布的隨機(jī)數(shù))(xfy dxxfxF)()()(1xfFdxdxFdF即正比于即分布密度試驗(yàn)粒子數(shù)值模擬研究粒子在電磁場中的運(yùn)動研究粒子加速等全粒子試驗(yàn)粒