隨機(jī)過(guò)程實(shí)驗(yàn)

PAGEPAGE3隨機(jī)過(guò)程實(shí)驗(yàn)講義劉繼成華中科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院2011-2012年上半年為華中科技大學(xué)數(shù)學(xué)系本科生講授隨機(jī)過(guò)程課程參考資料目錄TOC\o"1-1"\h\z\u前言 1第一章Matlab簡(jiǎn)介 2第二章簡(jiǎn)單分布的模擬 6第三章基本隨機(jī)過(guò)程 9第四章Markov過(guò)程 12第五章模擬的應(yīng)用和例子 16附錄各章的原程序 51參考文獻(xiàn) 81前言若想檢驗(yàn)數(shù)學(xué)模型是否反映客觀現(xiàn)實(shí)，最自然的方法是比較由模型計(jì)算的理論概率和由客觀試驗(yàn)得到的經(jīng)驗(yàn)頻率。不幸的是，這兩件事都往往是費(fèi)時(shí)的、昂貴的、困難的，甚至是不可能的。此時(shí)，計(jì)算機(jī)模擬在這兩方面都可以派上用場(chǎng)：提供理論概率的數(shù)值估計(jì)與接近現(xiàn)實(shí)試驗(yàn)的模擬。模擬的第一步自然是在計(jì)算機(jī)程序的算法中如何產(chǎn)生隨機(jī)性。程序語(yǔ)言，甚至計(jì)算器，都提供了“隨機(jī)”生成[0，1]區(qū)間內(nèi)連續(xù)數(shù)的方法。因?yàn)槊看芜\(yùn)行程序常常生成相同的“隨機(jī)數(shù)”，因此這些數(shù)被稱(chēng)為偽隨機(jī)數(shù)。盡管如此，對(duì)于多數(shù)的具體問(wèn)題這樣的隨機(jī)數(shù)已經(jīng)夠用。我們將假定計(jì)算機(jī)已經(jīng)能夠生成[0，1]上的均勻隨機(jī)數(shù)。也假定這些數(shù)是獨(dú)立同分布的，盡管它們常常是周期的、相關(guān)的、……。……本講義的安排如下，第一章是Matlab簡(jiǎn)介，從實(shí)踐動(dòng)手角度了解并熟悉Matlab環(huán)境、命令、幫助等，這將方便于Matlab的初學(xué)者。第二章是簡(jiǎn)單隨機(jī)變量的模擬，只給出了常用的Matlab模擬語(yǔ)句，沒(méi)有堆砌同一種變量的多種模擬方法。對(duì)于沒(méi)有列舉的隨機(jī)變量的模擬，以及有特殊需求的讀者應(yīng)該由這些方法得到啟發(fā)，或者參考更詳細(xì)的其他文獻(xiàn)資料。第三章是基本隨機(jī)過(guò)程的模擬。主要是簡(jiǎn)單獨(dú)立增量過(guò)程的模擬，多維的推廣是直接的。第四章是Markov過(guò)程的模擬。包括服務(wù)系統(tǒng)，生滅過(guò)程、簡(jiǎn)單分支過(guò)程等。第五章是這些模擬的應(yīng)用。例如，計(jì)算概率、估計(jì)積分、模擬現(xiàn)實(shí)、誤差估計(jì)，以及減小方差技術(shù)，特別給讀者提供了一些經(jīng)典問(wèn)題的模擬，通過(guò)這些問(wèn)題的模擬將會(huì)更加牢固地掌握實(shí)際模擬的步驟。平穩(wěn)過(guò)程的模擬、以及利用平穩(wěn)過(guò)程來(lái)預(yù)測(cè)的內(nèi)容并沒(méi)有包含在本講義之內(nèi)，但這絲毫不影響該內(nèi)容的重要性，這也是將會(huì)增補(bǔ)進(jìn)來(lái)的主要內(nèi)容之一。希望讀者碰到類(lèi)似的問(wèn)題時(shí)能夠查閱相關(guān)資料解決之。各章的內(nèi)容包括了模擬的基本思路和Matlab代碼。源程序包展示了對(duì)各種隨機(jī)過(guò)程與隨機(jī)機(jī)制的有效模擬和可視化的基本技術(shù)，試圖強(qiáng)調(diào)matlab自然處理矩陣和向量的方法，目標(biāo)是為涉及應(yīng)用隨機(jī)模擬的讀者在準(zhǔn)備自己的程序代碼時(shí)找到靈感和想法。建議讀者在了解了模擬的基本方法之后就著手解決自己感興趣的實(shí)際問(wèn)題。對(duì)實(shí)際具體問(wèn)題的解決有助于更深刻理解模擬的思想、也會(huì)在具體應(yīng)用中拓展現(xiàn)有的模擬方法。

第一章Matlab簡(jiǎn)介若你想在計(jì)算機(jī)上運(yùn)行Matlab，點(diǎn)擊：開(kāi)始/程序/MATLAB6.5，這樣將會(huì)出現(xiàn)有三個(gè)窗口的交互界面。如果你是初學(xué)者，可以先瀏覽一下Matlab的指導(dǎo)材料，點(diǎn)擊：Help/MATLABHelp,打開(kāi)窗口左邊的“MATLAB”一節(jié)即可看到相關(guān)的內(nèi)容。就自己而言，我學(xué)習(xí)Matlab更喜歡的方式是：輸入并運(yùn)行一些命令、觀察出現(xiàn)的結(jié)果，然后查閱想了解的幫助文件。這也正是本節(jié)的方法。在“commandwindow”窗口中顯示有提示符“>>”,在提示符后輸入下面的命令，按回車(chē)鍵即可運(yùn)行并顯示相應(yīng)的結(jié)果。當(dāng)然，不要輸入行號(hào)、也不必輸入后面的注釋。在這個(gè)部分討論的Matlab文件有:rando.m，vrando.m，show.m。一、Matlab初步1：2*9Matlab當(dāng)作計(jì)算器用。2：sin(1)Matlab僅顯示四位小數(shù)，但保存的更多！3：formatlong顯示更多位小數(shù)。4：sin(1)5：2^9996：formatshort7：x=sin(1)將計(jì)算結(jié)果存在變量中。8：x顯示x的值。9：x=rand(10,1)x是包含有10個(gè)[0，1]上均勻分布隨機(jī)數(shù)的集合，它是一個(gè)列向量或者是10×1的矩陣。10：x+5x的每個(gè)分量都加5。11：1000*xx的每個(gè)分量都乘以1000。12：x=rand(10,7)10×7的隨機(jī)數(shù)矩陣。若想重復(fù)此命令或其他命令，按住向上的光標(biāo)鍵直至看到想重復(fù)的命令。13：x=rand(1000,1)將1000換成更大的數(shù)試試。14：x=rand(1000,1);用分號(hào)“；”可以不顯示結(jié)果。15：help顯示標(biāo)準(zhǔn)的幫助列表。16：helpelmat顯示關(guān)于初等矩陣的幫助，包括命令“rand”。17：helprand直接顯示“rand”的幫助。18：x(1:20,1)取出x的第一列中的1-20個(gè)數(shù)。19：helppunctMatlab中關(guān)于標(biāo)點(diǎn)符號(hào)的用法。20：max(x)21：mean(x)22：sum(x)23：median(x)x的中位數(shù)。24：cumsum(x)x的分量累計(jì)和向量。25：y=sort(rand(10,1))由小到大排序后的向量。26：hist(x)作出x的直方圖。27：hist(x,30)用30個(gè)方柱代替缺省的10個(gè)。28：y=-log(x)對(duì)x的分量取自然對(duì)數(shù)。29：hist(y,30)多數(shù)的y的分量只是接近0的，但有些是和6差不多大的，y中的數(shù)被稱(chēng)為指數(shù)分布隨機(jī)數(shù)。30：z=randn(1000,1);生成1000個(gè)標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。31：hist(z,30)直方圖是鐘形的。對(duì)大于1000的數(shù)試試結(jié)果。二、獲取更多幫助32：如果你想查找不會(huì)使用的命令，可以點(diǎn)擊：：Help/MATLABHelp，打開(kāi)左邊的“MATLAB”節(jié)，選擇“Functions–CategoricalList”即可。據(jù)我所知，這是尋求幫助的最好方法。三、畫(huà)出數(shù)據(jù)點(diǎn)33：plot(x(1:10),’*’);用“*”描出x的前10個(gè)點(diǎn)。注意兩個(gè)單引號(hào)為英文的單引號(hào)，下同。34：plot(x-0.5);向下平移0.5，描出述據(jù)點(diǎn)，且將其連成線。35：holdon將下面的圖形加到上面的圖形中。36：plot(cumsum(x-0.5),’r’);將這個(gè)結(jié)果圖畫(huà)到上面的圖形中?！啊痳’”表示用紅色的線繪出，而缺省的顏色為藍(lán)色。37：zoomon用鼠標(biāo)點(diǎn)擊可放大圖形，雙擊回到原始的尺寸。38：clf清除當(dāng)前的圖形。39：z=randn(1000,1);生成1000個(gè)標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。40：w=z+randn(1000,1);生成依賴(lài)z的隨機(jī)數(shù)。41：plot(z,w,’*’);作出（z,w）的圖形。42：axis([-33-44]);顯示x在[-3，3]與y在[-4，4]范圍的圖形。四、作圖函數(shù)43：clf44：ezplot(’sin(x)’,[03*pi]);畫(huà)出正弦函數(shù)的圖形。45：holdon46：t=0:0.01:3*pi;定義一個(gè)時(shí)間點(diǎn)向量，間隔為0.01。47：tt為一行向量。48：t=t’現(xiàn)在t為一列向量。49：plot(t,sin(5*t),’r’);用紅色畫(huà)sin(t)關(guān)于t的函數(shù)。顯然，函數(shù)ezplot不能設(shè)置圖形的顏色。50：title(’sin(x)andsin(5x)’)給圖形加上更恰當(dāng)?shù)臉?biāo)題。五、運(yùn)行現(xiàn)有的Matlab程序51：上網(wǎng)下載或者拷貝一些編輯好的Matlab程序到自己的電腦中。52：如果在你電腦的某個(gè)文件夾中有現(xiàn)成的Matlab程序（*.m）,可以設(shè)置“CurrentDirectory”(CommandWindow窗口的上面)為該文件夾即可運(yùn)行這些程序。53：如果在你電腦里的幾個(gè)文件夾里都有Matlab程序，點(diǎn)擊菜單中：File/SetPath/AddFolder,加入所有這些文件夾，最后選擇“Save”。當(dāng)你在CommandWindow窗口鍵入一命令后，Matlab會(huì)在所有的這些文件夾中查找這個(gè)命令名。六、拋硬幣54：3<5不等式滿足結(jié)果為1。55：5<3結(jié)果是0。56：x=rand(20,1)前面已輸入過(guò)類(lèi)似的命令。輸入“x=”，然后用向上的光標(biāo)鍵往回翻看，找到后將1000改為20。57：x>0.5對(duì)x的所有分量檢查該不等式。58：z=1+(x>0.5)z的值為1或者2。這有點(diǎn)像拋硬幣，1為正面，2為反面。59：show(z,’正反’)這是一個(gè)名字為show的程序，有兩個(gè)變量，一個(gè)是自然數(shù)向量，一個(gè)是用來(lái)與每個(gè)數(shù)相對(duì)應(yīng)顯示的字符串。它是自己編制的程序，保存在：d:\MATLAB6p5\work\show.m。60：show(1+(rand(1500,1)>0.5),’正反’)生成1500個(gè)拋硬幣的結(jié)果?，F(xiàn)在按下向上的光標(biāo)鍵/回車(chē)，就會(huì)得到很多拋硬幣的結(jié)果。你找到連續(xù)出現(xiàn)正面最多的個(gè)數(shù)了嗎？61：show(1+(rand(1500,1)>0.5),’O-’)可以通過(guò)改變顯示的字符來(lái)簡(jiǎn)化剛才的問(wèn)題。用向上的光標(biāo)鍵很容易更改前面的命令來(lái)實(shí)現(xiàn)它。這些語(yǔ)句對(duì)拋硬幣的問(wèn)題當(dāng)然是足夠了，因?yàn)樗挥袃蓚€(gè)結(jié)果。但對(duì)其他，像擲色子，的隨機(jī)試驗(yàn)，“rando.m”將更加有用，這也是自己編制的程序，保存在：d:\MATLAB6p5\work\show.m。62：d=[111111]/6擲色子的結(jié)果概率是一個(gè)行向量（或者1×6矩陣）。63：sum(d)確認(rèn)它們的和為1！64：rando(d)用這些概率去模擬擲色子的每個(gè)結(jié)果。用向上的光標(biāo)鍵重復(fù)這個(gè)命令幾次。模擬擲色子的另一個(gè)簡(jiǎn)單的方法是放大均勻分布隨機(jī)數(shù)后取整，floor(1+6*rand(1))。65：vrando([1111]/4,20)程序rando的向量版本。每個(gè)數(shù)是等概率出現(xiàn)的。66：show(vrando([1111]/4,100),’BGSU’)隨機(jī)地生成字符B、G、S和U。出現(xiàn)BUGS之前，BGSU出現(xiàn)了嗎？七、寫(xiě)一個(gè)Matlab程序你將創(chuàng)建一個(gè)新的Matlab程序，名字為mywalk.m，用它來(lái)模擬100步的隨機(jī)游動(dòng)。在“file”菜單下有一個(gè)空白的按鈕，按下它即打開(kāi)一個(gè)新的編輯窗口。在那個(gè)窗口里，分行輸入下面的命令，然后保存該程序?yàn)閙ywalk.m。如果你保存在新的文件夾里，請(qǐng)確認(rèn)這個(gè)文件夾是否已加入到Path中或者改變?yōu)镃urrentDirectory。67：n=100;選取步數(shù)。68：x=rand(n,1);生成均勻分布隨機(jī)數(shù)。69：y=2*(x>0.5)-1;轉(zhuǎn)換這些數(shù)到為-1和+1。70：z=cumsum(y);計(jì)算y的累積和。71：clf72：plot(z)畫(huà)出z的第1,2,3,...等的值。在commandwindow窗口中輸mywalk，運(yùn)行（按回車(chē)）該程序，然后用光標(biāo)鍵多次重復(fù)它。如果有錯(cuò)誤提示，檢查你的輸入是否是正確的。73：運(yùn)行幾次后，你或許想一次就生成一個(gè)更長(zhǎng)的字符串。到此目的的一個(gè)好的方法是將mywalk.m改為帶參數(shù)的Matlabfunction，這樣就可以調(diào)用它。74：在你的程序中，將行“n=100;”替換為function[z]=mywalk(n)這樣，mywalk是一個(gè)帶參數(shù)n的函數(shù)(生成序列的長(zhǎng)度)，返回變量z。函數(shù)里面的變量（比如y）是內(nèi)部變量，它的值不被帶到函數(shù)外面。就像sin和rand一樣，函數(shù)mywalk返回一個(gè)值（向量z）?；氐絚ommandwindow窗口輸入：75：mywalk(1000);運(yùn)行參數(shù)為1000的程序mywalk。八、矩陣76：M=rand(6,6)6×6的隨機(jī)數(shù)矩陣。77：M(2,:)取出矩陣M的第2行。78：M(:,4)取出矩陣M的第4列。79：diag(M)取出矩陣M的對(duì)角線元素。80：sum(M)矩陣列求和。81：sum(M’)’對(duì)矩陣M的行求和。“’”表示轉(zhuǎn)置。九、Markov鏈在第66行中，序列中字母的出現(xiàn)是相互獨(dú)立的。我們將建立下面的一種情形，B通常跟隨在U之后，但決不跟在G之后。出現(xiàn)B后，依概率向量[0.20.60.20]選擇下一個(gè)字母。G出現(xiàn)后，又以另一不同的概率向量出現(xiàn)下一個(gè)字母，以此類(lèi)推。為此，我們將創(chuàng)建名字為BGSU_markov.m的新程序。打開(kāi)一個(gè)新的編輯窗口，輸入下面的命令，然后再命令窗口輸入BGSU_markov運(yùn)行之。82：P=[[0.20.60.20];[00.20.60.2];[0.200.20.6];[0.60.200.2]];P是一個(gè)4×4矩陣。每一行表明將以多大的概率選擇下一個(gè)字母。第一行即是數(shù)字1之后（對(duì)應(yīng)字母B）的概率，第二行是數(shù)字2之后（G）的概率等等。83：x(1)=rando([1111]/4);隨機(jī)地選擇第一個(gè)狀態(tài)。84：fori=1:399,85：x(i+1)=rando(P(x(i),:));這是非常明智的：無(wú)論在哪個(gè)時(shí)刻，x(i)的值是多少，P(x(i),:)總是矩陣P的第x(i)行。該行的概率作為rando的參數(shù)來(lái)生成下一個(gè)狀態(tài)。86：end87：show(x,’BGSU’);88：hist(x,4);

第二章簡(jiǎn)單分布的模擬Matlab里生成[0，1]上的均勻隨機(jī)數(shù)的語(yǔ)句是：rand(1,1);rand(n,m)。一旦有了[0，1]上均勻隨機(jī)數(shù)，則我們就能夠做下面的事情。在這個(gè)部分討論的Matlab文件有:simexp.m,simpareto.m，simparetonrm.m,simdiscr.m,simbinom.m,simgeom.m,distrmu.m,distrstat.m。一、一般連續(xù)分布（逆變換法、拒絕法、Hazard率方法）生成有連續(xù)分布函數(shù)隨機(jī)數(shù)的一般方法是用反函數(shù)法。設(shè)G(y)=F^{-1}(y)，如果u(1)...,u(n)是服從（0，1））上均勻分布的隨機(jī)數(shù)，那么G(u(1)),...,G(u(n))就是分布函數(shù)為F(x)的隨機(jī)數(shù)。比如，指數(shù)分布，Pareto分布等。1、指數(shù)分布simexp.m事件以強(qiáng)度lambda的時(shí)間隨機(jī)地發(fā)生，即事件在[t，t＋h]時(shí)間內(nèi)發(fā)生的可能性是lambda×h，令t為事件發(fā)生前的等待時(shí)間。t=-log(rand)/lambda;%服從參數(shù)為lambda的指數(shù)分布Exp(lambda)的隨機(jī)數(shù)。t=-log(rand(1,m))./lambda;%服從Exp(lambda)的m維行向量。2、Pareto分布simpareto.m概率密度函數(shù):f(x)=alpha/(1+x)^(1+alpha),x>0累積分布函數(shù):F(x)=1-(1+x)^(-alpha)。這是帶有所謂重尾分布中最簡(jiǎn)單的分布列子。產(chǎn)生一個(gè)均勻分布的樣本，并用分布函數(shù)的反函數(shù)：sample=(1-rand(1,npoints)).^(-1/alpha)-1;3、標(biāo)準(zhǔn)Pareto分布simparetonrm.m概率密度函數(shù):f(x)=gamma*alpha/(1+gamma*x)^(1+alpha),x>0累積分布函數(shù):F(x)=1-(1+gamma*x)^(-alpha)其中，參數(shù)gamma是用來(lái)控制期望值的。在分布有重尾的情況下，若1<alpha<2，期望值存在且等于1/(gamma*(alpha-1))。產(chǎn)生一個(gè)均勻分布樣本并用反函數(shù)法：sample=((1-rand(M,N)).^(-1/alpha)-1)./gamma;二、一般離散分布simdiscr.m（除了上面的外，還有Alias方法）假設(shè)給出n個(gè)概率p=[p(1)...p(n)],滿足sum(p)=1且分量p(j)是非負(fù)的。為產(chǎn)生m個(gè)服從這個(gè)分布的隨機(jī)數(shù)，可以想象將區(qū)間（0，1）以p(1)...,p(n)為長(zhǎng)度間隔做一個(gè)劃分.產(chǎn)生一個(gè)均勻隨機(jī)數(shù)，如果該數(shù)落在第j個(gè)間隔中，賦予此離散分布值j，重復(fù)m次。uni=rand(1,m);cumprob=[0cumsum(p)];sample=zeros(1,m);forj=1:nind=find((uni>cumprob(j))&(uni<=cumprob(j+1)));sample(ind)=j;end1、0-1分布(rand(1,m)<=p);%生成m個(gè)以概率p為1，概率1-p為0的隨機(jī)數(shù)（m維行向量）。三、特殊分布1、二項(xiàng)分布simbinom.m將每次成功的概率為p的試驗(yàn)獨(dú)立做n次，設(shè)x是成功的個(gè)數(shù)x=sum(rand(n,m)<=p);%x是服從Bin(n,p)的m維隨機(jī)數(shù)向量。2、幾何分布simgeom.m實(shí)驗(yàn)每次成功的概率為p，設(shè)x為第一次成功前失敗的次數(shù)。x=floor(-log(rand(1,m))./(-log(1-p)));%服從參數(shù)為p的幾何分布Ge(p)的m維隨機(jī)數(shù)行向量。floor是取小于它的最小整數(shù)的函數(shù)。3、泊松分布Matlab的統(tǒng)計(jì)工具箱含有產(chǎn)生泊松分布隨機(jī)數(shù)的命令，為poissrnd。poissrnd(lambda);poissrnd(lambda,n,m);%產(chǎn)生參數(shù)為lambda的泊松分布Po(lambda)隨機(jī)數(shù)的n×m矩陣。如果沒(méi)有上面的命令，也可以用如下的命令替代之。arrival=cumsum(-log(rand(1,5)./lambda));n=length(find(arrival<=lambda));%find是找出非0值所在的位置，length是它的維數(shù)。4、正態(tài)分布高斯分布，或正態(tài)分布的隨機(jī)數(shù)用Matlab生成的命令是randn。randn(1,m);%服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)的m維隨機(jī)數(shù)行向量。randn(n,m);%每個(gè)分量是服從N(0,1)的n×m矩陣。mu+sigma.*randn(1,m);%m個(gè)服從N(mu,sigma^2)分布的隨機(jī)數(shù)四、離散試驗(yàn)的模擬1、從{1,…,n}中任取一個(gè)。int(n*rand(1,m)+1);從{1,…,n}中任取不可重復(fù)兩個(gè)。a=int(n*rand(1)+1);b=int(n*rand(1)+1);while(a=b)b=int(n*rand(1)+1);end2、隨機(jī)子集模擬集合{1,…,n}的隨機(jī)子集，我們是定義序列S(j)={0,1}，S(j)=1即表示將j在S中。每個(gè)S(j)，j=1,…,n,以1/2的概率獨(dú)立選擇0或1。forj=1:ns(j)=int(rand(1)+1);j=j+1;end3、隨機(jī)排列假如我們向隨機(jī)地排列a(1),…,a(n)，一個(gè)快速的方法是每一次互換兩個(gè)數(shù)的位置，共n-1次。forj=n:2N=int(j*rand(1)+1);y=a(N);a(N)=a(j);a(j)=y;end五、外部參數(shù)的隨機(jī)數(shù)產(chǎn)生器distrmu.m1、在一些模擬程序中(比如更新過(guò)程)，把概率分布作為一個(gè)外部參數(shù)來(lái)傳遞是很方便的。這是通過(guò)創(chuàng)建一個(gè)MATLAB函數(shù)來(lái)實(shí)現(xiàn)。例如，@rand,@simpareto。分布的參數(shù)是作為數(shù)組來(lái)傳遞的(在rand中為空數(shù)組{}，simpareto中為參量alpha{1.4}。distrmu.m是一個(gè)表-查找函數(shù)，從它的參數(shù)列表中提取期望參數(shù)的外部隨機(jī)數(shù)發(fā)生器，如：mu=distrmu(@simparetonrm,{1.4,2.5});2、平穩(wěn)分布distrstat.m假設(shè)有一個(gè)分布函數(shù)為F(x)、期望值為mu的分布，則它的平穩(wěn)或均衡分布的分布函數(shù)是G(x)=1/mu*int_0^x(1-F(y))dy。例如，密度函數(shù)為f(x)=2-2*x,0<=x<=1的線性分布是（0，1）上均勻分布上的平穩(wěn)分布。參數(shù)為(alpha-1)Pareto分布是參數(shù)為alpha的Pareto分布，參數(shù)為lambda的指數(shù)分布的平穩(wěn)分布就是自己。這將出現(xiàn)在平穩(wěn)更新過(guò)程或者排隊(duì)系統(tǒng)的平穩(wěn)版本的例子中。distrstat.m是一個(gè)表-查找函數(shù)，給定一個(gè)外部隨機(jī)數(shù)生成器，返回它的平穩(wěn)分布隨機(jī)數(shù)生成器。兩個(gè)參數(shù)都以數(shù)組的形式給出。至于應(yīng)用，可參見(jiàn)平穩(wěn)更新計(jì)數(shù)過(guò)程。例如：[statdist,statpar]=distrstat(@rand,{});

第三章基本隨機(jī)過(guò)程兩個(gè)基本機(jī)制是離散時(shí)間的隨機(jī)游動(dòng)與連續(xù)時(shí)間的泊松過(guò)程。這些過(guò)程是基于獨(dú)立的簡(jiǎn)單模擬算法的原型。擴(kuò)展到二維和三維中的模擬是直接的。在這個(gè)部分討論的Matlab文件有：ranwalk.m,brownian.m,poissonti.m,poissonjp.m,ranwalk2d.m,ranwalk3d.m,bm3plot.m,poisson2d.m,poisson3d.m一、一維情形1、隨機(jī)游動(dòng)1).簡(jiǎn)單隨機(jī)游動(dòng)ranwalk.m“從0開(kāi)始，向前跳一步的概率為p，向后跳一步的概率為1－p”p=0.5;y=[0cumsum(2.*(rand(1,n-1)<=p)-1)];%n步。plot([0:n-1],y);%畫(huà)出折線圖。2).隨機(jī)步長(zhǎng)的隨機(jī)游動(dòng)選取任一零均值的分布為步長(zhǎng),比如，均勻分布。x=rand(1,n)-1/2;y=[0cumsum(x)-1)];plot([0:n],y);2、布朗運(yùn)動(dòng)brownian.m這是連續(xù)情形的對(duì)稱(chēng)隨機(jī)游動(dòng),每個(gè)增量W(s+t)-W(s)是高斯分布N(0,t)，不相交區(qū)間上的增量是獨(dú)立的。典型的模擬它方法是用離散時(shí)間的隨機(jī)游動(dòng)來(lái)逼近。n=1000;dt=1;y=[0cumsum(dt^0.5.*randn(1,n))];%標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。plot(0:n,y);3、泊松過(guò)程產(chǎn)生隨機(jī)事件，滿足:(i)事件彼此獨(dú)立發(fā)生,(ii)兩次或更多事件不會(huì)同時(shí)發(fā)生,(iii)事件以常數(shù)強(qiáng)度發(fā)生。[0，t]內(nèi)事件發(fā)生的次數(shù)是期望值為lambda*t的泊松分布。計(jì)數(shù)過(guò)程N(yùn)(t)是泊松過(guò)程。連續(xù)兩次發(fā)生的時(shí)間間隔服從參數(shù)為lambda的指數(shù)分布。1).固定步數(shù)poissonjp.m%模擬n個(gè)服從Exp(lambda)的間隔時(shí)間interarr=[0-log(rand(1,n))./lambda];stairs(cumsum(interarr),0:n);%stairs畫(huà)出的是水平線條。2).固定時(shí)間區(qū)間，一個(gè)過(guò)程固定時(shí)間區(qū)間[0tmax]。在該區(qū)間內(nèi)事件發(fā)生的總數(shù)是期望值為lambda*tmax的泊松分布。在給定事件發(fā)生次數(shù)的條件下,事件服從該區(qū)間上的均勻分布。%總點(diǎn)數(shù)是服從泊松分布的。npoints=poissrnd(lambda*tmax);%在點(diǎn)數(shù)為N的條件下，點(diǎn)是均勻分布的。if(npoints>0)arrt=[0;sort(rand(npoints,1)*tmax);elsearrt=0;end%畫(huà)出計(jì)數(shù)過(guò)程stairs(arrt,0:npoints);3).固定時(shí)間區(qū)間，N個(gè)過(guò)程poissonti.m它被復(fù)雜化為前面算法的向量形式。到達(dá)時(shí)間間隔為指數(shù)分布的更新過(guò)程也將使用相同的算法。%將0賦給到達(dá)時(shí)間。tarr=zeros(1,nproc);%將指數(shù)分布的時(shí)間間隔求和作為矩陣的列。i=1;while(min(tarr(i,:))<=tmax)tarr=[tarr;tarr(i,:)-log(rand(1,nproc))/lambda];i=i+1;end%畫(huà)出計(jì)數(shù)過(guò)程stairs(tarr,0:size(tarr,1)-1);二、高維情形1、二維隨機(jī)游動(dòng)ranwalk2d.m在(u,v)坐標(biāo)平面上畫(huà)出點(diǎn)(u(k),v(k)),k=1:n,其中(u(k))和(v(k))是一維隨機(jī)游動(dòng)。例子程序是用四種不同顏色畫(huà)了同一隨機(jī)游動(dòng)的四條軌道。n=100000;colorstr=['b''r''g''y'];fork=1:4z=2.*(rand(2,n)<0.5)-1;x=[zeros(1,2);cumsum(z')];col=colorstr(k);plot(x(:,1),x(:,2),col);holdonendgrid2、三維隨機(jī)游動(dòng)ranwalk3d.m三維空間和上面的一樣。p=0.5;n=10000;colorstr=['b''r''g''y'];fork=1:4z=2.*(rand(3,n)<=p)-1;x=[zeros(1,3);cumsum(z')];col=colorstr(k);plot3(x(:,1),x(:,2),x(:,3),col);holdonendgrid3、三維布朗運(yùn)動(dòng)bm3plot.mnpoints=5000;dt=1;bm=cumsum([zeros(1,3);dt^0.5*randn(npoints-1,3)]);figure(1);plot3(bm(:,1),bm(:,2),bm(:,3),'k');pcol=(bm-repmat(min(bm),npoints,1))./...repmat(max(bm)-min(bm),npoints,1);holdon;scatter3(bm(:,1),bm(:,2),bm(:,3),...10,pcol,'filled');gridon;holdoff;4、二維和三維空間中的泊松點(diǎn)poisson2d.m,poisson3d.m這是在空間中隨機(jī)、獨(dú)立地放置點(diǎn)的通用模型。在任何給定的空間集合中，將放置強(qiáng)度與其容量成比例的泊松分布的點(diǎn)數(shù)。在任意兩個(gè)不相交的集合中的點(diǎn)數(shù)是獨(dú)立的。%單位體積的泊松點(diǎn)數(shù)強(qiáng)度為lambdalambda=100;nmb=poissrnd(lambda)x=rand(1,nmb);y=rand(1,nmb);z=rand(1,nmb);gridscatter3(x,y,z,5,5.*rand(1,nmb));

第四章Markov過(guò)程Markov性是隨機(jī)序列x(1),x(2)...的一種特殊形式的依賴(lài)。如果我們知道過(guò)去x(1)...,x(n-1)和現(xiàn)在x(n)，這種信息可能會(huì)或者可能不會(huì)影響未來(lái)x(n+1),x(n+1)...。Markov過(guò)程(Markov鏈)沒(méi)有從過(guò)去獲得額外的信息。對(duì)于未來(lái)的一切都由我們現(xiàn)在的信息所決定。隨機(jī)游動(dòng)和泊松過(guò)程是特殊的簡(jiǎn)單Markov過(guò)程的例子。對(duì)于那些模擬，我們已經(jīng)使用了他們結(jié)構(gòu)中內(nèi)蘊(yùn)的特殊的獨(dú)立性質(zhì)。在這個(gè)部分討論的Matlab文件有：simgeod1.m,simmm1.m,simmd1.m,simmg1.m,simmginfty.m,simstmginfty.m,birthdeath.m,moran.m,galtonwatson.m一、離散服務(wù)系統(tǒng)中的緩沖動(dòng)力學(xué)simgeod1.m假設(shè)時(shí)間是離散的，并且顧客按照一個(gè)獨(dú)立序列a(1),a(2)...到達(dá)服務(wù)中心，其中a(k)是在第k期到達(dá)的數(shù)量。一名顧客被服務(wù)一期(單服務(wù)系統(tǒng))。其他的顧客在一個(gè)緩沖區(qū)域等候，直到可以被服務(wù)。因此，在k時(shí)刻系統(tǒng)中的顧客數(shù)量為n(k)=n(k-1)+a(k)-I{n(k-1)+a(k)>=1},k>=2。加上初始條件n(1)=0，遞歸定義一個(gè)Markov鏈n(k)，k>=1.試試參數(shù)p的不同取值。從長(zhǎng)遠(yuǎn)看，會(huì)發(fā)生什么呢?m=200;p=0.2;N=zeros(1,m);%初始化緩沖區(qū)A=geornd(1-p,1,m);%生成到達(dá)序列模型,比如，幾何分布forn=2:mN(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);endstairs((0:m-1),N);二、M/M/1模型simmm1.m這是一個(gè)連續(xù)時(shí)間的單服務(wù)緩沖模型。系統(tǒng)的到達(dá)由強(qiáng)度為lambda的泊松過(guò)程決定。服務(wù)員為每位顧客的服務(wù)時(shí)間服從指數(shù)分布，均值是1/mu。由此得到的系統(tǒng)規(guī)模N(t),t>=0，是一個(gè)連續(xù)時(shí)間的Markov過(guò)程，其演變?nèi)缦隆腘(0)=n_0開(kāi)始。等待強(qiáng)度為lambda+mu的指數(shù)分布時(shí)間(如果n_0=0，強(qiáng)度為lambda),然后以可能性lambda/(lambda+mu)向前跳躍和以可能性mu/(lambda+mu)向后跳躍。如此循環(huán)。在N=0時(shí)動(dòng)力學(xué)的改變是在開(kāi)始用一個(gè)短循環(huán)來(lái)實(shí)現(xiàn)的：ifi==0mutemp=0;elsemutemp=mu;end主循環(huán)僅僅檢查向前跳或者向后跳：ifrand<=lambda/(lambda+mutemp)i=i+1;%&向前跳：一個(gè)客戶(hù)到達(dá)elsei=i-1;%向后跳：一個(gè)客戶(hù)離開(kāi)endx(k)=i;%在i時(shí)刻的系統(tǒng)大小有一個(gè)避免所有循環(huán)的方法，見(jiàn)下面的M/G/1系統(tǒng)。三、M/D/1系統(tǒng)simmd1.m與M/M/1一樣，這個(gè)系統(tǒng)的到達(dá)為泊松過(guò)程，但每次服務(wù)時(shí)間是固定的長(zhǎng)度1(例如，在緩沖環(huán)節(jié)中，固定大小的數(shù)據(jù)包的傳輸時(shí)間)。這不是Markov過(guò)程?？梢宰C明,顧客離開(kāi)這個(gè)系統(tǒng)發(fā)生在u_k=k+max(t_1,t_2-2+1..,t_k-k+1)時(shí)刻,其中t_1,t_2...是泊松過(guò)程的到達(dá)時(shí)刻。因此，系統(tǒng)規(guī)模過(guò)程N(yùn)(t),t>=0,在t_k時(shí)向前跳躍，在u_k時(shí)向后跳躍。假設(shè)我們有長(zhǎng)度為n到達(dá)時(shí)向量t_k。這里是得到離開(kāi)時(shí)間的一種方法：arrsubtr=arrtime-(0:n-1)';%t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptime=(1:n)+max(triu(arrmatrix))現(xiàn)在想畫(huà)出N(t):B=[ones(n,1)arrtime;-ones(n,1)deptime'];Bsort=sortrows(B,2);%按次序?qū)⑻S分類(lèi)jumps=Bsort(:,1);jumptimes=[0;Bsort(:,2)];X=[0;cumsum(jumps)];四、M/G/1系統(tǒng)simmg1.m這是將M/D/1推廣到一般服務(wù)時(shí)間分布S，其均值為1/mu。一個(gè)相似的計(jì)算離開(kāi)時(shí)間的技術(shù)在此情形也是可以的。特別地,取Exp(mu)分布的服務(wù)時(shí)間就是M/M/1系統(tǒng)。因而，我們有了模擬M/M/1的另一種方法。如果lambda/mu<1，則系統(tǒng)處于穩(wěn)定狀態(tài)。五、M/G/infinity系統(tǒng)simmginfty.m這里每個(gè)顧客得到他自己的服務(wù)。沒(méi)有排隊(duì)。模擬比M/G/1簡(jiǎn)單。產(chǎn)生到達(dá)時(shí)間加上服務(wù)時(shí)間。然后，就像上面的M/D/1系統(tǒng)一樣，不管系統(tǒng)規(guī)模的改變時(shí)間為何時(shí)，總標(biāo)記+1或-1。在示例代碼中，將演示如何得到Pareto分布的服務(wù)時(shí)間：alpha=1.5;%Pareto服務(wù)時(shí)間servtimes=rand^(-1/(alpha-1))-1;%平穩(wěn)更新過(guò)程servtimes=[servtimes;rand(npoints-1,1).^(-1/alpha)-1];六、M/G/infinity系統(tǒng),平穩(wěn)情形,任意服務(wù)時(shí)間simstmginfty.m(依賴(lài)distrmu.m,distrstat.m)當(dāng)我們從時(shí)間0觀察一個(gè)平穩(wěn)系統(tǒng)時(shí)，它在負(fù)時(shí)間一直是活躍到"永遠(yuǎn)"的。因此，在時(shí)刻0，服務(wù)時(shí)間服從平穩(wěn)分布的系統(tǒng)中有參數(shù)為(lambda*mu)的泊松分布個(gè)顧客。這在simstmginfty.m中得以實(shí)現(xiàn)，它是simmginfty.m的一個(gè)修改版本。simstmginfty.m也允許用一個(gè)服從服務(wù)時(shí)間分布的隨機(jī)數(shù)生成器作為輸入?yún)?shù)。參數(shù)是帶一個(gè)適當(dāng)?shù)耐獠亢瘮?shù)和含有分布參數(shù)數(shù)組的MATLAB函數(shù)句柄,參見(jiàn)作為外部參數(shù)的RNG。例子．產(chǎn)生平穩(wěn)M/G/infinity隊(duì)列中[0,5)時(shí)間內(nèi)的系統(tǒng)規(guī)模過(guò)程，到達(dá)強(qiáng)度為lambda=2，服務(wù)時(shí)間服從alpha=1.6,gamma=2的標(biāo)準(zhǔn)Pareto分布。[jmptimes,syssize]=simstmginfty(5,2,@simparetonrm,{1.6,2},1);stairs(jmptimes,syssize);加入simmginfty.m得到平穩(wěn)版本的步驟：1、產(chǎn)生在0時(shí)刻的"負(fù)"到達(dá)和它們的平穩(wěn)服務(wù)時(shí)間。用表查找函數(shù)distrstat.m來(lái)得到平穩(wěn)分布隨機(jī)數(shù)生成器的句柄。2、%在時(shí)刻0，平穩(wěn)服務(wù)時(shí)間的系統(tǒng)中有泊松顧客數(shù)。3、nstart=poissrnd(lambda*servmu);%泊松隨機(jī)變量4、if(nstart>0)5、[statdist,statpar]=distrstat(servdist,servpar);%平穩(wěn)分布句柄6、rndpar1={nstart,1,statpar{:}};%隨機(jī)數(shù)生成器參數(shù)7、stattimes=feval(statdist,rndpar1{:});%平穩(wěn)服務(wù)時(shí)間8、arrtimes=zeros(size(stattimes));%在t＝0時(shí)刻前到達(dá)的顧客按到達(dá)時(shí)間為0來(lái)看待。9、end10、一旦創(chuàng)建了計(jì)數(shù)過(guò)程,就刪除開(kāi)始時(shí)"負(fù)"到達(dá)額外的0點(diǎn)。增加maxtime到跳躍，以得到正確的圖形。七、生滅過(guò)程1、一般的生滅強(qiáng)度birthdeath.m作為例子，我們選擇在水平i上出生強(qiáng)度為lambda_i=lambda/(1+i)，死亡強(qiáng)度為mu_i=mu*i的模型,其中l(wèi)ambda和mu為固定的常數(shù)。要求循環(huán)滿足直到下次跳躍，跳躍強(qiáng)度和等待時(shí)間才被更新,即lambda_i=lambda/(1+i);ifi==0mu_i=0;elsemu_i=mu*i;endtime=-log(rand)./(lambda_i+mu_i);2、Moran模型moran.m另一個(gè)起源于遺傳學(xué)的生滅過(guò)程。有限狀態(tài)空間,吸收界限。x=1:N+1;lambda(x)=(x-1).*(1-(x-1)./N);%出生率。mu(x)=(x-1).*(1-(x-1)./N);%死亡率。q(x)=lambda(x)+mu(x);%兩次跳躍時(shí)間間隔的指數(shù)分布率。八、分支過(guò)程Galton-Watson過(guò)程galtonwatson.m離散時(shí)間，生命長(zhǎng)度為1。死亡的每個(gè)個(gè)體產(chǎn)生隨機(jī)的后代個(gè)數(shù)。函數(shù)offspring(k)給出從人口規(guī)模為k開(kāi)始的祖先向量。p=[1/201/2];z=[cumsum(p)];n=length(p);%可能的子孫數(shù)量offmu=dot(0:n-1,p);%子孫的平均個(gè)數(shù)u1=sort(rand(1,k));forj=1:nu(j)=length(find(u1<z(j)));endu=diff([0u]);nu=u*(0:n-1)';九、計(jì)數(shù)過(guò)程計(jì)數(shù)過(guò)程N(yùn)(t)記錄的是實(shí)值隨機(jī)點(diǎn)過(guò)程{T_k}在區(qū)間[0,t)內(nèi)的點(diǎn)數(shù)。泊松過(guò)程及排隊(duì)系統(tǒng)中遇到的系統(tǒng)規(guī)模過(guò)程都是計(jì)數(shù)過(guò)程。十、更新過(guò)程更新過(guò)程，是一列獨(dú)立同分布的正隨機(jī)變量的部分和序列。這個(gè)過(guò)程可以被想象為：當(dāng)同種生物的一些個(gè)體生命期結(jié)束，同時(shí)他們也被新的生命所代替的時(shí)間點(diǎn)序列。更新計(jì)數(shù)過(guò)程記錄的是在時(shí)間區(qū)間[0，t）內(nèi)更新的次數(shù)。它是一個(gè)隨機(jī)階梯函數(shù)。

第五章模擬的應(yīng)用和例子大數(shù)定律表明：1、經(jīng)驗(yàn)均值收斂到它的期望值；2、統(tǒng)計(jì)物理中，軌道平均與總體平均是漸進(jìn)相同的；3、為隨機(jī)模擬提供了理論基礎(chǔ)，并建立了事件頻率和概率的聯(lián)系。一、計(jì)算積分I=\int_a^bf(x)dx1、I=(b-a)E(f(a+(b-a)U))，U\sim(0,1)。模擬X_j=E(f(a+(b-a)U_j))，用平均1/n\sum_{j=1}^nX_j逼近I/(b-a)。大數(shù)定律說(shuō)明，可以用獨(dú)立試驗(yàn)的頻率來(lái)近似期望值。2、選取一個(gè)包含函數(shù)f圖形的矩形，比如[a,b][min{f},max{f}]。生成該矩形上n對(duì)均勻分布隨機(jī)數(shù)(X,Y),記錄事件“Y<f(X)”發(fā)生的個(gè)數(shù)m。利用該事件發(fā)生的比例m/n來(lái)近似I。這利用了積分的幾何含義、幾何型概率、以及大數(shù)定律。3、I=\int_a^bf(x)dx=\int_a^b\frac{f(x)}{g(x)}g(x)dx=E(\frac{f(X)}{g(X)})，然后用1/N\sum_{j=1}^Nf(X_j)/g(X_j)來(lái)近似I，其中X_j為獨(dú)立的密度函數(shù)為g(x)的隨機(jī)數(shù)，\int_a^bg(x)dx=1。通過(guò)選擇合適的函數(shù)g可以減小方差,這被稱(chēng)為重要樣本法。這是MonteCarlo方法的基礎(chǔ)，它是一類(lèi)計(jì)算積分的概率方法。n次近似的方差的階是n^{-1/2}，這比光滑函數(shù)時(shí)的梯度法差些。但作為回報(bào)，該近似對(duì)維數(shù)、被積函數(shù)的光滑性不敏感。因此MonteCarlo方法應(yīng)用于不正則區(qū)域上的多重積分非常有效。另外，MonteCarlo方法的可靠性、誤差的上界依賴(lài)隨機(jī)數(shù)生成器的質(zhì)量。在需要大量隨機(jī)化的問(wèn)題中使用不知道的隨機(jī)數(shù)生成器是很草率的。二、誤差估計(jì)1、Chebyshev不等式。由P(|X-EX|>t)<\frac{1}{t^2}Var(X)，得P(|1/n\sum_{j=1}^nX_j-EX|>t)<\frac{1}{nt^2}Var(X).然而，Var(X)往往也是不知道的，這是由1/n\sum_{j=1}^n(X_j-\bar{X})^2來(lái)近似的，其中\(zhòng)bar{X}=1/n\sum_{j=1}^nX_j。2、中心極限定理。由中心極限定理知道，近似地有\(zhòng)frac{1/n\sum_{j=1}^nX_j-EX}{Var(X)/n}服從標(biāo)準(zhǔn)正態(tài)。因此，P(|1/n\sum_{j=1}^nX_j-EX|<t)=2\phi(\frac{t}{{Var(X)/n}})-1=q。即以q的概率誤差小于\frac{t}{{Var(X)/n}}。當(dāng)然，Var(X)也如上面一樣來(lái)近似。中心極限定理的另一個(gè)用處是模擬某種連續(xù)過(guò)程的軌道，即X_n(t)=\frac{1}{\sqrt{n}}\sum_{k<=nt}\xi_k，其中\(zhòng)xi_k為獨(dú)立的0均值隨機(jī)變量。當(dāng)然，它還可以用來(lái)近似模擬標(biāo)準(zhǔn)正態(tài)隨機(jī)數(shù)。三、減小方差技術(shù)（對(duì)立變量、條件期望、控制變量、重要樣本）上面提到，方差往往是不知道,它也是通過(guò)產(chǎn)生的隨機(jī)數(shù)來(lái)估計(jì)。由估計(jì)的誤差分析知,該方差越小越好，因此給出幾種減小方差的一般方法?！摹⒛M的例子（一）概率問(wèn)題的模擬問(wèn)題一車(chē)和羊的游戲；問(wèn)題二蒲豐投針問(wèn)題；問(wèn)題三擲骰子問(wèn)題；問(wèn)題四無(wú)記憶性的例子；問(wèn)題五生日問(wèn)題；問(wèn)題六Galton釘板實(shí)驗(yàn)；問(wèn)題七趕火車(chē)問(wèn)題；問(wèn)題一車(chē)和羊的游戲假設(shè)你在進(jìn)行一個(gè)游戲節(jié)目。現(xiàn)給三扇門(mén)供你選擇：一扇門(mén)后面是一輛轎車(chē)，另兩扇門(mén)后面分別都是一頭山羊。你的目的當(dāng)然是要想得到比較值錢(qián)的轎車(chē)，但你卻并不能看到門(mén)后面的真實(shí)情況。主持人先讓你作第一次選擇。在你選擇了一扇門(mén)后，知道其余兩扇門(mén)后面是什么的主持人，打開(kāi)了另一扇門(mén)給你看，而且，當(dāng)然，那里有一頭山羊?，F(xiàn)在主持人告訴你，你還有一次選擇的機(jī)會(huì)。那么，請(qǐng)你考慮一下，你是堅(jiān)持第一次的選擇不變，還是改變第一次的選擇，更有可能得到轎車(chē)？《廣場(chǎng)雜志》刊登出這個(gè)題目后，竟引起全美大學(xué)生的舉國(guó)辯論，許多大學(xué)的教授們也參與了進(jìn)來(lái)。真可謂盛況空前。據(jù)《紐約時(shí)報(bào)》報(bào)道，這個(gè)問(wèn)題也在中央情報(bào)局的辦公室內(nèi)和波斯灣飛機(jī)駕駛員的營(yíng)房里引起了爭(zhēng)論，它還被麻省理工學(xué)院的數(shù)學(xué)家們和新墨哥州洛斯阿拉莫斯實(shí)驗(yàn)室的計(jì)算機(jī)程序員們進(jìn)行過(guò)分析。問(wèn)題分析在一次實(shí)驗(yàn)中，如果第一次選擇選中了轎車(chē)（概率為1/3），那么主持人打開(kāi)一扇門(mén)后，如果堅(jiān)持原來(lái)的選擇，則能得到轎車(chē)，反之，改變第一次選擇則不能得到轎車(chē)。如果第一次沒(méi)有選中轎車(chē)（概率為2/3），那么其余兩扇門(mén)后面必有一個(gè)是轎車(chē)，主持人只能打開(kāi)有山羊有那扇門(mén)，則剩下的一扇門(mén)后面是轎車(chē)，此時(shí)堅(jiān)持原來(lái)的選擇不能得到轎車(chē)，改變第一次的選擇必能得到轎車(chē)。因此，經(jīng)過(guò)分析，堅(jiān)持第一次的選擇不變得到轎車(chē)的概率為1/3，改變第一次的選擇得到轎車(chē)的概率為2/3。實(shí)際上，在只有三扇門(mén)的情況下，那么改不改變選擇效果并不明顯。如果有100扇門(mén)，參與的嘉賓選擇了其中的一扇，而主持人隨后把剩下的99扇門(mén)中間的98扇門(mén)都打開(kāi)，這98扇門(mén)后面都沒(méi)有獎(jiǎng)品，這時(shí)應(yīng)該改變選擇，畢竟最開(kāi)始自己選擇的那扇門(mén)中獎(jiǎng)的概率只是1%而已。需要注意的是，主持人是在知道其他兩扇門(mén)后面都有什么的情況下選擇一個(gè)門(mén)打開(kāi)的。這種情況下三個(gè)門(mén)后是轎車(chē)的概率因?yàn)橹鞒秩酥澜Y(jié)果并參與其中而關(guān)聯(lián)在一起，而不是孤立等同的。如果打開(kāi)門(mén)的不是主持人，而是另一個(gè)參與者，并且當(dāng)他打開(kāi)門(mén)時(shí)發(fā)現(xiàn)什么也沒(méi)有，那么，剩下的兩個(gè)門(mén)后是轎車(chē)的概率才是相等的。計(jì)算機(jī)模擬為了驗(yàn)證這一結(jié)果，我們就要比較不改變選擇中獎(jiǎng)的幾率和改變選擇中獎(jiǎng)的幾率。模擬方法是：我們從0，1，2這3個(gè)數(shù)中隨機(jī)一個(gè)為轎車(chē)（即中獎(jiǎng)號(hào)碼），另隨機(jī)一個(gè)數(shù)為你的選擇。如果你的選擇與中獎(jiǎng)號(hào)相同，則計(jì)這次為不改變選擇中獎(jiǎng)；如果你的選擇不對(duì)，則是改變選擇中獎(jiǎng)。分別累積出不改變選擇中獎(jiǎng)和改變選擇中獎(jiǎng)的次數(shù)，就可以得到不改變選擇中獎(jiǎng)的幾率和改變選擇中獎(jiǎng)的幾率了。為了將結(jié)果表示的明顯，我們可以假設(shè)有100扇門(mén)，參與的嘉賓選擇了其中的一扇，而主持人隨后把剩下的99扇門(mén)中間的98扇門(mén)都打開(kāi)，這98扇門(mén)后面都沒(méi)有獎(jiǎng)品，然后模擬并比較不改變選擇中獎(jiǎng)的幾率和改變選擇中獎(jiǎng)的幾率。此時(shí)的情況也是相同的，只是每次隨即都是從0到99中隨機(jī)數(shù)而已。結(jié)果及分析下面兩幅圖分別是3個(gè)門(mén)時(shí)不改變選擇中獎(jiǎng)的概率在N次模擬結(jié)果下的概率分布（第二幅是為了便于觀察特意畫(huà)在固定坐標(biāo)軸上的）。下面則是100個(gè)門(mén)的情況下，不改變選擇中獎(jiǎng)的概率分布：可以顯然地看出在主持人幫助的情況下，改變選擇是可以大大增加自己中獎(jiǎng)的幾率的。通過(guò)這樣一個(gè)例子，我并不是想說(shuō)明什么概率意義上的問(wèn)題。只是通過(guò)這么一個(gè)模擬過(guò)程來(lái)學(xué)習(xí)計(jì)算機(jī)隨機(jī)模擬的一些基本方法與技巧。像在主持人不知道內(nèi)幕隨機(jī)的打開(kāi)一個(gè)發(fā)現(xiàn)是山羊這，我們可以通過(guò)同樣的隨機(jī)模擬過(guò)程來(lái)模擬這種情況。并可以驗(yàn)證改變選擇與否對(duì)自己中獎(jiǎng)的影響是相同的。當(dāng)模擬的次數(shù)逐漸的增多時(shí)，其模擬值越接近理論值，這說(shuō)明模擬的效果越好。在隨機(jī)事件的大量重復(fù)出現(xiàn)中，往往呈現(xiàn)幾乎必然的規(guī)律，這個(gè)規(guī)律就是大數(shù)定律。通俗地說(shuō)，這個(gè)定理就是，在試驗(yàn)不變的條件下，重復(fù)試驗(yàn)多次，隨機(jī)事件的頻率近似于它的概率。偶然必然中包含著必然。此次模擬試驗(yàn)也正好用實(shí)際的模擬例子說(shuō)明了大數(shù)定理的正確性和應(yīng)用性。Matlab程序1、編寫(xiě)函數(shù)n=10000;%實(shí)驗(yàn)次數(shù)stick=0;%堅(jiān)持選擇的獲獎(jiǎng)次數(shù)alter=0;%改變選擇的獲獎(jiǎng)次數(shù)fori=1:nx=floor(3*rand(1)+1);%在第x扇門(mén)后有轎車(chē)y=floor(3*rand(1)+1);%選擇第y扇門(mén)if(x==y)%第一次選中的情況stick=stick+1;else%第一次沒(méi)選中的情況alter=alter+1;endendstick=stick/n%轉(zhuǎn)化為百分比alter=alter/n保存為myone.m，在commandwindow輸入myone,結(jié)果為stick=0.3349alter=0.66512、moni.m這是模擬不同次數(shù)并繪圖的程序functionmoni()N=[10:100:10000];fori=1:length(N)bubian(i)=bubianmoni(N(i));endgridonaxis([0,10000,0,1])holdon(這3行為固定坐標(biāo)軸和柵格的)plot(N,bubian)3、bubianmoni.m這是模擬并返回不改變中獎(jiǎng)的概率的程序function[x]=bubianmoni(n)A=0;B=0;fori=1:na=randint(1,1,[0,99]);b=randint(1,1,[0,99]);ifa==bA=A+1;%不變的累加elseB=B+1;%改變的累加endendx=A/n;%不變中獎(jiǎng)的概率其中對(duì)改變選擇的累計(jì)在此處是可以加以省略的。問(wèn)題二蒲豐投針問(wèn)題 1777年法國(guó)科學(xué)家蒲豐提出的一種計(jì)算圓周率的方法隨機(jī)投針?lè)?，這種方法的具體步驟是：1、取一張白紙，在上面畫(huà)上許多條間距為d的平行直線。2、取一根長(zhǎng)為l(l<d)的針，隨機(jī)向畫(huà)有平行直線的紙上擲n次，觀察針與直線相交的次數(shù)，記為m。3、計(jì)算針與直線相交的概率。問(wèn)題分析針與平行線相交的充要條件是：，其中 建立直角坐標(biāo)系(，x)，上述條件在坐標(biāo)系下將是曲線所圍成的曲邊梯形區(qū)域：由幾何概率知：經(jīng)統(tǒng)計(jì)實(shí)驗(yàn)估計(jì)出概率，有上式即，故。Matlab程序n=input('請(qǐng)輸入投針的次數(shù)n:')l=1;d=2;m=0;fork=1:n%模擬n次x=unifrnd(0,d/2);%x為0到d/2之間的隨機(jī)變量，服從均勻分布y=unifrnd(0,pi);%y為0到pi之間的隨機(jī)變量，服從均勻分布ifx<0.5*l*sin(y);m=m+1;elseendendp=m/n;pi_m=1/p;fprintf('所模擬出的圓周率為：%f\n',pi_m)運(yùn)行結(jié)果與分析運(yùn)行了4次程序，并賦予不同的n值，得如下結(jié)果：n=100 所模擬出的圓周率為：3.225806n=1000 所模擬出的圓周率為：3.257329n=10000 所模擬出的圓周率為：3.176620n=100000 所模擬出的圓周率為：3.136664可見(jiàn)，模擬的次數(shù)越多，所得的結(jié)果就越準(zhǔn)確。這與理論所知的結(jié)果是一樣的。問(wèn)題三擲骰子問(wèn)題 擲三只相同的色子，求有兩個(gè)點(diǎn)數(shù)相同的概率。問(wèn)題分析 構(gòu)造3個(gè)整數(shù)隨機(jī)變量，使之服從1到6之間的均勻分布，比較每次模擬出的3個(gè)值，如果至少有兩個(gè)相同，則將成功的次數(shù)加1，分析可知，成功的頻率即為所求的概率。Matlab程序n=input('請(qǐng)輸入模擬的次數(shù):')l=0;fori=1:nx=unidrnd(6,1,3);%x中含有3個(gè)元素，相互獨(dú)立，且服從0到6的整數(shù)均勻分布ifx(1)==x(2)|x(2)==x(3)|x(1)==x(3)l=l+1;%比較模擬出的結(jié)果，l為成功數(shù)end;end;p=l/n;%計(jì)算成功的頻率，即為概率fprintf('至少兩個(gè)相同的概率是%f\n',p)運(yùn)行結(jié)果與分析 以下是取4個(gè)不同模擬次數(shù)所得的結(jié)果：n=100 至少兩個(gè)相同的概率是0.470000n=1000 至少兩個(gè)相同的概率是0.456000n=10000 至少兩個(gè)相同的概率是0.437300n=100000 至少兩個(gè)相同的概率是0.444390由概率知識(shí)可知：P=1-(4*5*6)/(6*6*6)=5/9=0.4445。這與以上模擬出的值基本上是一樣的。問(wèn)題四無(wú)記憶性的例子考慮有兩個(gè)營(yíng)業(yè)員的郵局。假設(shè)當(dāng)Smith先生進(jìn)入郵局的時(shí)候，他看到兩個(gè)營(yíng)業(yè)員分別在為Jones先生和Brown先生服務(wù)，并且被告知先服務(wù)完的營(yíng)業(yè)員即將為他服務(wù)。再假設(shè)為每顧客服務(wù)的時(shí)間都服從參數(shù)為.的指數(shù)分布。1).求這三個(gè)顧客中Smith最后離開(kāi)郵局的概率？Jones和Brown呢？2).求Smith在郵局花費(fèi)的平均時(shí)間E[T]？問(wèn)題分析當(dāng)Smith接受服務(wù)的時(shí)候，Jones和Brown只有一個(gè)人離開(kāi)，另一人仍在接受服務(wù)。然而，由指數(shù)分布的無(wú)記憶性，從現(xiàn)在起剩下的這個(gè)人被服務(wù)的時(shí)間是參數(shù)為.的指數(shù)分布，即與Smith是同分布的。再由對(duì)稱(chēng)性，他與Smith先離開(kāi)的概率一樣，都為1/2。第1問(wèn)答案為1/2,1/4,1/4。E[T]=E[minX;Y+S]=1/(2*λ)+1/λ,其中S為Smith服務(wù)花費(fèi)的時(shí)間,X;Y分別為Jones和Brown服務(wù)所花費(fèi)的時(shí)間。Matlab程序編寫(xiě)函數(shù)functionmytwo(lambda)n=100000;%實(shí)驗(yàn)次數(shù)jones=0;%三個(gè)變量記錄各人最后離開(kāi)次數(shù)brown=0;smith=0;sumT=0;%統(tǒng)計(jì)n次實(shí)驗(yàn)smith花費(fèi)時(shí)間總和fori=1:nt1=-log(rand)/lambda;%t1,t2,t3分別為三人服務(wù)時(shí)間t2=-log(rand)/lambda;t3=-log(rand)/lambda;if(t1+t3<t2)%brown最后離開(kāi)brown=brown+1;t=t1+t3;elseif(t2+t3<t1)%jones最后離開(kāi)jones=jones+1;t=t2+t3;else%smith最后離開(kāi)smith=smith+1;t=t3+(t1<t2)*t1+(t1>=t2)*t2;%t=t3+min(t1,t2)endsumT=sumT+t;endsmith=smith/n%轉(zhuǎn)化為概率jones=jones/nbrown=brown/nET=sumT/n%smith平均花費(fèi)時(shí)間ETreal=3/(2*lambda)%分析出的E[T],用來(lái)與模擬結(jié)果對(duì)比保存為mytwo.m，在commandwindow輸入mytwo(3),結(jié)果為smith=0.5019jones=0.2495brown=0.2486ET=0.5009ETreal=0.5000問(wèn)題五生日問(wèn)題有n個(gè)人的班級(jí)中至少有兩個(gè)人生日是同一天的概率，模擬之。問(wèn)題分析假設(shè)：所有人生日是相互獨(dú)立的。設(shè)一年有y天，為求至少兩個(gè)人生日在同一天的概率，我們可以先求n個(gè)人生日全不同的概率：n個(gè)人生日所有可能的組合（不考慮次序）種數(shù)有：而n個(gè)人生日都不一樣的可能的組合方式種數(shù)有：故任意兩人生日不同的理論概論是：（這樣寫(xiě)可以方便編程計(jì)算）從而有至少兩人生日同一天的概率是注：顯然當(dāng)學(xué)生人數(shù)n比一年的天數(shù)y大時(shí)，有抽屜原理易知：p=1.編程思想主程序：在程序初始位置可以定義學(xué)生人數(shù)n（默認(rèn)為50），隨機(jī)模擬次數(shù)tn（默認(rèn)為100），一年的天數(shù)y（默認(rèn)為365）。通過(guò)調(diào)用子函數(shù)計(jì)算隨機(jī)模擬的概率以及理論概率，將二者畫(huà)在同一個(gè)圖中，便于對(duì)比分析。子程序1：func1.m計(jì)算隨機(jī)試驗(yàn)頻率。在每一次隨機(jī)試驗(yàn)中，首先隨機(jī)生成1到y(tǒng)的一個(gè)隨機(jī)數(shù)，對(duì)應(yīng)于一個(gè)人的生日，根據(jù)人數(shù)存放于數(shù)組a中。使用unique指令，將a的相同項(xiàng)合并,并按照從小到大順序排列。得到的新數(shù)組與原來(lái)的數(shù)組a進(jìn)行長(zhǎng)度的比較，如果新數(shù)組的長(zhǎng)度小于原數(shù)組，說(shuō)明有兩人或者以上同一天生日，此時(shí)成功次數(shù)加1。總共經(jīng)過(guò)tn次隨機(jī)試驗(yàn)，即可求出至少兩個(gè)人同一天生日的試驗(yàn)的頻率。子程序2：func2.m計(jì)算理論概率通過(guò)上面的分析表達(dá)式分子是一個(gè)1*y的向量的連乘，分母也是一個(gè)1*y維向量的連乘。對(duì)應(yīng)分量做除法后求新向量個(gè)元素乘積即可。P就可以通過(guò)p=1-p0來(lái)計(jì)算了。結(jié)果分析圖一程序運(yùn)行所花費(fèi)時(shí)間為4.380000e-001秒圖二程序運(yùn)行所花費(fèi)時(shí)間為8.750000e-001秒圖三程序運(yùn)行所花費(fèi)時(shí)間為6.688000e+000秒圖四程序運(yùn)行所花費(fèi)時(shí)間為2.438000e+000秒圖五程序運(yùn)行所花費(fèi)時(shí)間為2.232800e+001秒圖六程序運(yùn)行所花費(fèi)時(shí)間為2.219210e+002秒（1）比較上面幾幅圖，對(duì)應(yīng)的學(xué)生人數(shù)均為50人。而隨機(jī)模擬的試驗(yàn)次數(shù)在增加，當(dāng)tn比較小的時(shí)候，如圖一tn=30，頻率線一直在理論概率線附近上下波動(dòng)，而且波動(dòng)比較大，十分顯著。當(dāng)試驗(yàn)次數(shù)逐漸增加的時(shí)候，頻率線也逐漸趨近于概率線，當(dāng)tn=1000時(shí)，兩者幾乎重合了。當(dāng)試驗(yàn)次數(shù)tn越大時(shí)，所得的頻率曲線就越精確接近于概率線，但是同時(shí)付出的時(shí)間代價(jià)也是巨大的。當(dāng)tn=1000時(shí)，運(yùn)行時(shí)間需要約三分鐘。（2）觀察每一幅圖，都有一個(gè)共同的特征，隨著人數(shù)的增加，至少兩個(gè)人同一天生日的概率也逐漸增大，而且通過(guò)概率線可以看出。當(dāng)人數(shù)為0和1時(shí)，p=0，當(dāng)人數(shù)達(dá)到40的時(shí)候，至少兩人同一天生日的概率就可以達(dá)到90%，當(dāng)人數(shù)達(dá)到60或者70時(shí)候，概率已經(jīng)達(dá)到99%以上，這是幾乎可以肯定的說(shuō)至少有兩個(gè)人是同一天生日的。（3）注：由于數(shù)字是隨機(jī)生成的，所以每次運(yùn)行的結(jié)果都不一樣，但是大體上的趨勢(shì)是如此的。另外，這個(gè)運(yùn)行時(shí)間也僅供對(duì)比參考，他還與及其運(yùn)行的狀態(tài)有關(guān)系，因此沒(méi)有太大的實(shí)際意義。Matlab程序程序說(shuō)明：主程序名suiji.m子程序1：func1.m子程序2：func2.m程序如下：【suiji.m】%肖延淵%2008年5月%每年的天數(shù)year%人數(shù)n[number]%隨機(jī)模擬次數(shù)tn[testnumber]clear;clc;c1=clock;n=50;%學(xué)生人數(shù)tn=100;%隨機(jī)模擬次數(shù)y=365;%一年的天數(shù),一般計(jì)算取365freq=[];%用以存放不同n對(duì)應(yīng)的頻率prob=[];%用以存放不同n對(duì)應(yīng)的概率fori=1:nfreq(i)=func1(i,tn,y);%求頻率prob(i)=func2(i,y);%求概率endplot(1:i,freq,1:i,prob);%畫(huà)圖title(['頻率線接近概率線隨機(jī)模擬次數(shù)',int2str(tn)]);legend('頻率線','概率線',2);xlabel('學(xué)生人數(shù)n');ylabel('頻率或者概率');c2=clock;c=c2-c1;fprintf('程序運(yùn)行所花費(fèi)時(shí)間為%d秒',c(4)*3600+c(5)*60+c(6));【func1.m】%求隨機(jī)試驗(yàn)頻率函數(shù)functionf=func1(n,tn,y)f=0;fort=1:tna=[];%存放一個(gè)人生日f(shuō)ork=1:n%向上取整，隨機(jī)取1到y(tǒng)的一個(gè)整數(shù)，對(duì)應(yīng)于一個(gè)人的生日b=ceil(y*rand(1));a=[a,b];endc=unique(a);%合并a中的相同的項(xiàng)，并按照從小到大排序iflength(a)~=length(c)%a和c的長(zhǎng)度不相等，說(shuō)明至少有兩人同一天生日f(shuō)=f+1;%頻數(shù)加一endendf=f/tn;%頻數(shù)轉(zhuǎn)換成頻率%fprintf('任兩人不在同一天生日的頻率為：%f\n',f);【func2.m】%求理論概率函數(shù)functionp=func2(n,y)if(n>y)p=1;%學(xué)生人數(shù)比y（一年的天數(shù)大時(shí)），p衡為1.elsep1=1:y;p2=[1:y-n,y*ones(1,n)];p=p1./p2;%所有人都不在同一天生日的概率p=1-prod(p);%至少兩人同一天生日的概率end%fprintf('至少兩人同一天生日的概率為：%f\n',p);問(wèn)題六Galton釘板實(shí)驗(yàn)我們每次向釘板投一個(gè)小球時(shí)，并不知道它會(huì)落在哪個(gè)格子里，但通過(guò)大量的重復(fù)試驗(yàn)可知，小球的堆積形狀成單蜂形狀，而且左右對(duì)稱(chēng)，大部分落在中間的格子里，兩頭的格子里較少。模擬之。實(shí)驗(yàn)步驟確定釘子的位置：將釘子的橫，縱坐標(biāo)存儲(chǔ)在兩個(gè)矩陣X和Y中。在Galton釘板實(shí)驗(yàn)中，小球每次碰到釘子是都具有兩中可能情況，設(shè)向右的概率p，向左的概率為q=1-p,這里p=0.5,表示向左和向右的機(jī)會(huì)是相同的。模擬過(guò)程中，需要用到隨機(jī)數(shù)發(fā)生器指令rand(m,n)用以表示小球碰到釘子時(shí)是向左還是向右，每次調(diào)用rand(m,n)的結(jié)果都會(huì)不同，如果想保持一致，可與rand(‘seed’,s)配合使用，其中s是一個(gè)正整數(shù)，例如【rand(‘seed’,1),u=rand(1,6)】{u=0.95010.23110.60680.48600.89130.7621}而再次運(yùn)行該指令時(shí)結(jié)果保持不變。除非重新設(shè)置種子seed的值。如【rand(‘seed’,2),u=rand(1,6)】{u=0.02580.92100.70080.19010.86730.4185}模擬小球堆積的情況，輸入扔球次數(shù)m(例如，m=50,100,等),計(jì)算落在第i個(gè)格子的小球數(shù)在總球數(shù)m中的比例，這樣當(dāng)模擬結(jié)束是，就得到了頻率，用頻率反映小球的堆積形狀。所用的動(dòng)畫(huà)指令如下Moviein(n):創(chuàng)建動(dòng)畫(huà)矩陣Getframe:拷貝動(dòng)畫(huà)矩陣；Movie(Mat,m):播放動(dòng)畫(huà)矩陣m次Matlab程序clear,clfm=100;n=5;y0=2;%設(shè)置參數(shù)ballnum=zeros(1,n+1);p=0.5;q=1-p;fori=n+1:-1:1%創(chuàng)建釘子坐標(biāo)x,y;x(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;forj=2:ix(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);endendmm=moviein(m)