一蒙特卡洛隨機模擬

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

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

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