蒙特卡羅隨機模擬投點法

蒙特卡羅隨機模擬投點法在數(shù)字積分中的使用數(shù)學(xué)和使用數(shù)學(xué)0901班：張瑞宸指導(dǎo)老師：任明慧摘要：本文首先介紹了蒙特卡羅方法的產(chǎn)生和發(fā)展，然后分析了蒙特卡羅方法計算數(shù)值積分的理論原理，最后給出了蒙特卡羅方法計算數(shù)值積分的MATLAB編程實現(xiàn)，全文主要是討論了蒙特卡羅方法在定積分計算的使用。而蒙特卡羅的優(yōu)點：可以計算被積函數(shù)非常復(fù)雜的定積分、重積分，并且維數(shù)沒有限制，這是別的數(shù)值積分方法還未達到的。蒙特卡羅的缺點：收斂速度慢，誤差一般較大，且是概率的誤差，不是真正的誤差。關(guān)鍵詞：蒙特卡羅方法，均值估計法，數(shù)值積分， Matlab編程Abstract：ThispaperfirstintroducestheemergenceanddevelopmentoftheMonteCarlomethod,andthenanalyzethetheoreticalprinciplesofMonteCarlonumericalintegrationmethod,Full-textmainlydiscussedtheapplicationoftheMonteCarlomethodinthedefiniteintegral.TheadvantagesofMonteCarlo:canbecalculatedtheintegrablefunctionsverycomplexdefiniteintegral,Multipleintegrals,anddimensionnolimit,othernumericalintegrationmethodshavenotyetreached.MonteCarloDisadvantages:slowconvergencespeed,errorgenerallyhigher,andtheprobabilityoferror,notarealerror.Keywords:MonteCarlomethod，Meanestimationmethod，numericalintegral，Matlabprogramming0引言歷史上有記載的蒙特卡羅試驗始于十八世紀末期（約1777年），當時布豐（Buffon）為了計算圓周率，設(shè)計了一個“投針試驗”，后文會給出。雖然方法已經(jīng)存在了200多年，此方法命名為蒙特卡羅則是在二十世紀四十年，美國原子彈計劃的一個子項目需要使用蒙特卡羅方法模擬中子對某種特殊材料的穿透作用。出于保密緣故，每個項目都要一個代號，傳聞命名代號時，項目負責(zé)人之一vonNeumann靈犀一點選擇摩洛哥著名賭城蒙特卡羅（MonteCarlo）作為該項目名稱，自此這種方法也就被命名為MonteCarlo方法廣為流傳。

蒙特卡羅方法，又名隨機模擬法或統(tǒng)計實驗法它是以概率統(tǒng)計理論為基礎(chǔ)，依據(jù)大數(shù)定律(樣本均值替代總體均值)利用電子計算機數(shù)字模擬技術(shù)，解決一些很難直接用數(shù)學(xué)運算求解或用其他方法不能解決的復(fù)雜問題的一種近似計算法。本世紀40年代電子計算機的出現(xiàn)，特別是近年來高速電子計算機的出現(xiàn)，使得用數(shù)學(xué)方法在計算機上大量、快速地模擬這樣的試驗成為可能。通常蒙特卡羅方法通過構(gòu)造符合一定規(guī)則的隨機數(shù)來解決數(shù)學(xué)上的各種問題。對于那些由于計算過于復(fù)雜而難以得到分析解或者根本沒有分析解的問題，蒙特卡羅方法是一種有效的求出數(shù)值解的方法。一般蒙特卡羅法在數(shù)學(xué)中最常見的使用就是蒙特卡羅隨機投點和蒙特卡羅數(shù)值積分。1蒙特卡羅方法的產(chǎn)生和發(fā)展蒙特卡羅方法是在二戰(zhàn)期間產(chǎn)生和發(fā)展起來的他的奠基者是美籍匈牙利人數(shù)學(xué)家馮諾伊曼(J.VonNeumann1903-1957)由于通常計算量相當大而電子計算機在當時還沒有出現(xiàn)，所有運算只能用手工進行，故而相當長的時間里蒙特卡羅方法難以推廣。1.1蒙特卡羅方法的產(chǎn)生作為蒙特卡羅方法的最初使用，是解決蒲豐氏問題1777年，法國數(shù)學(xué)家Buffon提出利用投針實驗求解的問題:設(shè)平面上有無數(shù)多條距離為1的等距平行線，現(xiàn)向該平面隨機的投擲一根長度為Z(/<1)的針。隨機投針是指針的中心點于最近的平行線間的距離x均勻分布在［0,1/2］上，針和平行線的夾角中(不管相交和否)均勻分布在【0,兀］上，如圖一所示，故而，故其x?U(0,1)，中~U(0,兀)概率密度函數(shù)分別為2,、1 ,小、1相交的充P(x)=-,p(里)=—相交的充們得到針和線要條件是x<們得到針和線要條件是x<lsin中2圖1.1針和線相交的幾種情況

則針和線相交的概率是p(x<—sin中)2="J—sin9p(x)p(中)="J—sin92dxd^TOC\o"1-5"\h\z00 00丸J2—sin9d9=22兀2 兀一，一，一一、 2/所以得到圓周率兀=——2——\o"CurrentDocument"p(x<+sin9) p假如我們能做大量的投針實驗并記錄下針和線的相交次數(shù)，則可以根據(jù)大數(shù)定律估計出針線相交的概率P。投針實驗N次可能有n次使針和任意平行線相交,那么p--，顯然，試驗次數(shù)N越多，P的近似程度越好。歷史上曾有幾位學(xué)者N做過這樣的投針試驗，并用手工計算出兀值，結(jié)果參見表1。表1實驗者時間(年)針長/投針次數(shù)相交次數(shù)兀估計值Wolf18500.80500025323.15956Smith18550.60302412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592922一一應(yīng)該指出，上述試驗的精度一般不會很高。譬如，假設(shè)/=1,則p=-。則由兀中心極限定理，如果試驗次數(shù)為N，則p的估計值p漸進服從N(p,四二0)，近N似為N(0.6366,02313)。因此，若要以95%的概率保證p的精確到三位有效數(shù)字，N即p和p的差距小于0.001，則N必須滿足N>1.962x0.2313/0.0012"8.89x105。重復(fù)進行上千次的投針實驗和手工計算，要消耗大量的人力、財力，蒙特卡羅方法雖然能解決此類問題，但得不到推廣使用。1.2蒙特卡羅方法的發(fā)展20世紀40年代以后，隨著電子計算機的出現(xiàn)和發(fā)展，人們有可能用計算機來模擬這類實驗和計算。計算機具有計算速度高和存儲容量大的特點，采用數(shù)字模擬技術(shù)可以代替許多實際上非常龐大而復(fù)雜的實驗，并迅速將實驗結(jié)果進行運算

處理，于是MonteCarlo方法重新被提起，引起世人重視，使用日漸廣泛。實際上，采用MonteCarlo方法在計算機上建立模型來解決Buffon問題是非常簡單的。我們在計算機上進行模擬試驗，給定，，我們可以在計算機上隨機產(chǎn)生x和'，然后判斷，＜-是否成立。若成立，則針線相交，否則不交。假如我們在sin中2計算機上獨立的產(chǎn)生N對這樣的x和中，并記錄下，＜L成立的次數(shù)，記為n。sin中 2則兀估計值可取為"，這就是隨機模擬計算的結(jié)果。借助Matlab在計算機上進行投針實驗，（相應(yīng)程序見附錄1）所取得的計算結(jié)果參見表2。表2N^-^nl=0.4l=0.6l=0.8l=0.9l=1.0500003.12524423.15839343.13922463.146423093.14643511000003.12414573.12119783.12695433.13245043.12866642000003.14144353.13963533.13734723.13662623.13718115000003.15193923.14707873.14554433.14228263.14254928000003.14833583.14656093.14408253.14230043.14302943.14352963.14590723.14126463.14193703.14233683.14971973.14629073.14187223.14105703.14073063.13933303.13987353.14107843.14168363.1412106其中兀=3.1416。不難看出，利用計算機運算不僅結(jié)果精確，而且迅速。隨著電子計算機的普及，蒙特卡羅方法作為一種獨到的方法得到開發(fā)，并首先使用到核武器的試驗和研制中。尤其是各種可是編程方法的不斷涌現(xiàn)，更顯示出蒙特卡羅方法的最獨到的優(yōu)點，即形象直觀地用數(shù)學(xué)方法在電子計算機上實現(xiàn)數(shù)字模擬實驗。2 蒙特卡羅方法計算數(shù)值積分的理論原理

隨機模擬是一種隨機試驗的方法，也稱為蒙特卡羅方法（MonteCarlo）［1］。這種方法利用隨機試驗，根據(jù)頻率和概率、平均值和期望值等之間的關(guān)系，推斷出預(yù)期的結(jié)果。2.1蒙特卡羅算法的理論原理隨機模擬法的基本原理非常簡單，下面給以直觀的說明。首先構(gòu)造一個定積分4j\;1-x2dx=兀，其中積分S=j\/1-x2dx（即1/4單0圖2.100圖2.1位圓的面積）用隨機模擬方法求得。圖2.1中粗線是1/4單位圓，如果向圖2.1中邊長為1的正方形里隨機投n塊小石頭，當n很大時小石頭會大致均勻地分布在正方形中，數(shù)一下落在1/4單位圓內(nèi)的小石頭，假定有k個，那么k/n就可以看作1/4單位圓面積n/4的近似值。小石頭的位置坐標可以用產(chǎn)生均勻分布隨機數(shù)的程序得到，記為x,J（i=1,2,...,n）是［0,1］區(qū)間均勻分布隨機數(shù)（x,y相互獨立），記錄滿足ii iix2+y2<1（i=1,2,...,n）的數(shù)量k，即得n=4k/n。我們用概率論中的大數(shù)定律來說明這個直觀認識的原理。大數(shù)定律（伯努利（Bernoull）定理）⑵設(shè)k是n次獨立重復(fù)試驗中事件A發(fā)生的次數(shù)，p是事件A每次試驗中發(fā)生的概率，則對任意的正數(shù)￡，有TOC\o"1-5"\h\z「 "k 〕 「limP〈一-p<￡>=1 （1）nT8 〔|n J若規(guī)定“向圖中正方形隨機投一塊小石頭落在四分之一單位圓里”為事件A發(fā)生，則A發(fā)生的概率p應(yīng)該等于四分之一單位圓面積，隨機投n塊石頭就是獨立重復(fù)做n次試驗，事件A發(fā)生k次，由（1）式，n無限變大時k/n和p之差小于任意一個數(shù)￡的概率趨于1。用計算機算一下就會發(fā)現(xiàn)，即使n很大結(jié)果也不好，并且很不穩(wěn)定，遠不

如常規(guī)數(shù)值積分的幾種方法。實際上，隨機模擬法很少用來做定積分，它的特點是能夠方便地推廣到計算多重積分，而不少多重積分是其他方法很難或者根本無法計算的。2.2蒙特卡羅方法計算定積分的理論原理通常有兩類辦法計算定積分。(1)隨機投點法在“投石算面積”的例子中，事件A在每次試驗中發(fā)生的概率p是四分之一單位圓面積，即(2)p=f1ji2dydx=f\1-x2dx(2)00 0n次試驗由計算機完成，采用［0,1］區(qū)間上的均勻分布產(chǎn)生相互獨立的隨機數(shù) 記這文樣產(chǎn)生的n個點的坐標為3，y=L2，...,n) 事件A發(fā)生等價于3,y)^X。記這樣樣產(chǎn)生的 |點的坐標、/為 ii 。事件A^發(fā)生等價I^T ii滿足y.V<1-X2,A發(fā)生的次數(shù)是滿足y.VJ1-X2(i=1,2,...,n)點的個數(shù)k。由伯努利定理，p可以用k/n近似替代。這種方法可以推廣如下。對任意區(qū)間 ［a,b］內(nèi)的連續(xù)函數(shù)f(x)，滿足c<f(x)<d(c,d>0)，為計算定積分jbf(x)dx，先由計算機產(chǎn)生n個點的坐標a(x,y)(i=1,2,...,n)，其中x,y分別為［a,b］和［c,d］區(qū)間上的均勻分布隨機數(shù)，然ii ii后記錄n個點中滿足y,<f(x.)的數(shù)目k，則(3)k(3)jbf(x)dxr(b一a)(d一c)—+(b一a)ca n(2)均值估計法這種方法依據(jù)概率論的一下兩個定理:大數(shù)定理(辛欽定理)［2］若隨機變量［,Y2,...,七相互獨立，服從同一個分布，且具有數(shù)學(xué)期望EY=r(i=1,2,...,n),則對任意的正數(shù)e，有i(4)limP<J_Z"Y-日<8>=1isHn.(4)I=1

隨機變量函數(shù)的期望［2］若隨機變量X的概率分布密度是p(x)(a<x<b)，則隨機變量的函數(shù)Y=f(X)的期望為E(f(X))=1bf(x)p(x)dx (5)a于是，當X為［a,b］區(qū)間均勻分布的隨機變量時，p(x)=1/(b一a)(a<x<b)，(5)式給出(6)jbf(x)dx=(b-a)E(f((6)a只要產(chǎn)生［a,b］區(qū)間相互獨立、均勻分布的隨機數(shù)x.(i=1,2,...,n)，y.=f(x)就相互獨立。根據(jù)辛欽定理，當n很大時期望EY=E(f(X))就可以用f(x)的平均值i近似，所以由(6)式得到(7)jbf(x)dx-—_a)Ef(x)(7)a ni=1.和隨機投點法相比，均值估計法沒有c<f(x)<d(c,d>0)的限制，只需計算f(x)，不需要產(chǎn)生隨機數(shù)y，也不需要作y<f(x)的比較，顯然大為方便。i i i i2.3 蒙特卡羅方法計算重積分的理論原理用隨機模擬作數(shù)值積分的優(yōu)點不僅在于計算簡單，尤其是它可以方便地推廣到計算多重積分，而不少多重積分是其他方法很難或者根本無法計算的。而通過上文對蒙特卡羅計算定積分的兩種方法的比較，我們可知均值估計法比投點法更好操作，經(jīng)過計算還可知均值估計法精確度(若選擇的密度函數(shù)適合)會更高。重積分本文討論均值估計法?？紤]k維積分：I=j...jf(x,x,...,x)dxdx...dx1 2k1 2kQ其中Q為k維積分區(qū)域，。匚V(V:a<x<b，i=1,2,...,k)。iiin卜的一個;IW灤宓聲屎新rt(YV OuV(V:a<x<b選取Q上的I概率密度函數(shù)p(x,x2,...,x^)，一(iii，i=1,2,...,k)，且p(x,x,...,x)滿足j…jp(x,x,...,x)dxdx...dx=1，1 2k 1 2 k1 2kQ

(X1,X2,...,x「GQ,故：p(x,x,...,x)dxdx...dxTOC\o"1-5"\h\z1 2k1 2kf(x,x，…,x)其中,(氣,x2,...,x「G。其中,(氣,x2,...,x「G其中,(氣,x2,...,x「G。其中,(氣,x2,...,x「G。(8)p(x,x,...,x)1 mf(x,x,...,x)\o"CurrentDocument"注—> k_m.]p(x,x,...,x)這里m是落在。區(qū)域中的隨機投點數(shù)。對于二重積分：I=jjf(x,y)dxdy其中。:a<x<b，c<g1(x)<y<g2(x)<d。記U=[a,b]x[c,d]，設(shè)S。，S^分1力別是。，V的測度，右選取g(x,y)是。上的均勻分布，則g(x,y)=一,(x,y)g。。S由幾何概率得以=m，即七r七，于是：SnmnI=jjf(x,y)g(x,y)dxdy。g(x,y)=11mf(x,y)m.=1g(x,y)=*'f(x,y)=土丈f(x,y)i=1 i=1綜上可得，1=斗丈f(x,y) (9)i=13蒙特卡羅方法MATLAB編程實例用隨機模擬作數(shù)值積分主要利用均勻分布產(chǎn)生的隨機數(shù)及相應(yīng)的判斷、計算和簡單運算。MATLAB提供的命令unifrnd(a,b,m,n)產(chǎn)生m行n列［a,b］區(qū)間上均勻分布隨機數(shù)。當a=0，b=1時，可用rand(m,n)。

彳列1用蒙特卡羅方法計算n。解(1)隨機投點法。MATLAB程序如下：n=10000;x=rand(2,n);k=0;fori=1:nifx(1,i)A2+x(2,i)A2<=1k=k+1;endendp=4*k/n重復(fù)計算4次，計算結(jié)果分別為：p=3.1244 P=3.1420 P=3.1780 P=3.1592當n提高到50000時，重復(fù)計算4次，計算結(jié)果為：p=3.1372 P=3.1456 P=3.1442 P=3.1426解(2)均值估計法，此時f(x)=x：1-x2,MATLAB程序如下：n=50000;x=rand(1,n);y=0;fori=1:ny=y+sqrt(1-x(i)八2);endp=4*y/n重復(fù)計算4次，計算結(jié)果分別為：p=3.1439 P=3.1510 P=3.1474 P=3.1449當n提高到50000時，重復(fù)計算4次，計算結(jié)果為：p=3.1471 P=3.1419 P=3.1423 P=3.1430可以看出它的缺點是計算量大，結(jié)果具有隨機性，精度較低。一般說來精度為n-1/2階，n增加時精度提高較慢。但是用隨機模擬方法可以計算被積函數(shù)非常復(fù)雜的定積分、重積分，并且維數(shù)沒有限制。如下介紹多重積分的MALTAB實現(xiàn)。3.1多重積分的蒙特卡羅方法MATLAB編程實現(xiàn)以二重積分為例，運用算法(8)編程，一個是通常的for循環(huán)編程，文件名為mtc.m，另一個是向量化編程，文件名為mtcx.m，后者運行速度要快很多。functions=mtc(f,fai1,fai2,a,b,c,d,n)ifnargin<8n=10000;endx=unifrnd(a,b,1,n);y=unifrnd(c,d,1,n);s=0;fork=1:n

ifand(y(k)>=feval(fai1,x(k)),y(k)<=feval(fai2,x(k)))s=s+feval(f,x(k),y(k));endends=s/n*(b-a)*(d-c);returnfunctions=mtcx(f,fai1,fai2,a,b,c,d,n)ifnargin<8n=10000;endx=unifrnd(a,b,1,n);y=unifrnd(c,d,1,n);s=sum(feval(f,x,y)).*and(y>=feval(fai1,x),y<=feval(fai2,x)));s=s/n*(b-a)*(d-c);return三重或三重以上積分只須增加若十個輸入?yún)?shù)，即新增積分變量的上、下限，方法不變，對程序稍作修改即可。例2計算A』'匕e-：sin(x2+y)dydx- X22-"一2積分的精確值為0.41192954617630，在命令窗口調(diào)用函數(shù)M文件mtc.m和mtcx.m情況如下：>>tic;s1=mtc(inline(‘exp(-x.A2/2).*sin(x.A2+y)5),inline(‘-sqrt(1-x.A2/2),),inline(‘sqrt(1-x.A2/2),),-0.5,1,-1,1,100000),tocs1=0.40901034347183Elapsedtimeis40.094000seconds.>>tic;s2=mtcx(inline(‘exp(-x.A2/2).*sin(x.A2+y)’),inline(‘-sqrt(1-x.A2/2)’),inline(‘sqrt(1-x.A2/2)’),-0.5,1,-1,1,100000),tocs2=0.41183962563018Elapsedtimeis0.047000seconds.在運用了向量化編程后運行時間從40.094s縮短到0.047s，是原來的1/853,

效果好；誤差為0.00008992054612000366,也比前者好。若將上面程序稍作如下修改，即輸入?yún)?shù)減少2個，那么當n=1000000時，運行時間大約增加0.45s，精度不變。functions=mtcxg(f,fai1,fai2,a,b,c,d,n)ifnargin<6n=10000;endx=unifrnd(a,b,1,n);c=min(feval(fai1,x));d=max(feval(fai2,x));y=unifrnd(c,d,1,n);s=sum(feval(f,x,y)).*and(y>=feval(fai1,x),y<=feval(fai2,x)));s=s/n*(b-a)*(d-c);return>>tic;s3=mtcxg(inline(‘exp(-x.A2/2).*sin(x.A2+y)5),inline(‘-sqrt(1-x.A2/2),),inline(‘sqrt(1-x.A2/2),),-0.5,1,-1,1,100000),tocs3=0.41200550284008Elapsedtimeis0.922000seconds.對比常規(guī)數(shù)值積分和蒙特卡羅方法數(shù)值積分，同樣數(shù)量的n值一一也就意味這幾乎相同的計算量一一常規(guī)數(shù)值積分結(jié)果的精確度要高于蒙特卡羅數(shù)值積分的結(jié)果。那么，我們?yōu)楹芜€需要用蒙特卡羅來算數(shù)值積分呢？答案的關(guān)鍵在于，常規(guī)數(shù)值積分的精度直接取決于每個維度上取點數(shù)量，維度增加了，但是每個維度上要取的點卻不能減少。在多重積分中，隨著被積函數(shù)維度增加，需要計算的函數(shù)值數(shù)量以指數(shù)速度遞增。例如在一重積分jaf(a>(x)dx中，只要沿著x軸取n個點；要達到相同大b小的精確度，在k重積分jj...jf(xx...xN(xx...x)dxdx...dx中，仍然需要在12k12n1 2k每個維度上取n個點，k個緯度的坐標相組合，共需要計算nXk個坐標對應(yīng)的f(x)函數(shù)值。取點越多，會占用計算機大量內(nèi)存，也需要更長運算時間，最終導(dǎo)致這種計算方法不可行！蒙特卡羅方法卻不同，不管是積分有多少重，取n個點計算的結(jié)果精確度都差不多。因此，即使在一重積分的情形下，蒙特卡羅方法的效率比不過常規(guī)數(shù)值

積分，但隨著積分維度增加，常規(guī)數(shù)值積分的速度呈指數(shù)下降，蒙特卡羅方法的效率卻基本不變。經(jīng)驗表明，當積分重數(shù)達到4重積分甚至更高時，蒙特卡羅方法將遠遠優(yōu)于常規(guī)數(shù)值積分方法。4總結(jié)歸納起來本論文主要論述了：蒙特卡羅算法在數(shù)值積分中的運用。本篇論文給出了蒙特卡羅方法使用在數(shù)值積分中的原理，討論蒙特卡羅算法運用在數(shù)值積分中的兩種計算方法，隨機投點法和均值估計法，具體蒙特卡羅算法的實現(xiàn)給出了MATLAB編程的實例。最后對蒙特卡羅算法的優(yōu)缺點進行了討論。蒙特卡羅理論的使用非常廣泛，兩種蒙特卡羅數(shù)值計算方法的使用也很全面,由于作者的知識水平，研究能力有限，蒙特卡羅理論在各行各業(yè)中都有廣泛使用，思想、理論等方面仍存在欠缺之處，但這些都是以后的努力研究的方