一蒙特卡洛隨機(jī)模擬

系列一蒙特卡洛隨機(jī)模擬實(shí)驗(yàn)?zāi)康模簩W(xué)會(huì)用計(jì)算機(jī)隨機(jī)模擬方法來(lái)解決隨機(jī)性問(wèn)題蒙特卡洛模擬法簡(jiǎn)介蒙特卡洛(MonteCarlo)方法是一種應(yīng)用隨機(jī)數(shù)來(lái)進(jìn)行計(jì)算機(jī)摸擬的方法。此方法對(duì)研究對(duì)象進(jìn)行隨機(jī)抽樣，通過(guò)對(duì)樣本值的觀察統(tǒng)計(jì)，求得所研究系統(tǒng)的某些參數(shù)。作為隨機(jī)模擬方法，起源可追溯到18世紀(jì)下半葉蒲峰實(shí)驗(yàn)。蒙特卡洛模擬法的應(yīng)用領(lǐng)域蒙特卡洛模擬法的應(yīng)用領(lǐng)域主要有：直接應(yīng)用蒙特卡洛模擬：應(yīng)用大規(guī)模的隨機(jī)數(shù)列來(lái)模擬復(fù)雜系統(tǒng)，得到某些參數(shù)或重要指標(biāo)。蒙特卡洛積分：利用隨機(jī)數(shù)列計(jì)算積分，維數(shù)越高，積分效率越高。蒙特卡洛模擬法求解步驟應(yīng)用此方法求解工程技術(shù)問(wèn)題可以分為兩類(lèi)：確定性問(wèn)題和隨機(jī)性問(wèn)題。解題步驟如下：1.根據(jù)提出的問(wèn)題構(gòu)造一個(gè)簡(jiǎn)單、適用的概率模型或隨機(jī)模型，使問(wèn)題的解對(duì)應(yīng)于該模型中隨機(jī)變量的某些特征(如概率、均值和方差等)，所構(gòu)造的模型在主要特征參量方面要與實(shí)際問(wèn)題或系統(tǒng)相一致2.根據(jù)模型中各個(gè)隨機(jī)變量的分布，在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，實(shí)現(xiàn)一次模擬過(guò)程所需的足夠數(shù)量的隨機(jī)數(shù)。通常先產(chǎn)生均勻分布的隨機(jī)數(shù)，然后生成服從某一分布的隨機(jī)數(shù)，方可進(jìn)行隨機(jī)模擬試驗(yàn)。根據(jù)概率模型的特點(diǎn)和隨機(jī)變量的分布特性，設(shè)計(jì)和選取合適的抽樣方法，并對(duì)每個(gè)隨機(jī)變量進(jìn)行抽樣(包括直接抽樣、分層抽樣、相關(guān)抽樣、重要抽樣等)。按照所建立的模型進(jìn)行仿真試驗(yàn)、計(jì)算，求出問(wèn)題的隨機(jī)解。統(tǒng)計(jì)分析模擬試驗(yàn)結(jié)果，給出問(wèn)題的概率解以及解的精度估計(jì)。在可靠性分析和設(shè)計(jì)中，用蒙特卡洛模擬法可以確定復(fù)雜隨機(jī)變量的概率分布和數(shù)字特征，可以通過(guò)隨機(jī)模擬估算系統(tǒng)和零件的可靠度，也可以模擬隨機(jī)過(guò)程、尋求系統(tǒng)最優(yōu)參數(shù)等。一.預(yù)備知識(shí)：隨機(jī)數(shù)的產(chǎn)生提示：均勻分布U(0,1)的隨機(jī)數(shù)可由C語(yǔ)言或Matlab自動(dòng)產(chǎn)生，在此基礎(chǔ)上可產(chǎn)生其他分布的隨機(jī)數(shù).逆變換法：設(shè)隨機(jī)變量U服從(0,1)上的均勻分布，則X=F-'(U)的分布函數(shù)為F(x)步驟：(1)產(chǎn)生U(0J)的隨機(jī)數(shù)U；②計(jì)算X=F-1(U),則X服從F(x)分布.問(wèn)題：練習(xí)用此方法產(chǎn)生常見(jiàn)分布隨機(jī)數(shù)例如“指數(shù)分布，均勻分布U(a,b)”.還有其它哪種常見(jiàn)分布的隨機(jī)數(shù)可用此方法方便產(chǎn)生？產(chǎn)生離散分布隨機(jī)數(shù)己知離散隨機(jī)變量X的概率分布:P(X=xk)=I\,(K=1,2…)，產(chǎn)生隨機(jī)變量X的隨機(jī)數(shù)可采用如下算法：將區(qū)間［0.1］依次分為長(zhǎng)度為Pi,p?,???的小區(qū)間L，L，???；產(chǎn)生［0,1］均勻分布隨機(jī)數(shù)R,若Rclk則令X=xk,重復(fù)(b),即得離散隨機(jī)變量X的隨機(jī)數(shù)序列.問(wèn)題：(1)下表給出了離散分布X的概率分布表，試產(chǎn)生100個(gè)隨機(jī)數(shù)X的概率分布表Xk13579Pk0.50.050.20.050.2用此方法給出100個(gè)二項(xiàng)分布B(20,0.1)的隨機(jī)數(shù)及10個(gè)泊松分布P(l)的隨機(jī)數(shù).正態(tài)分布的抽樣提示：設(shè)U］,U2是獨(dú)立同分布的U(0Q變量，令X］=(-21nU］)”2cos(2^u2)X2=(-21nU1)1/2sin(2MJ2)則X.與X,獨(dú)立，均服從標(biāo)準(zhǔn)正態(tài)分布.步驟：(1)由U(0J)獨(dú)立抽取Ui=g=U2(2)用(*)式計(jì)算^,X2.用此方法可同時(shí)產(chǎn)生兩個(gè)標(biāo)準(zhǔn)正忐分布的隨機(jī)數(shù)問(wèn)題：有關(guān)隨機(jī)數(shù)產(chǎn)生方法很多，查閱相關(guān)材料進(jìn)行系統(tǒng)總結(jié).二.隨機(jī)決策問(wèn)題某小販每天以一元的價(jià)格購(gòu)進(jìn)一種鮮花，賣(mài)出價(jià)為b元/束，當(dāng)天賣(mài)不出去的花全部損失，顧客一天內(nèi)對(duì)花的需求量是隨機(jī)變量，服從泊松分布，P(X=k)=e-4—,k=0,1,2,...,,其中常數(shù);I由多口銷(xiāo)傳量的平均值來(lái)估計(jì)，問(wèn)小販每天應(yīng)購(gòu)進(jìn)多少束鮮花?(準(zhǔn)則:期望收入，(①最局)問(wèn)題：(1)在給定b=1.25,2=50的值后，畫(huà)出目標(biāo)函數(shù)S(u)連線散點(diǎn)圖，觀察單調(diào)性，給出最優(yōu)決策U*：。選取其他的b,4,再觀察S(u)的單調(diào)性；用計(jì)算機(jī)模擬方法來(lái)求出最優(yōu)決策對(duì)固定的U,例如，u=40,對(duì)隨機(jī)變量X模擬100次，每次模擬得到一個(gè)收入，求出100個(gè)收入的平均值，即得到在決策u=40情況下的可能收入；對(duì)所有的可能的u,重復(fù)(3),從中找最大的，并與⑴的結(jié)果相比較.蒙特卡羅計(jì)算兀值思路：圖1表示一個(gè)內(nèi)接于正方形的圓的半徑R.圓的面積是刀R，,正方形的面積是(2R)2.圓和正方形的面枳的比值就是將上述比值乘以4,就能獲得〃值.(2R)-4

2R2R圖1例用蒙特卡羅方法計(jì)算力值一重定積分的蒙特卡羅算法問(wèn)題描述：假設(shè)函數(shù)f(x)在［a,b］內(nèi)有界連續(xù)，且f(x)ZO,求解定積分I=「f(x)dxJa為計(jì)算出其值，可構(gòu)造概率模型如下：取一個(gè)邊長(zhǎng)分別為b-a和c的矩形D,使曲邊梯形在矩形域之內(nèi)，如圖2,并在矩形內(nèi)隨機(jī)投點(diǎn)，假設(shè)隨機(jī)點(diǎn)均勻地落在整個(gè)矩形之內(nèi),則落在圖中灰色區(qū)域內(nèi)的隨機(jī)點(diǎn)數(shù)k與投點(diǎn)總數(shù)N之比k/N就近似地等于曲線下方面枳(即陰影面枳)與矩形面積之比，從而得出近似積分a)c京京e1.——求例由于廠子是非初等函數(shù)，我們很難求出其原函數(shù)，所以用牛頓■萊布尼茨公式無(wú)法求解，但可以運(yùn)用蒙特卡羅方法求出其近似值.將上述方法推廣到一般情況：假設(shè)函數(shù)f(x)在［a,b］內(nèi)有界連續(xù)，對(duì)于定枳分I=［f(x)dx，為計(jì)算出其值，可構(gòu)造如下?概率模型:取一個(gè)邊長(zhǎng)分別為b-a和c-d的矩形D,使曲線［a,b］段的值在矩形域之內(nèi)，如圖3,并在矩形內(nèi)隨機(jī)投點(diǎn)，假設(shè)隨機(jī)點(diǎn)均勻地落在整個(gè)矩形之內(nèi)，則落在圖中x軸上下灰色區(qū)域內(nèi)的隨機(jī)點(diǎn)數(shù)m與n的差與投點(diǎn)總數(shù)p之比(m-n)/P就近似地等于曲線上下方面枳之差(即陰影面積之差)與矩形面枳之比，從而得出近似積分I/%2(b-a)(c-d).圖3二重積分的蒙特卡羅算法問(wèn)題描述：實(shí)際計(jì)算中常常要遇到如j］f(x,y)dxdy的二重積分，發(fā)現(xiàn)被積函數(shù)的原函數(shù)往*D往很難求出，或者原函數(shù)根本就不是初等函數(shù)，對(duì)于這樣的重積分，蒙特卡羅方法也有成熟的計(jì)算方法.方法1:步驟：取一個(gè)包含D的矩形區(qū)域Q：a<x<b,c<y<d,面積A=(b-a)(d-c)；(無(wú),乂)，i=1,2,???,n,為Q上的均勻分布隨機(jī)數(shù)列，不妨設(shè)(&乂),(i= 為落Ak在D中的ii個(gè)隨機(jī)數(shù)，則11充分大時(shí)，有f(x,y)dxdyb—習(xí)f(*,y)d n1=1方法2：對(duì)二重積分1= 假設(shè)f(&y)為區(qū)域A上的有界函數(shù)，且f(x,y)20,幾何意義對(duì)應(yīng)的是以f(x,y)為曲面頂，A為底的曲頂柱體C的體枳.因此，用均勻隨機(jī)數(shù)計(jì)算二重枳分的蒙特卡羅方法基本思路為：假設(shè)曲頂柱體C包含在己知體積為Vd的幾何體D的內(nèi)部，在D內(nèi)產(chǎn)生N個(gè)均勻隨機(jī)點(diǎn)，統(tǒng)計(jì)出在C內(nèi)部的隨機(jī)點(diǎn)數(shù)目Nc,則I=^Nc..

例：計(jì)算JJ(1+小一x^-y2一jL+y^dxdy,其中A={(x,y)|x2+y2<1}A分析：該二重積分可以看作以l+小一妒—y2-jV+y?為頂?shù)那斨黧w的體積,此曲頂柱體在一個(gè)邊長(zhǎng)為2的立方體內(nèi)，用數(shù)學(xué)分析方法可計(jì)算出其精確值為兀.5用蒙特卡洛方法計(jì)算體積z>y/x2+y2,z<1+1—x2—y2冰激凌錐含于體積=8的六面體Q={(x,y,z)-1<x<l,-l<y<l,0<z<2)蒙特卡羅求一元函數(shù)的最值數(shù)學(xué)分析中求一元函數(shù)f(x)在閉區(qū)間［a,b］的最值的方法為：找出函數(shù)f(x)的所有駐點(diǎn)和不可導(dǎo)點(diǎn)，然后計(jì)算出其相應(yīng)函數(shù)值，并將其與區(qū)間端點(diǎn)函數(shù)值比較得出函數(shù)的最大值或最小值，但實(shí)際問(wèn)題中方程f(x)=0很難求出其根，所以駐點(diǎn)往往很難求出.若利用蒙特卡羅方法，在［a,b］上產(chǎn)生n個(gè)隨機(jī)數(shù)，計(jì)算出這些數(shù)的函數(shù)值并作比較,則可得出最值的近似值.例求f(x)=(l-x3)sin3x在［-盆,2勿］上的最大值-6-4 -2 0 2 4 6x(1-x°)sin(3x)-6-4 -2 0 2 4 6x(1-x°)sin(3x)150100<-v=(1-x3<-v=(1-x3)*sin(3*x)0-50-100圖4蒙特卡羅方法求一元函數(shù)的根根據(jù)數(shù)學(xué)分析相關(guān)知識(shí)，求一元函數(shù)f(x)=0的根有很多方法，如一般迭代法，牛頓切線法等等，這些方法討論根的收斂性，與初始迭代值相密切相關(guān)通常要給出一個(gè)合適的初始迭代值知是比較困難的，利用蒙特卡羅方法可以擺脫根的收斂性對(duì)初始值Xg的依賴(lài)性.具體方法如下：要解方程f(x)=0,XG［a,b］,其中函數(shù)f(x)為連續(xù)函數(shù)，為指定精度令Xd=m=a,X*=b,k=l,F_(1)=FO,步驟如下：(1)^k=k+l,在［X頃，擔(dān)］內(nèi)產(chǎn)生n個(gè)隨機(jī)數(shù)*(i=l,2,...n)，計(jì)算并比較出這n個(gè)隨機(jī)數(shù)的函數(shù)絕對(duì)值的最小值f(&J=得事|1'(而)|｝，業(yè)為i的某個(gè)取值，令Ff(k)=min｛Fxi林(k-l),f(&)｝⑵若寫(xiě)_?)<言成立，則終止計(jì)算，令。碩=跖，根就是捻ot；若Fx^(k)>^fiF_(k)=F_(k-l),則令k=k—1,轉(zhuǎn)至(1)；若E_(k)>S且％頗?)<岌心(卜1)，則令。皿=自，轉(zhuǎn)至(3).(3)令d=cIq/k,Xdami=fgjt—d，紈？=。碩+d，轉(zhuǎn)至(1)?說(shuō)明：a)F°是人為給定的一個(gè)很大的正數(shù)，do?b-afido>Ob)k=k+1表示重新賦值給k,使k的值增加1,對(duì)k=k-1同理.C)區(qū)間［Xm'XJ一定在定義域［a,b］之內(nèi).此迭代步驟能使函數(shù)值序列Fxmi(l>Fx >.Fxmk,最終使岌心時(shí)成立，得出函數(shù)f(x)=0達(dá)到精度要求的根板心例求f(x)=e-°—tanx+800=0在(0,乙)上的實(shí)根，(假定|f(x)|<10'\則認(rèn)為x為根)2分析：對(duì)上面的方程，如取初值05,0.7,0.9,1.1,1.3,1.5,用牛頓迭代公式在區(qū)間內(nèi)找不到根，若用蒙特卡羅方法，則不需要給定初值.蒙特卡羅方法解線性規(guī)劃用蒙特卡羅方法可解約束規(guī)劃問(wèn)題：minf(X),XeEngi(X)Z0,i=l,2,...,ms.t.( ,?■■a,Sbj,j=l,2,?..,n基本思想：在估計(jì)的區(qū)域｛(改，也，…，％Xj￡［aj,bj］,j=1,2,...n)｝內(nèi)隨機(jī)取若干個(gè)試驗(yàn)點(diǎn)，然后從試驗(yàn)點(diǎn)中找出可行點(diǎn)，再?gòu)目尚悬c(diǎn)中選擇出使函數(shù)值最小的點(diǎn).符號(hào)假設(shè)和說(shuō)明：設(shè)試驗(yàn)點(diǎn)的第j個(gè)分量Xj服從［apbj內(nèi)的均勻分布；P：試驗(yàn)點(diǎn)總數(shù)；MAX{P)：最大試驗(yàn)點(diǎn)總數(shù)；K：可行點(diǎn)總數(shù)：MAX{K}：最大可行點(diǎn)數(shù)；X*：迭代產(chǎn)生的最優(yōu)點(diǎn)；R：在［0,1］上的均勻隨機(jī)數(shù)：Q：迭代產(chǎn)生的最小值f(X*),其初始