數(shù)學(xué)建模專題三MonteCarlo模擬課件

數(shù)學(xué)建模專題三MonteCarlo模擬幻燈片本課件PPT僅供大家學(xué)習(xí)使用學(xué)習(xí)完請(qǐng)自行刪除，謝謝！本課件PPT僅供大家學(xué)習(xí)使用學(xué)習(xí)完請(qǐng)自行刪除，謝謝！內(nèi)容提綱1.引言2.MonteCarlo模擬根本思想3.隨機(jī)數(shù)生成函數(shù)4.應(yīng)用實(shí)例舉例5.排隊(duì)論模擬6.MonteCarlo模擬求解規(guī)劃問(wèn)題引言(Introduction)MonteCarlo方法：蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來(lái)的。亦稱統(tǒng)計(jì)模擬方法，statisticalsimulationmethod利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法MonteCarlo名字的由來(lái)：MonteCarlo是摩納哥（monaco)的首都，該城以賭博聞名NicholasMetropolis(1915-1999)Monte-Carlo,Monaco二十世紀(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è)簡(jiǎn)單實(shí)例MonteCarlo方法的根本思想舉例Buffon投針實(shí)驗(yàn)1768年，法國(guó)數(shù)學(xué)家ComtedeBuffon利用投針實(shí)驗(yàn)估計(jì)的值dLSolutionThepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:Therandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegral例1.蒲豐投針問(wèn)題利用關(guān)系式：求出π值其中Ｎ為投計(jì)次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問(wèn)題。一些人進(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根本思想由上面的例子可以看出，當(dāng)所求問(wèn)題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的數(shù)學(xué)期望，或者是與之有關(guān)的量時(shí)，通過(guò)某種試驗(yàn)的方法，得出該事件發(fā)生的頻率，再通過(guò)它得到問(wèn)題的解。這就是蒙特卡羅方法的根本思想。蒙特卡羅方法的基本思想雖然早已被人們提出，卻很少被使用。直到電子計(jì)算機(jī)出現(xiàn)后，使得人們可以通過(guò)電子計(jì)算機(jī)來(lái)模擬巨大數(shù)目的隨機(jī)試驗(yàn)過(guò)程，使得蒙特卡羅方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中發(fā)揮應(yīng)有的作用。例1在我方某前沿防守地域，敵人以一個(gè)炮排（含兩門火炮）為單位對(duì)我方進(jìn)行干擾和破壞．為躲避我方打擊，敵方對(duì)其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn)．經(jīng)過(guò)長(zhǎng)期觀察發(fā)現(xiàn)，我方指揮所對(duì)敵方目標(biāo)的指示有50％是準(zhǔn)確的，而我方火力單位，在指示正確時(shí)，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部毀傷敵人火炮．現(xiàn)在希望能用某種方式把我方將要對(duì)敵人實(shí)施的20次打擊結(jié)果顯現(xiàn)出來(lái)，確定有效射擊的比率及毀傷敵方火炮的平均值。分析：這是一個(gè)概率問(wèn)題，可以通過(guò)理論計(jì)算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動(dòng)的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動(dòng)的動(dòng)態(tài)過(guò)程.為了能顯示我方20次射擊的過(guò)程，現(xiàn)采用模擬的方式。舉例需要模擬出以下兩件事：

[2]當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況[1]觀察所對(duì)目標(biāo)的指示正確與否模擬試驗(yàn)有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2．因此，可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時(shí)為指示正確，反之為不正確．模擬試驗(yàn)有三種結(jié)果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒(méi)能毀傷敵火炮的可能性為1/2(即3/6)．這時(shí)可用投擲骰子的方法來(lái)確定：如果出現(xiàn)的是１、２、３三個(gè)點(diǎn)：則認(rèn)為沒(méi)能擊中敵人；如果出現(xiàn)的是４、５點(diǎn)：則認(rèn)為毀傷敵人一門火炮；若出現(xiàn)的是６點(diǎn)：則認(rèn)為毀傷敵人兩門火炮．問(wèn)題分析i：要模擬的打擊次數(shù)；k1：沒(méi)擊中敵人火炮的射擊總數(shù)；k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)；E：有效射擊比率；E1：20次射擊平均每次毀傷敵人的火炮數(shù)．符號(hào)說(shuō)明模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點(diǎn)數(shù)?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=0*k1/20+1*k2/20+2*k3/20停止硬幣正面?YNNY1，2，34，56模擬結(jié)果理論計(jì)算結(jié)果比較

雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過(guò)程．用蒙特卡洛方法進(jìn)行計(jì)算機(jī)模擬的步驟：[1]設(shè)計(jì)一個(gè)邏輯框圖，即模擬模型．[2]根據(jù)流程圖編寫程序，模擬隨機(jī)現(xiàn)象．可通過(guò)具有各種概率分布的模擬隨機(jī)數(shù)來(lái)模擬隨機(jī)現(xiàn)象．[3]分析模擬結(jié)果，計(jì)算所需要結(jié)果.注：rand(n)=rand(n,n)Matlab中的隨機(jī)數(shù)生成函數(shù)randperm(m)生成一個(gè)由1:m

組成的隨機(jī)排列randn(m,n)生成一個(gè)滿足正態(tài)m

n

的隨機(jī)矩陣rand(m,n)

生成一個(gè)滿足均勻分布的m

n

隨機(jī)矩陣，矩陣的每個(gè)元素都在(0,1)

之間。perms(1:n)生成由1:n

組成的全排列，共n!

個(gè)name

的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......random('name',A1,A2,A3,M,N)Matlab中的隨機(jī)數(shù)生成函數(shù)fix(x):

截尾取整，直接將小數(shù)部分舍去floor(x):

不超過(guò)x

的最大整數(shù)ceil(x):

不小于x

的最小整數(shù)round(x):

四舍五入取整Matlab中的取整函數(shù)

隨機(jī)投擲均勻硬幣，驗(yàn)證國(guó)徽朝上與朝下的概率是否都是1/2

n=10000;%給定試驗(yàn)次數(shù)m=0;fori=1:nx=randperm(2)-1;y=x(1);ify==0%0表示國(guó)徽朝上，1表示國(guó)徽朝下m=m+1;endendfprintf('國(guó)徽朝上的頻率為：%f\n',m/n);小實(shí)例一：投擲硬幣隨機(jī)投擲骰子，驗(yàn)證各點(diǎn)出現(xiàn)的概率是否為1/6

n=10000;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;fori=1:nx=randperm(6);y=x(1);switchycase1,m1=m1+1;case2,m2=m2+1;case3,m3=m3+1;case4,m4=m4+1;case5,m5=m5+1;otherwise,m6=m6+1;endend...%

輸出結(jié)果小實(shí)例二：投擲骰子

用蒙特卡羅(MonteCarlo)投點(diǎn)法計(jì)算

的值n=100000;a=2;m=0;fori=1:nx=rand(1)*a/2;y=rand(1)*a/2;if(x^2+y^2<=(a/2)^2)m=m+1;endendfprintf('計(jì)算出來(lái)的pi為：%f\n',4*m/n);小實(shí)例三：蒙特卡羅投點(diǎn)法

在畫有許多間距為d

的等距平行線的白紙上，隨機(jī)投擲一根長(zhǎng)為l(l

d)的均勻直針，求針與平行線相交的概率，并計(jì)算的值小實(shí)例四：蒲豐投針實(shí)驗(yàn)n=100000;l=0.5;d=1;m=0;fori=1:nalpha=rand(1)*pi;y=rand(1)*d/2;ify<=l/2*sin(alpha)m=m+1;endendfprintf('針與平行線相交的頻率為：%f\n',m/n);fprintf('計(jì)算出來(lái)的pi為：%f\n’,2*n*l/(m*d));小實(shí)例四源程序

設(shè)某班有m

個(gè)學(xué)生，則該班至少有兩人同一天生日的概率是多少？小實(shí)例五：生日問(wèn)題解：設(shè)一年為365天，且某一個(gè)學(xué)生的生日出現(xiàn)在一年中的每一天都是等可能的，則班上任意兩個(gè)學(xué)生的生日都不相同的概率為：所以，至少有兩個(gè)學(xué)生同一天生日的概率為：n=1000;p=0;m=50;%設(shè)該班的人數(shù)為50fort=1:na=[];q=0;fork=1:mb=randperm(365);a=[a,b(1)];endc=unique(a);iflength(a)~=length(c)p=p+1;endendfprintf(‘至少兩人同一天生日的概率為：%f\n',p/n);試驗(yàn)五源程序clear;m=50;p1=1:365;p2=[1:365-m,365*ones(1,m)];p=p1./p2;p=1-prod(p);fprintf('至少兩人同一天生日的概率為：%f\n',p);小實(shí)例五的理論值計(jì)算排隊(duì)問(wèn)題隨機(jī)模擬排隊(duì)論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過(guò)程。在排隊(duì)系統(tǒng)中，服務(wù)對(duì)象的到達(dá)時(shí)間和服務(wù)時(shí)間都是隨機(jī)的。排隊(duì)論通過(guò)對(duì)隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計(jì)研究，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律指標(biāo)，如排隊(duì)的長(zhǎng)度、等待的時(shí)間及服務(wù)利用率。[1]系統(tǒng)的假設(shè)：（1）顧客源是無(wú)窮的；（2）排隊(duì)的長(zhǎng)度沒(méi)有限制；（3）到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)，“先到先服務(wù)”。在某商店有一個(gè)售貨員，顧客陸續(xù)來(lái)到，售貨員逐個(gè)地接待顧客．當(dāng)?shù)絹?lái)的顧客較多時(shí)，一部分顧客便須排隊(duì)等待，被接待后的顧客便離開商店．設(shè)：

1．顧客到來(lái)間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布．２．對(duì)顧客的服務(wù)時(shí)間服從［４，１５］上的均勻分布．３．排隊(duì)按先到先服務(wù)規(guī)則，隊(duì)長(zhǎng)無(wú)限制．假定一個(gè)工作日為8小時(shí)，時(shí)間以分鐘為單位。[1]模擬一個(gè)工作日內(nèi)完成服務(wù)的個(gè)數(shù)及顧客平均等待時(shí)間t．[2]模擬100個(gè)工作日，求出平均每日完成服務(wù)的個(gè)數(shù)及每日顧客的平均等待時(shí)間。單服務(wù)員的排隊(duì)模型模擬w：總等待時(shí)間；ci：第i個(gè)顧客的到達(dá)時(shí)刻；bi：第i個(gè)顧客開始服務(wù)時(shí)刻；ei：第i個(gè)顧客服務(wù)結(jié)束時(shí)刻；xi：第i-1個(gè)顧客與第i個(gè)顧客之間到達(dá)的間隔時(shí)間；yi：對(duì)第i個(gè)顧客的服務(wù)時(shí)間。符號(hào)說(shuō)明c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xiei=bi+yibi=max(ci,ei-1)t思路分析初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時(shí)間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=xi,bi=xi

產(chǎn)生服務(wù)時(shí)間隨機(jī)數(shù)yi~[4,15]的均勻分布ei=bi+yi累計(jì)等待時(shí)間：w=w+bi-ci準(zhǔn)備下一次服務(wù)：i=i+1產(chǎn)生間隔時(shí)間隨機(jī)數(shù)xi~參數(shù)為0.1的指數(shù)分布ci=ci-1+xi

確定開始服務(wù)時(shí)間：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i輸出結(jié)果：完成服務(wù)個(gè)數(shù)：m=i

平均等待時(shí)間：t停止[1]模擬一日ToMatlab(simu1)[2]模擬100日ToMatlab(simu2)流程框圖用蒙特卡洛法解非線性規(guī)劃問(wèn)題基本假設(shè)試驗(yàn)點(diǎn)的第j個(gè)分量xj服從[aj,bj]內(nèi)的均勻分布．符號(hào)假設(shè)求解過(guò)程先產(chǎn)生一個(gè)隨機(jī)數(shù)作為初始試驗(yàn)點(diǎn),以后則將上一個(gè)試驗(yàn)點(diǎn)的第j個(gè)分量隨機(jī)產(chǎn)生，其它分量不變而產(chǎn)生一新的試驗(yàn)點(diǎn)．這樣，每產(chǎn)生一個(gè)新試驗(yàn)點(diǎn)只需一個(gè)新的隨機(jī)數(shù)分量．當(dāng)K>MAXK或P>MAXP時(shí)停止迭代．框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數(shù)xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?Y