隨機(jī)過程試驗(yàn)講義

1、隨機(jī)過程實(shí)驗(yàn)講義劉繼成華中科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院2011-2012年上半年為華中科技大學(xué)數(shù)學(xué)系本科生講授隨機(jī)過程課程參考資料 TOC o 1-5 h z 刖 日1 HYPERLINK l bookmark5 o Current Document 第一章Matlab簡(jiǎn)介2 HYPERLINK l bookmark11 o Current Document 第二章簡(jiǎn)單分布的模擬6 HYPERLINK l bookmark36 o Current Document 第三章基本隨機(jī)過程9 HYPERLINK l bookmark66 o Current Document 第四章Markov過程12 H

2、YPERLINK l bookmark93 o Current Document 第五章模擬的應(yīng)用和例子16附錄各章的原程序51參考文獻(xiàn)75若想檢驗(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ù)被稱為偽隨機(jī)數(shù)。盡管如此，

3、對(duì)于多數(shù)的具體問題這樣的隨機(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ǔ)句，沒有堆砌同一種變量的多種模擬方法。對(duì)于沒有列舉的隨機(jī)變量的模擬，以及有特 殊需求的讀者應(yīng)該由這些方法得到啟發(fā)，或者參考更詳細(xì)的其他文獻(xiàn)資料。第三章是基本隨機(jī) 過程的模擬。主要是簡(jiǎn)單獨(dú)立增量過程的模擬，多維的推廣是直接的。第四章是Mark

4、ov過程的模擬。包括服務(wù)系統(tǒng)，生滅過程、簡(jiǎn)單分支過程等。第五章是這些模擬的應(yīng)用。例如，計(jì)算概 率、估計(jì)積分、模擬現(xiàn)實(shí)、誤差估計(jì)，以及減小方差技術(shù)，特別給讀者提供了一些經(jīng)典問題的 模擬，通過這些問題的模擬將會(huì)更加牢固地掌握實(shí)際模擬的步驟。平穩(wěn)過程的模擬、以及利用 平穩(wěn)過程來預(yù)測(cè)的內(nèi)容并沒有包含在本講義之內(nèi)，但這絲毫不影響該內(nèi)容的重要性，這也是將 會(huì)增補(bǔ)進(jìn)來的主要內(nèi)容之一。希望讀者碰到類似的問題時(shí)能夠查閱相關(guān)資料解決之。各章的內(nèi)容包括了模擬的基本思路和Matlab代碼。源程序包展示了對(duì)各種隨機(jī)過程與隨機(jī)機(jī)制的有效模擬和可視化的基本技術(shù)，試圖強(qiáng)調(diào)matlab自然處理矩陣和向量的方法，目標(biāo)是為涉及應(yīng)用

5、隨機(jī)模擬的讀者在準(zhǔn)備自己的程序代碼時(shí)找到靈感和想法。建議讀者在了解了模擬的 基本方法之后就著手解決自己感興趣的實(shí)際問題。對(duì)實(shí)際具體問題的解決有助于更深刻理解模 擬的思想、也會(huì)在具體應(yīng)用中拓展現(xiàn)有的模擬方法。第一章 Matlab簡(jiǎn)介若你想在計(jì)算機(jī)上運(yùn)行 Matlab，點(diǎn)擊：開始/程序/MATLAB 6.5 ,這樣將會(huì)出現(xiàn)有三個(gè)窗口 的交互界面。如果你是初學(xué)者，可以先瀏覽一下 Matlab的指導(dǎo)材料，點(diǎn)擊：Help/ MATLABHelp,打開窗口左邊的“ MATLAB一節(jié)即可看到相關(guān)的內(nèi)容。就自己而言，我學(xué)習(xí) Matlab更喜歡的方式是：輸入并運(yùn)行一些命令、觀察出現(xiàn)的結(jié)果，然后查閱想了解的幫助文

6、件。這也正是本節(jié)的方法。在command window”窗口中顯示有提示符“”，在提示符后輸入下面的命令，按回車鍵即可運(yùn)行并顯示相應(yīng)的結(jié)果。當(dāng)然，不要輸入行號(hào)、也不必輸入后面的注釋。在這個(gè)部分討論的 Matlab 文件有:rando.m , vrando.m , show.m。、Matlab 初步1:2*92： sin3: format long4： sin2A999format short： x=sinMatlabMatlab當(dāng)作計(jì)算器用。僅顯示四位小數(shù)，但保存的更多! 顯示更多位小數(shù)。將計(jì)算結(jié)果存在變量中。x顯示x的值。x=rand(10,1) x是包含有10個(gè)0 , 1上均勻分布隨機(jī)數(shù)的

7、集合， 它是一個(gè)列向量或者是 10 x 1的矩陣。10：x+511 : 1000*x的每個(gè)分量都加5。 的每個(gè)分量都乘以1000。12： x=rand(10,7) 10 X 7的隨機(jī)數(shù)矩陣。 若想重復(fù)此命令或其他命令，按住向上的光標(biāo)鍵直至看到想重復(fù)的命令。13:x=rand(1000,1)x=rand(1000,1);helphelp elmat： help rand 18:x(1:20,1): help punct: max(x)： mean(x): sum(x): median(x): cumsum(x)25 :y=sort(rand(10,1)將1000換成更大的數(shù)試試。用分號(hào)“；”可以

8、不顯示結(jié)果。顯示標(biāo)準(zhǔn)的幫助列表。顯示關(guān)于初等矩陣的幫助， 包括命令“rand”, 直接顯示rand ”的幫助。取出x的第一列中的1-20個(gè)數(shù)。Matlab中關(guān)于標(biāo)點(diǎn)符號(hào)的用法。的中位數(shù)。的分量累計(jì)和向量。 由小到大排序后的向量。26 : hist(x)27：hist(x,30)28 :y=-log(x)作出x的直方圖。用30個(gè)方柱代替缺省的10個(gè)。 又x的分量取自然對(duì)數(shù)。29： hist(y,30)多數(shù)的y的分量只是接近0的，但有些是和6差不多大的，y中的數(shù)被稱為指30 :z=randn(1000,1);31 :hist(z,30)數(shù)分布隨機(jī)數(shù)。生成1000個(gè)標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。直方圖是鐘形的

9、。對(duì)大于1000的數(shù)試試結(jié)果。二、獲取更多幫助32:如果你想查找不會(huì)使用的命令，可以點(diǎn)擊： Help/ MATLAB Help ,打開左邊的“ MATLAB 節(jié)，選擇Functions - Categorical List ”即可。據(jù)我所知，這是尋求幫助的最好方法。三、畫出數(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)，且將其連成線。: hold on將下面的圖形加到上面的圖形中。: plot(cumsum(x- 0.5), r);將這個(gè)結(jié)果圖畫到上面的圖形中。r表示用

10、紅色的線繪出，而缺省的顏色為藍(lán)色。： zoom on用鼠標(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);生成依賴z的隨機(jī)數(shù)。: plot(z,w, *);作出(z,w)的圖形。: axis(-3 3 -4 4);顯示 x在-3 , 3與y在-4 , 4范圍的圖形。四、作圖函數(shù): clf： ezplot( sin(x) ,0 3*pi);畫出正弦函數(shù)的圖形。: hold on46:t=0:0.01:3*pi;定義一個(gè)時(shí)間點(diǎn)向量，間隔為0.01。47 :tt為一行向

11、量。48:t=t 現(xiàn)在t為一列向量。: plot(t,sin(5*t), r);用紅色畫sin(t) 關(guān)于t的函數(shù)。顯然，函數(shù)ezplot不能設(shè)置圖形的顏色。:title( sin(x) and sin(5x)給圖形加上更恰當(dāng)?shù)臉?biāo)題。五、運(yùn)行現(xiàn)有的Matlab程序:上網(wǎng)下載或者拷貝一些編輯好的Matlab程序到自己的電腦中。：如果在你電腦的某個(gè)文件夾中有現(xiàn)成的Matlab程序(*.m)，可以設(shè)置Current Directory (Command WindoWB 口的上面)為該文件夾即可運(yùn)行這些程序。53:如果在你電腦里的幾個(gè)文件夾里都有Matlab程序，點(diǎn)擊菜單中：File/ Set Pat

12、h/Add Folder,加入所有這些文件夾，最后選擇Save。當(dāng)你在CommandWindo喇口鍵入一命令后，Matlab會(huì)在所有的這些文件夾中查找這個(gè)命令名。六、拋硬幣54：35不等式滿足結(jié)果為1。55:50.5Xx的所有分量檢查該不等式。58 : z=1+(x0.5)z的值為1或者2。這有點(diǎn)像拋硬幣，1為正面，2為反面。59 : show(z,正反)這是一個(gè)名字為show的程序，有兩個(gè)變量，一個(gè)是自然數(shù)向量，一個(gè)是用來與每個(gè)數(shù)相對(duì)應(yīng)顯示的字符串。它是自己編制的程序，保存在：d:MATLAB6p5workshow.m 。60: show(1+(rand(1500,1)0.5),正反)生成1

13、500個(gè)拋硬幣的結(jié)果。 現(xiàn)在按下向上的光標(biāo)鍵/回車，就會(huì)得到很多拋硬幣的結(jié)果。你找到連續(xù)出現(xiàn)正面最多的個(gè)數(shù)了嗎？61 ：show(1+(rand(1500,1)0.5),0 -)可以通過改變顯示的字符來簡(jiǎn)化剛才的問題。用向上的光標(biāo)鍵很容易更改前面的命令來實(shí)現(xiàn)它。這些語(yǔ)句對(duì)拋硬幣的問題當(dāng)然是足夠了，因?yàn)樗挥袃蓚€(gè)結(jié)果。但對(duì)其他，像擲色子，的隨機(jī)試驗(yàn)， rando.m ”將更加有用，這也是自己編制的程序，保存在：d:MATLAB6p5workshow.m 。62： d=1 1 1 1 1 1/6擲色子的結(jié)果概率是一個(gè)行向量(或者1X6矩陣)。63 ：sum(d)確認(rèn)它們的和為1!64： rando

14、(d)用這些概率去模擬擲色子的每個(gè)結(jié)果。用向上的光標(biāo)鍵重復(fù)這個(gè)命令幾次。模擬擲色子的另一個(gè)簡(jiǎn)單的方法是放大均勻分布隨機(jī)數(shù)后取整，floor(1+6*rand(1)。: vrando(1 1 1 1/4,20)程序rando的向量版本。每個(gè)數(shù)是等概率出現(xiàn)的。: show(vrando(1 1 1 1/4,100), BGSU )隨機(jī)地生成字符 B、G 臣口 U。出現(xiàn) BUGS之前，BGS曲現(xiàn)了嗎？七、寫一個(gè) Matlab程序你將創(chuàng)建一個(gè)新的 Matlab程序，名字為mywalk.m,用它來模擬100步的隨機(jī)游動(dòng)。在“file 菜單下有一個(gè)空白的按鈕，按下它即打開一個(gè)新的編輯窗口。在那個(gè)窗口里，分

15、行輸入下面的命令，然后保存該程序?yàn)閙ywalk.m。如果你保存在新的文件夾里，請(qǐng)確認(rèn)這個(gè)文件夾是否已加入到Path中或者改變?yōu)镃urrent Directory 。：n = 100;選取步數(shù)。:x = rand(n,1);生成均勻分布隨機(jī)數(shù)。: y = 2*(x 0.5) - 1;轉(zhuǎn)換這些數(shù)到為-1和+1。：z = cumsum(y);計(jì)算y的累積和。71 ： clf72: plot(z)畫出z的第1, 2, 3,.等的值。在command window口中輸mywalk,運(yùn)行(按回車)該程序，然后用光標(biāo)鍵多次重復(fù)它。如果有錯(cuò)誤提示，檢查你的輸入是否是正確的。73:運(yùn)行幾次后，你或許想一次就生

16、成一個(gè)更長(zhǎng)的字符串。到此目的的一個(gè)好的方法是將mywalk.m改為帶參數(shù)的Matlab function ,這樣就可以調(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)?；氐絚ommand window口輸入：75 ： mywalk(1000);運(yùn)行參數(shù)為 1000的程序 mywalk。76 : M=rand(6,6)77:M(2,:)78:

17、M(:,4)79:diag(M)80：sum(M)81 ： sum(M )八、矩陣6x 6的隨機(jī)數(shù)矩陣。取出矩陣M勺第2行。取出矩陣M勺第4歹U。取出矢I陣M勺對(duì)角線元素。矩陣列求和。對(duì)矩陣M勺行求和?！北硎巨D(zhuǎn)置。九、Markov 鏈在第66行中，序列中字母的出現(xiàn)是相互獨(dú)立的。我們將建立下面的一種情形，B通常跟隨在U之后，但決不跟在 G之后。出現(xiàn)B后，依概率向量0.2 0.6 0.2 0選擇下一個(gè)字母。GH現(xiàn)后，又以另一不同的概率向量出現(xiàn)下一個(gè)字母，以此類推。為此，我們將創(chuàng)建名字為 BGSU_markov.m勺新程序。打開一個(gè)新的編輯窗口，輸入下面的命令，然后再命令窗口輸入BGSU_marko

18、運(yùn)行之。82 : P=0.2 0.6 0.2 0; 0 0.2 0.6 0.2; 0.2 0 0.2 0.6; 0.6 0.2 0 0.2; P是一個(gè)4X4矩陣。每一行表明將以多大的概率選擇下一個(gè)字母。第一行即是數(shù)字1之后(對(duì)應(yīng)字母B)的概率，第二行是數(shù)字 2之后(G的概率等等。x(1) = rando(1 1 1 1/4);隨機(jī)地選擇第一個(gè)狀態(tài)。fo門=1:399,x(i+1) = rando(P(x(i),:);這是非常明智的：無論在哪個(gè)時(shí)刻，x(i)的值是多少，P(x(i),:)總是矩陣P的第x(i)行。該行的概率作為rando的參數(shù)來生成下一個(gè)狀態(tài)。endshow(x, BGSU );

19、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)=FA-1(y),如果u(1)., u(n)是服從(0

20、, 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 xh,令t為事件發(fā)生前的等待時(shí)間。t=-log(rand)/lambda; %服從參數(shù)為lambda的指數(shù)分布Exp(lambda)的隨機(jī)數(shù)。t=-log(rand(1,m)./lambda; % 服從 Exp(lambda)的 ntt行向量。2、Pareto 分布 simpareto.m概率密度函數(shù)：f(x)=alpha/(

21、1+xF(1+alpha), x0累積分布函數(shù)：F(x) = 1-(1+x)A(-alpha)。這是帶有所謂重尾分布中最簡(jiǎn)單的分布列子。產(chǎn)生一個(gè)均勻分布的樣本，并用分布函數(shù)的 反函數(shù)：sample = (1-rand(1, npoints).A(-1/alpha)-1;3、標(biāo)準(zhǔn) Pareto 分布 simparetonrm.m概率密度函數(shù)：f(x)=gamma*alpha/(1+gamma*x)A(1+alpha), x0累積分布函數(shù)：F(x) = 1-(1+gamma*x)A(-alpha)其中，參數(shù)gammai用來控制期望值的。在分布有重尾的情況下，若 1alphacumprob(j) &

22、 (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)的唯！隨機(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 是取

23、小于它的最小整數(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ù)的nxnm巨陣。如果沒有上面的命令，也可以用如下的命令替代之。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。r

24、andn(1,m);%服從標(biāo)準(zhǔn)正態(tài)分布 N(0,1)的m維隨機(jī)數(shù)行向量。randn(n,m);%每個(gè)分量是服從N(0,1)的nxm巨陣。mu+sigma.*randn(1,m); % m 個(gè)服從 N(mu,sigmaA2)分布的隨機(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

25、在S中。每個(gè)S(j) , j=1, , ,n,以1/2的概率獨(dú)立選擇0或1。for j=1:ns(j)=int(rand(1)+1);j=j+1；end3、隨機(jī)排列假如我們向隨機(jī)地排列 a(1), , ,a(n), 一個(gè)快速的方法是每一次互換兩個(gè)數(shù)的位置，共 n-1 次。for j=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、在一些模才程序中(比如更新過程),把概率分布作為一個(gè)外部參數(shù)來傳遞是很方便的。這是通過創(chuàng)建一個(gè) MATLA函數(shù)來實(shí)現(xiàn)。例如， rand, simpareto。分布的參

26、數(shù)是作為數(shù)組來傳遞的 (在rand中為空數(shù)組 , simpareto中為參量alpha1.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_0Ax (1-F(y)dy。例如，密度函數(shù)為 f(x)=2-2*x, 0=x=1 的線性分布是(0, 1)上均勻分布上的平穩(wěn)分布。參數(shù)為(alpha-1) Pareto分布是參數(shù)為alpha

27、的Pareto分布，參數(shù)為lambda的指數(shù)分布的平穩(wěn)分布就是自己。這將出現(xiàn)在平穩(wěn)更新過程或者排隊(duì)系統(tǒng)的平穩(wěn)版本的 例子中。distrstat.m 是一個(gè)表-查找函數(shù)，給定一個(gè)外部隨機(jī)數(shù)生成器，返回它的平穩(wěn)分布隨機(jī)數(shù)生成器。兩個(gè)參數(shù)都以數(shù)組的形式給出。至于應(yīng)用，可參見平穩(wěn)更新計(jì)數(shù)過程。例如：statdist, statpar = distrstat(rand, );第三章基本隨機(jī)過程兩個(gè)基本機(jī)制是離散時(shí)間的隨機(jī)游動(dòng)與連續(xù)時(shí)間的泊松過程。這些過程是基于獨(dú)立的簡(jiǎn)單 模擬算法的原型。擴(kuò)展到二維和三維中的模擬是直接的。在這個(gè)部分討論的 Matlab 文件有：ranwalk.m, brownian.m,

28、 poissonti.m, poissonjp.m, ranwalk2d.m, ranwalk3d.m, bm3plot.m, poisson2d.m, poisson3d.m一、一維情形1、隨機(jī)游動(dòng).簡(jiǎn)單隨機(jī)游動(dòng)ranwalk.m“從0開始，向前跳一步的概率為p,向后跳一步的概率為1 p”p=0.5;y=0 cumsum(2.*(rand(1,n-1)0)arrt = 0; sort(rand(npoints, 1)*tmax);elsearrt = 0;end%i出計(jì)數(shù)過程stairs(arrt, 0:npoints);.固定時(shí)間區(qū)間，NHi程poissonti.m它被復(fù)雜化為前面算法的向

29、量形式。到達(dá)時(shí)間間隔為指數(shù)分布的更新過程也將使用相同的 算法。%各0賦給到達(dá)時(shí)間。tarr = zeros(1, nproc);%各指數(shù)分布的時(shí)間間隔求和作為矩陣的列。i = 1;while (min(tarr(i,:)=tmax)tarr = tarr; tarr(i, :)-10g(rand(1, nproc)/lambda;= i+1;end%i出計(jì)數(shù)過程stairs(tarr, 0:size(tarr, 1)-1);二、高維情形1、二維隨機(jī)游動(dòng)ranwalk2d.m在(u, v)坐標(biāo)平面上畫出點(diǎn)(u(k),v(k),k=1:n,其中(u(k)和(v(k)是一維隨機(jī)游動(dòng)。例子程序是用四種

30、不同顏色畫了同一隨機(jī)游動(dòng)的四條軌道。n=100000;colorstr=b r g y;for k=1:4z=2.*(rand(2,n)0.5)-1;x=zeros(1,2); cumsum(z);col=colorstr(k);plot(x(:,1),x(:,2),col);hold onendgrid 2、三維隨機(jī)游動(dòng) ranwalk3d.m三維空間和上面的一樣。10p=0.5;n=10000;colorstr=b r g y;for k=1:4z=2.*(rand(3,n)=1,k=2。加上初始條件n(1)=0 ,遞歸定義一個(gè) Markov鏈n(k) , k=1.試試參數(shù)p的不同取值。從

31、 長(zhǎng)遠(yuǎn)看，會(huì)發(fā)生什么呢 ？m=200;p=0.2;N=zeros(1,m); %初始化緩沖區(qū)A=geornd(1-p,1,m); %生成到達(dá)序列模型，比如，幾何分布for n=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的泊松過程決定。服務(wù)員為每位顧客的服務(wù)時(shí)間服從指數(shù)分布，均值是 1/mu。由此得到的系統(tǒng)規(guī)模 N(t), t=0,是一個(gè)連續(xù)時(shí)間的Markov過程，其演變?nèi)缦隆腘(0)=n_0開始。等待強(qiáng)度為lambda+

32、mU勺指數(shù)分布 時(shí)間(如果n_0=0,強(qiáng)度為lambda),然后以可能性lambda/(lambda+mu)向前跳躍和以可能性 mu/(lambda+mu)向后跳躍。如此循環(huán)。在N=0時(shí)動(dòng)力學(xué)的改變是在開始用一個(gè)短循環(huán)來實(shí)現(xiàn)的：if i=0mutemp=0; elsemutemp=mu; end主循環(huán)僅僅檢查向前跳或者向后跳：if rand=0, 在t_k時(shí)向前跳躍，在u_k時(shí)向后跳躍。假設(shè)我們有長(zhǎng)度為n到達(dá)時(shí)向量t_k。這里是得到離開時(shí)間的一種方法：arrsubtr=arrtime-(0:n-1); % t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptim

33、e=(1:n)+max(triu(arrmatrix)現(xiàn)在想畫出N(t):B=ones(n,1) arrtime ; -ones(n,1) deptime;Bsort=sortrows(B,2);%按次序?qū)⑻S分類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/muo 一個(gè)相似的計(jì)算離開時(shí)間的技術(shù)在此情形也是可以的。特別地，取Exp(mu)分布的服務(wù)時(shí)間就是 M/M/1系統(tǒng)。因而，我們有了模擬M/M/1的另一種方法。如果lambda/m

34、u0)statdist,statpar=distrstat(servdist,servpar);rndpar1=nstart,1,statpar:; %stattimes=feval(statdist, rndpar1:);%arrtimes=zeros(size(stattimes); %泊松隨機(jī)變量%平穩(wěn)分布句柄隨機(jī)數(shù)生成器參數(shù)平穩(wěn)服務(wù)時(shí)間在t = 0時(shí)刻前到達(dá)的顧客按到達(dá)時(shí)間為 0來看待。9、end 10、一旦創(chuàng)建了計(jì)數(shù)過程，就刪除開始時(shí)負(fù)到達(dá)額外的0點(diǎn)。增加maxtime到跳躍，以得到正確 的圖形。七、生滅過程1、一般的生滅強(qiáng)度 birthdeath.m作為例子，我們選擇在水平i上出生

35、強(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);if i=0mu_i=0;elsemu_i=mu*i;endtime=-log(rand)./(lambda_i+mu_i);2、Moran模型 moran.m另一個(gè)起源于遺傳學(xué)的生滅過程。有限狀態(tài)空間，吸收界限。x=1:N+1;lambda(x)=(x-1).*(1-(x-1)./N); %出生率。14mu(x)=(x-1).*(1-(x-1)./N);%死亡率。q(x)=

36、lambda(x)+mu(x); %兩次跳躍時(shí)間間隔的指數(shù)分布率。八、分支過程Galton-Watson 過程 galtonwatson.m離散時(shí)間，生命長(zhǎng)度為1。死亡的每個(gè)個(gè)體產(chǎn)生隨機(jī)的后代個(gè)數(shù)。函數(shù)offspring(k) 給出從人口規(guī)模為k開始的祖先向量。p=1/2 0 1/2;z=cumsum(p);n=length(p); %可能的子孫數(shù)量offmu=dot(0:n-1,p); %子孫的平均個(gè)數(shù)u1=sort(rand(1,k);for j=1:nu(j)=length(find(u1 z(j);endu=diff(0 u);nu=u*(0:n-1);九、計(jì)數(shù)過程計(jì)數(shù)過程N(yùn)(t)記錄

37、的是實(shí)值隨機(jī)點(diǎn)過程T_k在區(qū)間0, t) 內(nèi)的點(diǎn)數(shù)。泊松過程及排隊(duì)系 統(tǒng)中遇到的系統(tǒng)規(guī)模過程都是計(jì)數(shù)過程。十、更新過程更新過程，是一列獨(dú)立同分布的正隨機(jī)變量的部分和序列。這個(gè)過程可以被想象為：當(dāng)同 種生物的一些個(gè)體生命期結(jié)束，同時(shí)他們也被新的生命所代替的時(shí)間點(diǎn)序列。更新計(jì)數(shù)過程記 錄的是在時(shí)間區(qū)間0, t)內(nèi)更新的次數(shù)。它是一個(gè)隨機(jī)階梯函數(shù)。15第五章模擬的應(yīng)用和例子大數(shù)定律表明：1、經(jīng)驗(yàn)均值收斂到它的期望值；2、統(tǒng)計(jì)物理中，軌道平均與總體平均是漸進(jìn)相同的；3、為隨機(jī)模擬提供了理論基礎(chǔ)，并建立了事件頻率和概率的聯(lián)系。一、計(jì)算積分 I=int_aAb f(x)dx1、I=(b-a)E(f(a+(

38、b-a)U), Usim (0,1)。模擬 X_j=E(f(a+(b-a)U_j)，用平均1/nsum_j=1An X_j逼近I/(b-a)。大數(shù)定律說明，可以用獨(dú)立試驗(yàn)的頻率來近似期望值。2、選取一個(gè)包含函數(shù)f圖形的矩形，比如a,bminf,maxf。生成該矩形上n對(duì)均勻分布隨機(jī)數(shù)(X,Y),記錄事件“Yt)t)frac1n tA2Var(X).然而，Var(X)往往也是不知道的，這是由 1/nsum_j=1An(X_j-barX)A2來近似的，其中barX=1/nsum_j=1An X_j 。2、中心極限定理。由中心極限定理知道，近似地有frac1/nsum_j=1AnX_j-EXVar(

39、X)/n服從標(biāo)準(zhǔn)正態(tài)。因此， P(|1/nsum_j=1AnX_j-EX|t)=2phi(fractVar(X)/n)-1=q。即以 q 的概率誤差小于fractVar(X)/n 。當(dāng)然，Var(X)也如上面一樣來近似。中心極限定理的另一個(gè)用處是模擬某種連續(xù)過程的軌道，即X_n(t)=frac1sqrtnsum_k=ntxi_k,其中 xi_k 為獨(dú)立的 0均值隨機(jī)變量。當(dāng)然，它還可以用來近似模擬標(biāo)準(zhǔn)正態(tài)隨機(jī)數(shù)。三、減小方差技術(shù)(對(duì)立變量、條件期望、控制變量、重要樣本)上面提到，方差往往是不知道，它也是通過產(chǎn)生的隨機(jī)數(shù)來估計(jì)。由估計(jì)的誤差分析知，該方差越小越好，因此給出幾種減小方差的一般方法。

40、16四、模擬的例子（一）概率問題的模擬問題一車和羊的游戲；問題二蒲豐投針問題；問題三擲骰子問題； 問題四無記憶性的 例子；問題五 生日問題；問題六 Galton 釘板實(shí)驗(yàn)；問題七 趕火車問題；問題一車和羊的游戲假設(shè)你在進(jìn)行一個(gè)游戲節(jié)目。現(xiàn)給三扇門供你選擇：一扇門后面是一輛轎車，另兩扇門后 面分別都是一頭山羊。你的目的當(dāng)然是要想得到比較值錢的轎車，但你卻并不能看到門后面的 真實(shí)情況。主持人先讓你作第一次選擇。在你選擇了一扇門后，知道其余兩扇門后面是什么的 主持人，打開了另一扇門給你看，而且，當(dāng)然，那里有一頭山羊?，F(xiàn)在主持人告訴你，你還有 一次選擇的機(jī)會(huì)。那么， 請(qǐng)你考慮一下，你是堅(jiān)持第一次的選擇

41、不變，還是改變第一次的選擇，更有可能得到轎車？廣場(chǎng)雜志刊登出這個(gè)題目后，竟引起全美大學(xué)生的舉國(guó)辯論，許多大學(xué)的教授們也參 與了進(jìn)來。真可謂盛況空前。據(jù)紐約時(shí)報(bào)報(bào)道，這個(gè)問題也在中央情報(bào)局的辦公室內(nèi)和波 斯灣飛機(jī)駕駛員的營(yíng)房里引起了爭(zhēng)論，它還被麻省理工學(xué)院的數(shù)學(xué)家們和新墨哥州洛斯阿拉莫 斯實(shí)驗(yàn)室的計(jì)算機(jī)程序員們進(jìn)行過分析。問題分析在一次實(shí)驗(yàn)中，如果第一次選擇選中了轎車（概率為1/3）,那么主持人打開一扇門后，如果堅(jiān)持原來的選擇，則能得到轎車，反之，改變第一次選擇則不能得到轎車。如果第一次沒有 選中轎車（概率為 2/3）,那么其余兩扇門后面必有一個(gè)是轎車，主持人只能打開有山羊有那扇 門，則剩下的一

42、扇門后面是轎車，此時(shí)堅(jiān)持原來的選擇不能得到轎車，改變第一次的選擇必能 得到轎車。因此，經(jīng)過分析，堅(jiān)持第一次的選擇不變得到轎車的概率為1/3,改變第一次的選擇得到轎車的概率為 2/3。實(shí)際上，在只有三扇門的情況下，那么改不改變選擇效果并不明顯。如果有100扇門，參與的嘉賓選擇了其中的一扇， 而主持人隨后把剩下的 99扇門中間的98扇門都打開，這98扇門 后面都沒有獎(jiǎng)品，這時(shí)應(yīng)該改變選擇，畢竟最開始自己選擇的那扇門中獎(jiǎng)的概率只是1%而已。需要注意的是，主持人是在知道其他兩扇門后面都有什么的情況下選擇一個(gè)門打開的。這 種情況下三個(gè)門后是轎車的概率因?yàn)橹鞒秩酥澜Y(jié)果并參與其中而關(guān)聯(lián)在一起，而不是孤立等

43、 同的。如果打開門的不是主持人，而是另一個(gè)參與者，并且當(dāng)他打開門時(shí)發(fā)現(xiàn)什么也沒有，那 么，剩下的兩個(gè)門后是轎車的概率才是相等的。計(jì)算機(jī)模擬為了驗(yàn)證這一結(jié)果，我們就要比較不改變選擇中獎(jiǎng)的幾率和改變選擇中獎(jiǎng)的幾率。模擬方法是：我們從0, 1, 2這3個(gè)數(shù)中隨機(jī)一個(gè)為轎車（即中獎(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扇門，參與的嘉賓選擇了其中的一扇，而主持人隨后把剩下的99

44、扇門中間的98扇門都打開，這98扇門后面都沒有獎(jiǎng)品，然后模擬并比較不改 變選擇中獎(jiǎng)的幾率和改變選擇中獎(jiǎng)的幾率。此時(shí)的情況也是相同的， 只是每次隨即都是從0到99中隨機(jī)數(shù)而已。17結(jié)果及分析下面兩幅圖分別是3個(gè)門時(shí)不改變選擇中獎(jiǎng)的概率在 畋模擬結(jié)果下的概率分布（第二幅是 為了便于觀察特意畫在固定坐標(biāo)軸上的）。卜面則是100個(gè)門的情況下，不改變選擇中獎(jiǎng)的概率分布:可以顯然地看出在主持人幫助的情況下，改變選擇是可以大大增加自己中獎(jiǎng)的幾率的。通過這樣一個(gè)例子，我并不是想說明什么概率意義上的問題。只是通過這么一個(gè)模擬過程 來學(xué)習(xí)計(jì)算機(jī)隨機(jī)模擬的一些基本方法與技巧。像在主持人不知道內(nèi)幕隨機(jī)的打開一個(gè)發(fā)現(xiàn)是

45、 山羊這，我們可以通過同樣的隨機(jī)模擬過程來模擬這種情況。并可以驗(yàn)證改變選擇與否對(duì)自己 中獎(jiǎng)的影響是相同的。當(dāng)模擬的次數(shù)逐漸的增多時(shí)，其模擬值越接近理論值，這說明模擬的效果越好。在隨機(jī)事 件的大量重復(fù)出現(xiàn)中，往往呈現(xiàn)幾乎必然的規(guī)律，這個(gè)規(guī)律就是大數(shù)定律。通俗地說，這個(gè)定 理就是，在試驗(yàn)不變的條件下，重復(fù)試驗(yàn)多次，隨機(jī)事件的頻率近似于它的概率。偶然必然中 包含著必然。此次模擬試驗(yàn)也正好用實(shí)際的模擬例子說明了大數(shù)定理的正確性和應(yīng)用性。Matlab程序1、編寫函數(shù)n=10000; %實(shí)驗(yàn)次數(shù)stick=0; %堅(jiān)持選擇的獲獎(jiǎng)次數(shù)18alter=0; %改變選擇的獲獎(jiǎng)次數(shù)for i=1:nx=floor

46、(3*rand(1)+1);%y=floor(3*rand(1)+1) ;%if(x=y) %stick=stick+1;else %alter=alter+1;endendstick=stick/n %在第x扇門后有轎車選擇第y扇門第一次選中的情況第一次沒選中的情況轉(zhuǎn)化為百分比alter=alter/n保存為 myone.m, 在 command window輸入 myone,結(jié)果為 stick =0.3349alter =0.66512、moni.m這是模擬不同次數(shù)并繪圖的程序function moni()N=10:100:10000;for i=1:length(N)bubian(i)=

47、bubianmoni(N(i);endgrid onaxis(0,10000,0,1)hold on(這3行為固定坐標(biāo)軸和柵格的)plot(N,bubian)3、bubianmoni.m這是模擬并返回不改變中獎(jiǎng)的概率的程序function x=bubianmoni(n)A=0;B=0;for i=1:n19a=randint(1,1,0,99);b=randint(1,1,0,99);if a=bA=A+1; %不變的累加elseB=B+1; %改變的累加endendx=A/n; %不變中獎(jiǎng)的概率 其中對(duì)改變選擇的累計(jì)在此處是可以加以省略的。問題二蒲豐投針問題1777年法國(guó)科學(xué)家蒲豐提出的一種

48、計(jì)算圓周率的方法-隨機(jī)投針法，這種方法的具體步驟是：1、取一張白紙，在上面畫上許多條間距為d的平行直線。2、取一根長(zhǎng)為l(ld)的針，隨機(jī)向畫有平行直線的紙上擲n次，觀察針與直線相交的次數(shù)，記為 m 3、計(jì)算針與直線相交的概率。問題分析針與平行線相交的充要條件是：x Wsin巾，其中2d .0MxM ,0 2建立直角坐標(biāo)系(6, x),上述條件在坐標(biāo)系下將是曲線所圍成的曲邊梯形區(qū)域:20由幾何概率知:1 二.,2ln g的面積 2 0 sin d 2l 一 TT：；dmG的面積d ndji2經(jīng)統(tǒng)計(jì)實(shí)驗(yàn)估計(jì)出概率 p憶m,有上式即 m= NL ,故冗二迎。 nn 二 d dmMatlab程序n=

49、lnput（ 請(qǐng)輸入投針的次數(shù)n:）l=1;d=2;m=0;for k=1:n%模才n n 次x=unifrnd(0,d/2); %x為0到d/2之間的隨機(jī)變量，服從均勻分布y=unifrnd(0,pi); %y為0到pi之間的隨機(jī)變量，服從均勻分布if x0.5*l*sin(y);m=m+1;elseendendp=m/n;pi_m=1/p;fprintf（所模擬出的圓周率為：%fn,pi_m）運(yùn)行結(jié)果與分析運(yùn)行了 4次程序，并賦予不同的n值，得如下結(jié)果：n =100所模擬出的圓周率為：3.225806n =1000所模擬出的圓周率為：3.257329n =10000所模擬出的圓周率為：3.

50、176620n =100000所模擬出的圓周率為：3.136664可見，模擬的次數(shù)越多，所得的結(jié)果就越準(zhǔn)確。這與理論所知的結(jié)果是一樣的。問題三擲骰子問題21 擲三只相同的色子，求有兩個(gè)點(diǎn)數(shù)相同的概率。問題分析構(gòu)造3個(gè)整數(shù)隨機(jī)變量，使之服從1到6之間的均勻分布，比較每次模擬出的 3個(gè)值, 如果至少有兩個(gè)相同，則將成功的次數(shù)加1,分析可知，成功的頻率即為所求的概率。Matlab程序n=input(,請(qǐng)輸入模擬的次數(shù):)1=0;for i=1:nx=unidrnd(6,1,3); %x 中含有3個(gè)元素，相互獨(dú)立，且服從 0到6的整數(shù)均勻分布if x(1)=x(2)|x(2)=x(3)|x(1)=x(

51、3)l=l+1;%比較模擬出的結(jié)果，1為成功數(shù)end;end;p=1/n;%計(jì)算成功的頻率，即為概率fprintf(至少兩個(gè)相同的概率是%fn,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。這與以上模擬出的值基本上是一樣的。問題四 無記憶性的例子考慮有兩個(gè)營(yíng)業(yè)員白郵局。假設(shè)當(dāng) Smith先生進(jìn)入郵局的時(shí)候，

52、他看到兩個(gè)營(yíng)業(yè)員分別在 為Jones先生和Brown先生服務(wù)，并且被告知先服務(wù)完的營(yíng)業(yè)員即將為他服務(wù)。再假設(shè)為每顧 客服務(wù)的時(shí)間都服從參數(shù)為.的指數(shù)分布。.求這三個(gè)顧客中 Smith最后離開郵局的概率？ Jones和Brown呢？.求Smith在郵局花費(fèi)的平均時(shí)間 ET ?問題分析當(dāng)Smith接受服務(wù)的時(shí)候，Jones和Brown只有一個(gè)人離開，另一人仍在接受服務(wù)。 然而,22由指數(shù)分布的無記憶性，從現(xiàn)在起剩下的這個(gè)人被服務(wù)的時(shí)間是參數(shù)為.的指數(shù)分布，即與Smith是同分布的。再由對(duì)稱性，他與 Smith先離開的概率一樣，都為 1/2。第1問答案為 1/2,1/4,1/4。ET = Emin

53、X;Y+ S=1/(2* 入)+1/ 入,其中 S 為 Smith 服務(wù)花費(fèi)的時(shí)間,X;Y 分別為Jones和Brown服務(wù)所花費(fèi)的時(shí)間。Matlab程序編寫函數(shù) function mytwo(lambda) n=100000; % 實(shí)驗(yàn)次數(shù) jones=0; % brown=0; smith=0;sumT=0; % for i=1:nt1=-log(rand)/lambda; %t1,t2,t3 t2=-log(rand)/lambda;t3=-log(rand)/lambda;if(t1+t3t2) %brownbrown=brown+1;t=t1+t3;elseif(t2+t3t1) %

54、jones jones=jones+1; t=t2+t3; else%smithsmith=smith+1; t=t3+(t1=t2)*t2;end sumT=sumT+t;end smith=smith/n % jones=jones/n brown=brown/n ET=sumT/n %smith三個(gè)變量記錄各人最后離開次數(shù)統(tǒng)計(jì)n次實(shí)驗(yàn)smith花費(fèi)時(shí)間總和分別為三人服務(wù)時(shí)間最后離開最后離開最后離開%t=t3+min(t1,t2)轉(zhuǎn)化為概率平均花費(fèi)時(shí)間ETreal=3/(2*lambda) %分析出的ET,用來與模擬結(jié)果對(duì)比保存為 mytwo.m,在 command window輸入 my

55、two(3),結(jié)果為 smith =0.5019 jones = 0.2495 brown =0.248623ET =0.5009ETreal =0.5000問題五生日問題有n個(gè)人的班級(jí)中至少有兩個(gè)人生日是同一天的概率，模擬之。問題分析假設(shè)：所有人生日是相互獨(dú)立的。設(shè)一年有y天，為求至少兩個(gè)人生日在同一天的概率，我們可以先求n個(gè)人生日全不同的概率：nn個(gè)人生日所有可能的組合(不考慮次序)種數(shù)有：工n!而n個(gè)人生日都不一樣的可能的組合方式種數(shù)有：cy = y!(y -n)!* n!故任意兩人生日不同的理論概論是：(這樣寫可以方便編程計(jì)算)y! (y-n)!*n! y! 1*2*3*.*( y -

56、 n)*( y - n 1)*.* yP0 二=7 二乙y2(y-n)!* yn1*2*3*.*( yn)*y*y*.*yn!從而有至少兩人生日同一天的概率是p = 1 p0注：顯然當(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)。通過調(diào)用子函數(shù)計(jì)算隨機(jī)模擬的概率以及理論概率，將二者畫在 同一個(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中。使用

57、unique指令，將a的相同項(xiàng)合并，并按照從小到大順序排列。得到的新數(shù) 組與原來的數(shù)組 a進(jìn)行長(zhǎng)度的比較，如果新數(shù)組的長(zhǎng)度小于原數(shù)組，說明有兩人或者以上同一 天生日，此時(shí)成功次數(shù)加1?？偣步?jīng)過tn次隨機(jī)試驗(yàn)，即可求出至少兩個(gè)人同一天生日的試驗(yàn)的頻率。子程序2： func2.m計(jì)算理論概率通過上面的分析 p0表達(dá)式分子是一個(gè)1*y的向量的連乘，分母也是一個(gè)1*y維向量的連乘。對(duì)應(yīng)分量做除法后求新向量個(gè)元素乘積即可。P就可以通過p=1-p0來計(jì)算了。結(jié)果分析24987654321 O0.口口O.口O.Q程序運(yùn)行所花費(fèi)時(shí)間為4.380000e-001秒_ - _ i _ 一 _ - -C 9 0 7

58、 5 5 4 3 2 1 O O.O.O.Q cio.iDio.a.符蜜嵌糕格里程序運(yùn)行所花費(fèi)時(shí)間為8.750000e-001秒9B7B54-32 1 OD.D.0.口D.0.口0.0.科壹牘腎部曼圖三程序運(yùn)行所花費(fèi)時(shí)間為6.688000e+000秒25舟棗嘉M糜圖四程序運(yùn)行所花費(fèi)時(shí)間為2.438000e+000秒D W 2D 300505071080 9D WD學(xué)生人數(shù)n圖五程序運(yùn)行所花費(fèi)時(shí)間為2.232800e+001秒9 97654321 O D.D.o.ciD.o.o.QD.圖六程序運(yùn)行所花費(fèi)時(shí)間為2.219210e+002秒26tn(1)比較上面幾幅圖，對(duì)應(yīng)的學(xué)生人數(shù)均為50人。而隨

59、機(jī)模擬的試驗(yàn)次數(shù)在增加，當(dāng)比較小的時(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è)人同一天生日的概率也逐漸增大，而且通過概率線可以看出。當(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%

60、A上,這是幾乎可以肯定的說至少有兩個(gè)人是同一天生日的。(3)注：由于數(shù)字是隨機(jī)生成的，所以每次運(yùn)行的結(jié)果都不一樣，但是大體上的趨勢(shì)是如此的。另外，這個(gè)運(yùn)行時(shí)間也僅供對(duì)比參考，他還與及其運(yùn)行的狀態(tài)有關(guān)系，因此沒有太大的實(shí) 際意義。Matlab程序程序說明：主程序名suiji.m子程序1 ： func1.m子程序2 ： func2.m程序如下：suiji.m 陶延淵%2008年5月%每年的天數(shù)year%A數(shù) n number%隨機(jī)模才次數(shù)tn test number clear;clc;c1=clock;n=50; %學(xué)生人數(shù)tn=100; %隨機(jī)模擬次數(shù) y=365; %一年的天數(shù)，一般計(jì)算取3