隨機(jī)行為的模擬：隨機(jī)拋擲硬幣和骰子出現(xiàn)特定面的概率蒙特卡羅方法的計(jì)算機(jī)模擬1摘要對蒙特卡羅(MonteCarlo)方法的簡介并概述了蒙特卡羅方法的概念、應(yīng)用領(lǐng)域、求解步驟。以拋擲硬幣和骰子為例，論述了蒙特卡羅方法模擬隨機(jī)行為的基本思想和基本原理。給出了實(shí)現(xiàn)計(jì)算機(jī)模擬的MATLAB程序，并且通過最高達(dá)千萬次級別的計(jì)算機(jī)模擬試驗(yàn)，準(zhǔn)確地模擬了隨機(jī)拋擲硬幣和骰子出現(xiàn)特定面的概2關(guān)鍵詞蒙特卡羅(MonteCarlo)方法方法；計(jì)算機(jī)模擬；隨機(jī)行為；模擬；概率；MATLAB程序3弓I言3.1蒙特卡羅方法的概述：蒙特卡羅方法又稱統(tǒng)計(jì)模擬法、隨機(jī)抽樣技術(shù)，是一種隨機(jī)模擬方法，以概率和統(tǒng)計(jì)理論方法為基礎(chǔ)的一種計(jì)算方法是使用隨機(jī)數(shù)(或更常見的偽隨機(jī)數(shù))來解決很多計(jì)算問題的方法。將所求解的問題同一定的概率模型相聯(lián)系，用電子計(jì)算機(jī)實(shí)現(xiàn)統(tǒng)計(jì)模擬或抽樣，以獲得問題的近似解。為象征性地表明這一方法的概率統(tǒng)計(jì)特征，故借用賭城蒙特卡羅命名。3.2蒙特卡洛模擬法簡介：蒙特卡洛(MonteCarlo)模擬是一種通過設(shè)定隨機(jī)過程，反復(fù)生成時(shí)間序列，計(jì)算參數(shù)估計(jì)量和統(tǒng)計(jì)量，進(jìn)而研究其分布特征的方法。具體的，當(dāng)系統(tǒng)中各個(gè)單元的可靠性特征量已知，但系統(tǒng)的可靠性過于復(fù)雜，難以建立可靠性預(yù)計(jì)的精確數(shù)學(xué)模型或模型太復(fù)雜而不便應(yīng)用時(shí)，可用隨機(jī)模擬法近似計(jì)算出系統(tǒng)可靠性的預(yù)計(jì)值；隨著模擬次數(shù)的增多，其預(yù)計(jì)精度也逐漸增高。由于涉及到時(shí)間序列的反復(fù)生成，蒙特卡洛模擬法是以高容量和高速度的計(jì)算機(jī)為前提條件的，因此只是在近些年才得到廣泛推廣。3.3蒙特卡洛模擬法提出：蒙特卡羅方法于20世紀(jì)40年代美國在第二次世界大戰(zhàn)中研制原子彈的“曼哈頓計(jì)劃”計(jì)劃的成員S.M.烏拉姆和J.馮■諾伊曼首先提出。數(shù)學(xué)家馮■諾伊曼用馳名世界的賭城一摩納哥的MonteCarlo-來命名這種方法，為它蒙上了一層神秘色彩。在這之前，蒙特卡羅方法就已經(jīng)存在。1777年，法國Buffon提出用投針實(shí)驗(yàn)的方法求圓周率。3.4蒙特卡洛模擬法的應(yīng)用領(lǐng)域：、直接應(yīng)用蒙特卡洛模擬:應(yīng)用大規(guī)模的隨機(jī)數(shù)列來模擬復(fù)雜系統(tǒng)，得到某些參數(shù)或重要指標(biāo)。、蒙特卡洛積分:利用隨機(jī)數(shù)列計(jì)算積分，維數(shù)越高，積分效率越高。、MCMC:這是直接應(yīng)用蒙特卡洛模擬方法的推廣，該方法中隨機(jī)數(shù)的產(chǎn)生是采用的馬爾科夫鏈形式。、蒙特卡羅方法在金融工程學(xué)，宏觀經(jīng)濟(jì)學(xué)，生物醫(yī)學(xué)，計(jì)算物理學(xué)(如粒子輸運(yùn)計(jì)算、量子熱力學(xué)計(jì)算、空氣動(dòng)力學(xué)計(jì)算)等領(lǐng)域應(yīng)用廣泛。3.5蒙特卡羅解題歸結(jié)為三個(gè)主要步驟：、構(gòu)造或描述概率過程；、實(shí)現(xiàn)從已知概率分布抽樣;、建立各種估計(jì)量。4問題重述蒙特卡羅模擬的真正威力在于對隨機(jī)行為建模。從長期來看，一個(gè)事件的概率可以視為比值：概率P(a)有效的事件數(shù)事件的總數(shù)下面3個(gè)隨機(jī)模型：、拋擲一枚正規(guī)的硬幣、拋擲一個(gè)正規(guī)的骰子、拋擲一個(gè)不正規(guī)的骰子以剖析如何用蒙特卡羅方法模擬這些隨機(jī)行為，以及基于MATLAB軟件的計(jì)算機(jī)實(shí)現(xiàn)。5拋擲一枚正規(guī)的硬幣5.1過程分析拋擲一枚硬幣得到正面或反面的概率是1/2，我們可以把這種隨機(jī)事件和［0,1］內(nèi)的隨機(jī)數(shù)建立聯(lián)系。概率是長期平均值，于是拋很多次時(shí)出現(xiàn)次數(shù)的比例接近0.5。正面一［0,0.5］,反面f(0.5,1］。設(shè)x為［0,1］內(nèi)的隨機(jī)數(shù)，f(x)定義如下：f正面，0<x<0.5/(x)=<［反面，0.5<x<1f(x)將結(jié)果是正面或反面賦值到［0,1］內(nèi)的一個(gè)數(shù)，隨機(jī)賦值時(shí)我們可利用這個(gè)函數(shù)的累積性質(zhì)。拋很多次時(shí)能夠得到如下出現(xiàn)的百分比：隨機(jī)數(shù)區(qū)間出現(xiàn)的累積值出現(xiàn)的百分比x<000.000<x<00.5<x<1.010.505.2符號說明n:拋擲硬幣的次數(shù)COUNTER:記錄得到正面的次數(shù)的計(jì)數(shù)器COUNTER/n:隨機(jī)拋擲一枚正規(guī)的硬幣得到正面向上的概率，n的取值越大，得到的概率越接近理論值i:第i次拋擲硬幣5.3拋擲正規(guī)硬幣的蒙特卡羅算法輸入 模擬中生成的隨機(jī)拋擲硬幣的總次數(shù)n..輸出 拋擲硬幣時(shí)得到正面的概率.第1步 初始化：COUNTER=0.第2步 對于i=1,2,…,n，執(zhí)行3,4步.第3步 得到［0,1］內(nèi)的隨機(jī)數(shù)第4步若0<x.<0.5，則COUNTER=COUNTER+1.否則，COUNTER不變.第5步 計(jì)算P(正面)=COUNTER/n.第6步 輸出正面的概率P(正面)停止5.4MATLAB執(zhí)行代碼見附碼（注意：重復(fù)執(zhí)行代碼，若不更改初始值，則每次執(zhí)行，正面出現(xiàn)的次數(shù)會(huì)改變，正面出現(xiàn)的概率也隨次數(shù)改變而改變）5.5結(jié)果分析下表給出了對于不同的n由隨機(jī)數(shù)七得到的結(jié)果:拋硬幣次數(shù)正面出現(xiàn)的次數(shù)正面出現(xiàn)的百分比拋硬幣次數(shù)正面出現(xiàn)的次數(shù)正面出現(xiàn)的百分比100470.471000049690.496950024104820100000499640.499610005040.504010000005000450.5000隨著拋擲硬幣次數(shù)n的增大，正面出現(xiàn)的概率也逐漸的接近0.5,即次數(shù)的一半。當(dāng)實(shí)驗(yàn)次數(shù)達(dá)到百萬級以上時(shí)，模擬的概率值與理論值的誤差僅為土0.0001。6拋擲一個(gè)正規(guī)的骰子6.1過程分析一個(gè)骰子由點(diǎn)數(shù)｛1,2,3,4,5,6｝組成，拋擲一個(gè)正規(guī)的骰子，必須設(shè)計(jì)一種定義6個(gè)事件的方法。而每個(gè)數(shù)值出現(xiàn)的可能性相等，所以每個(gè)數(shù)值出現(xiàn)的概率是1/6。一個(gè)指定的數(shù)值出現(xiàn)的概率定義為：｛1，2,3,4,5,6｝中指定的數(shù)值出現(xiàn)的次 數(shù)試驗(yàn)的總數(shù)將區(qū)間［0,1］等分為6個(gè)子區(qū)間，隨機(jī)生成一個(gè)［0,1］上的隨機(jī)數(shù)，它等可能的屬于在6個(gè)子區(qū)間的某個(gè)區(qū)間上，這樣就可將隨機(jī)拋擲骰子出現(xiàn)的某個(gè)點(diǎn)數(shù)用隨機(jī)生成［0,1］上的隨機(jī)數(shù)屬于某個(gè)子區(qū)間來模擬。例如，隨機(jī)生成的一個(gè)隨機(jī)數(shù)屬于［0,1/6］，則認(rèn)為隨機(jī)拋擲一個(gè)正規(guī)的骰子出現(xiàn)的點(diǎn)數(shù)為1。6.2符號說明n:拋擲骰子的次數(shù)COUNTER:記錄次數(shù)的計(jì)數(shù)器COUNTER①/n:隨機(jī)拋擲一個(gè)正規(guī)骰子的概率，n的取值越大，得到的概率越接近理論值i:第i次拋擲骰子

j:骰子的點(diǎn)數(shù)6.3拋擲一個(gè)正規(guī)骰子的蒙特卡羅算法輸入輸出第1步第輸入輸出第1步第2步第3步擲出｛1,2,3,4,5,6｝的百分比或概率.將COUNTER1到COUNTER6初始化為0.對于i=1,2,…,n，執(zhí)行第3，4步得到一個(gè)隨機(jī)數(shù)，滿足0<x<1.0<xi<6得到一個(gè)隨機(jī)數(shù)，滿足0<x<1.0<xi<612—<x<—6i623—<x<—6i634—<x<—6i645-<x<—6i65<x<116 1COUNTER1=COUNTER1+1COUNTER2=COUNTER2+1COUNTER3=COUNTER3+1COUNTER4=COUNTER4+1COUNTER5=COUNTER5+1COUNTER6=COUNTER6+1第5步第6步第5步第6步停止輸出這些概率.6.4MATLAB執(zhí)行代碼見附碼(注意：重復(fù)執(zhí)行代碼，若不更改初始值，則每次執(zhí)行，概率也會(huì)改變)6.5結(jié)果分析下表給出了擲10、100、1000、10000和100000次正規(guī)骰子的結(jié)果，擲100000次時(shí)接近理論值1/6^0.1667骰子出現(xiàn)的次數(shù)10100100010000100000理論值10.10000.12000.17200.16330.16770.166720.20000.17000.16000.16210.16740.166730.20000.19000.16900.16150.16360.166740.10000.21000.17400.16910.16600.166750.30000.17000.15900.17580.16800.166760.10000.14000.16600.16820.16730.16677拋擲一個(gè)不正規(guī)的骰子7.1過程分析在拋擲一個(gè)正規(guī)的骰子的基礎(chǔ)上考慮每個(gè)事件不是等可能出現(xiàn)的模型：拋擲一個(gè)不正規(guī)的骰子。假定給骰子的某幾個(gè)面加重使結(jié)果發(fā)生偏移。將區(qū)間［0,1］分成6個(gè)不等長的子區(qū)間，子區(qū)間的長度代表骰子出現(xiàn)的某個(gè)值的概率。例如，隨機(jī)生成的一個(gè)隨機(jī)數(shù)屬于［0.4,0.7］，對應(yīng)著隨機(jī)拋擲一個(gè)不正規(guī)骰子出現(xiàn)的點(diǎn)數(shù)為4的概率是0.3。詳細(xì)見下表：拋擲一個(gè)不正規(guī)骰子的概率分布和函數(shù)關(guān)系不正規(guī)骰子的概率分布 算法中采用的函數(shù)關(guān)系骰子出現(xiàn)的值出現(xiàn)的概率七的值賦值10.1[0.0,0.1]120.1(0.1,0.2]230.2(0.2,0.4]340.3(0.4,0.7]450.2(0.7,0.9]560.1(0.9,1.0]67.2符號說明n:拋擲骰子的次數(shù)COUNTER:記錄次數(shù)的計(jì)數(shù)器COUNTER①/n:隨機(jī)拋擲一個(gè)不正規(guī)骰子的概率i:第i次拋擲骰子j：骰子的點(diǎn)數(shù)7.3拋擲一個(gè)不正規(guī)骰子的蒙特卡羅算法輸入 模擬中隨機(jī)拋擲骰子的總數(shù)n.輸出 擲出｛1,2,3,4,5,6｝的百分比或概率.第1步 將COUNTER1到COUNTER6初始化為0.第2步對于i=1,2,…,n，執(zhí)行第3,4步第3步 得到一個(gè)隨機(jī)數(shù)，滿足0<罕1.第4步若x屬于如下的區(qū)間，則相應(yīng)的COUNTER+1.0<x<0>1COUNTER1=COUNTER1+10>1<x.<0.2COUNTER2=COUNTER2+10.2<x<0.4 COUNTER3=COUNTER3+10.4<x<0.7COUNTER4=COUNTER4+10.7<x<0.9COUNTER5=COUNTER5+10.9<x<1.0COUNTER6=COUNTER6+1第5步 計(jì)算擲出j={1,2,3,4,5,6}的概率為COUNTER(j)./n.第6步 輸出這些概率.停止7.4MATLAB執(zhí)行代碼見附碼(注意：重復(fù)執(zhí)行代碼，若不更改初始值，則每次執(zhí)行，概率也會(huì)改變)7.5結(jié)果分析卜表給出了擲100、1000、5000、10000、40000和1000000次不正規(guī)骰子的結(jié)果骰子出現(xiàn)次數(shù)1001000500010000400001000000理論值10.06000.10600.10300.09810.09960.09950.120.10000.10500.09540.10350.09840.10020.130.19000.19200.19460.20600.19890.19980.240.35000.30600.29740.29700.29840.30020.350.23000.19400.20580.19590.20130.20050.260.07000.09700.10380.09950.10340.09990.1為使試驗(yàn)的結(jié)果接近理論值需要更多的試驗(yàn)次數(shù)。當(dāng)n=1000000，即重復(fù)進(jìn)行100萬次蒙特卡羅模擬試驗(yàn)后，得到的點(diǎn)數(shù)1，2,3,4，5，6的概率依次為P=0.0995 0.1002 0.1998 0.3002 0.2005 0.09998小結(jié)蒙特卡羅方法與一般計(jì)算方法有很大區(qū)別，一般計(jì)算方法對于解決多維或因素復(fù)雜的問題非常困難，而蒙特卡羅方法對于解決這方面的問題卻比較簡單。其特點(diǎn)如下：椅直接追蹤粒子，物理思路清晰，易于理解。椅采用隨機(jī)抽樣的方法，較真切的模擬粒子輸運(yùn)的過程，反映了統(tǒng)計(jì)漲落的規(guī)律。椅不受系統(tǒng)多維、多因素等復(fù)雜性的限制，是解決復(fù)雜系統(tǒng)粒子輸運(yùn)問題的好方法。MC程序結(jié)構(gòu)清晰簡單。椅研究人員采用MC方法編寫程序來解決粒子輸運(yùn)問題，比較容易得到自己想得到的任意中間結(jié)果，應(yīng)用靈活性強(qiáng)。椅MC方法主要弱點(diǎn)是收斂速度較慢和誤差的概率性質(zhì)，其概率誤差正比于，如果單純以增大抽樣粒子個(gè)數(shù)N來減小誤差，就要增加很大的計(jì)算量。蒙特卡羅算法已成為計(jì)算數(shù)學(xué)中不可缺少的組成部分，這主要是因?yàn)橐韵略颍阂蝹鹘y(tǒng)的分析方法受到了問題復(fù)雜性的限制。椅MC方法直觀，對實(shí)驗(yàn)者很有吸引力。椅計(jì)算機(jī)變得更快更便宜。椅量子理論的發(fā)展為我們提供了輻射與物質(zhì)相互作用的截面數(shù)據(jù)。9參考文獻(xiàn)葉其孝、姜啟源,FrankR..Gioradano（美）等著，數(shù)學(xué)建模（原書第3版）[M],北京機(jī)械工業(yè)出版社,2005MBA智庫百科?蒙特卡羅方法/wiki/%E8%92%99%E7%89%B9%E5%8D%A1%E7%BD%97%E6%A8%A1%E6%8B%9F百度文庫?蒙特卡洛模擬法/view/4d5d6307cc175527072208a1.html[4]百度百科?蒙特卡羅模擬/view/284709.htm[5]新浪博客-Dwass（1957）蒙特卡羅模擬/s/blog_4b700c4c0100i28q.html10附碼10.1拋擲一枚正規(guī)的硬幣COUNTER=0;

n=100; %硬幣拋擲次數(shù)，賦初始值，大小可調(diào)節(jié)x=rand(1,n); %產(chǎn)生了n個(gè)隨機(jī)數(shù)fori=1:nifx(i)>=0&&x(i)<=0.5COUNTER=COUNTER+1;endendCOUNTER %正面出現(xiàn)的次數(shù)P=COUNTER./n %正面出現(xiàn)的概率10.2拋擲一個(gè)正規(guī)的骰子forj=1:6COUNTER(j)=0;end%骰子拋擲次數(shù)，賦初始值，大小可調(diào)節(jié)%骰子拋擲次數(shù)，賦初始值，大小可調(diào)節(jié)%產(chǎn)生了n個(gè)隨機(jī)數(shù)x=rand(1,n);fori=1:nifx(i)>=0&&x(i)<=1/6COUNTER(1)=COUNTER(1)+1;COUNTER(2)=COUNTER(2)+1;COUNTER(3)=COUNTER(3)+1;COUNTER(2)=COUNTER(2)+1;COUNTER(3)=COUNTER(3)+1;COUNTER(4)=COUNTER(4)+1;COUNTER