隨機過程試驗講義

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

2、YPERLINK l bookmark93 o Current Document 第五章模擬的應(yīng)用和例子16附錄各章的原程序51參考文獻75若想檢驗數(shù)學模型是否反映客觀現(xiàn)實，最自然的方法是比較由模型計算的理論概率和由客 觀試驗得到的經(jīng)驗頻率。不幸的是，這兩件事都往往是費時的、昂貴的、困難的，甚至是不可 能的。此時，計算機模擬在這兩方面都可以派上用場：提供理論概率的數(shù)值估計與接近現(xiàn)實試 驗的模擬。模擬的第一步自然是在計算機程序的算法中如何產(chǎn)生隨機性。程序語言，甚至計算器，都 提供了 “隨機”生成0,1區(qū)間內(nèi)連續(xù)數(shù)的方法。因為每次運行程序常常生成相同的“隨機數(shù)”因此這些數(shù)被稱為偽隨機數(shù)。盡管如此，

3、對于多數(shù)的具體問題這樣的隨機數(shù)已經(jīng)夠用。我們將 假定計算機已經(jīng)能夠生成0,1上的均勻隨機數(shù)。也假定這些數(shù)是獨立同分布的，盡管它們常 常是周期的、相關(guān)的、，。,本講義的安排如下，第一章是Matlab簡介，從實踐動手角度了解并熟悉 Matlab環(huán)境、命令、 幫助等，這將方便于 Matlab的初學者。第二章是簡單隨機變量的模擬，只給出了常用的 Matlab 模擬語句，沒有堆砌同一種變量的多種模擬方法。對于沒有列舉的隨機變量的模擬，以及有特 殊需求的讀者應(yīng)該由這些方法得到啟發(fā)，或者參考更詳細的其他文獻資料。第三章是基本隨機 過程的模擬。主要是簡單獨立增量過程的模擬，多維的推廣是直接的。第四章是Mark

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

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

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

7、集合， 它是一個列向量或者是 10 x 1的矩陣。10：x+511 : 1000*x的每個分量都加5。 的每個分量都乘以1000。12： x=rand(10,7) 10 X 7的隨機數(shù)矩陣。 若想重復(fù)此命令或其他命令，按住向上的光標鍵直至看到想重復(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ù)試試。用分號“；”可以

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

9、。對大于1000的數(shù)試試結(jié)果。二、獲取更多幫助32:如果你想查找不會使用的命令，可以點擊： Help/ MATLAB Help ,打開左邊的“ MATLAB 節(jié)，選擇Functions - Categorical List ”即可。據(jù)我所知，這是尋求幫助的最好方法。三、畫出數(shù)據(jù)點33: plot(x(1:10), *);用“*”描出x的前10個點。注意兩個單引號為英文的單引號，下同。34： plot(x-0.5);向下平移0.5 ,描出述據(jù)點，且將其連成線。: hold on將下面的圖形加到上面的圖形中。: plot(cumsum(x- 0.5), r);將這個結(jié)果圖畫到上面的圖形中。r表示用

10、紅色的線繪出，而缺省的顏色為藍色。： zoom on用鼠標點擊可放大圖形，雙擊回到原始的尺寸。38：clf清除當前的圖形。39:z=randn(1000,1);生成1000個標準正態(tài)分布隨機數(shù)。40 ：w=z+randn(1000,1);生成依賴z的隨機數(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;定義一個時間點向量，間隔為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)給圖形加上更恰當?shù)臉祟}。五、運行現(xiàn)有的Matlab程序:上網(wǎng)下載或者拷貝一些編輯好的Matlab程序到自己的電腦中。：如果在你電腦的某個文件夾中有現(xiàn)成的Matlab程序(*.m)，可以設(shè)置Current Directory (Command WindoWB 口的上面)為該文件夾即可運行這些程序。53:如果在你電腦里的幾個文件夾里都有Matlab程序，點擊菜單中：File/ Set Pat

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

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

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

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

16、成一個更長的字符串。到此目的的一個好的方法是將mywalk.m改為帶參數(shù)的Matlab function ,這樣就可以調(diào)用它。74：在你的程序中，將行“ n = 100; ”替換為function z = mywalk(n)這樣，mywalk是一個帶參數(shù)n的函數(shù)(生成序列的長度)，返回變量z。函數(shù)里面的變量(比如 y)是內(nèi)部變量，它的值不被帶到函數(shù)外面。就像 sin和rand一樣， 函數(shù)mywalk返回一個值(向量 z)。回到command window口輸入：75 ： mywalk(1000);運行參數(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的隨機數(shù)矩陣。取出矩陣M勺第2行。取出矩陣M勺第4歹U。取出矢I陣M勺對角線元素。矩陣列求和。對矩陣M勺行求和。”表示轉(zhuǎn)置。九、Markov 鏈在第66行中，序列中字母的出現(xiàn)是相互獨立的。我們將建立下面的一種情形，B通常跟隨在U之后，但決不跟在 G之后。出現(xiàn)B后，依概率向量0.2 0.6 0.2 0選擇下一個字母。GH現(xiàn)后，又以另一不同的概率向量出現(xiàn)下一個字母，以此類推。為此，我們將創(chuàng)建名字為 BGSU_markov.m勺新程序。打開一個新的編輯窗口，輸入下面的命令，然后再命令窗口輸入BGSU_marko

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

19、hist(x,4);第二章簡單分布的模擬Matlab里生成0 , 1上的均勻隨機數(shù)的語句是：rand(1,1); rand(n,m) 。一旦有了 0, 1上均勻隨機數(shù)，則我們就能夠做下面的事情。在這個部分討論的 Matlab 文件有：simexp.m, simpareto.m , simparetonrm.m, simdiscr.m, simbinom.m, simgeom.m, distrmu.m, distrstat.m 。一、一般連續(xù)分布(逆變換法、拒絕法、Hazard率方法)生成有連續(xù)分布函數(shù)隨機數(shù)的一般方法是用反函數(shù)法。設(shè)G(y)=FA-1(y),如果u(1)., u(n)是服從(0

20、, 1)上均勻分布的隨機數(shù)，那么 G(u(1), ., G(u(n)就是分布函數(shù)為 F(x)的隨機數(shù)。比如，指數(shù)分布， Pareto分布等。1、指數(shù)分布 simexp.m事件以強度lambda的時間隨機地發(fā)生，即事件在 t , t + h時間內(nèi)發(fā)生的可能性是lambda xh,令t為事件發(fā)生前的等待時間。t=-log(rand)/lambda; %服從參數(shù)為lambda的指數(shù)分布Exp(lambda)的隨機數(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)。這是帶有所謂重尾分布中最簡單的分布列子。產(chǎn)生一個均勻分布的樣本，并用分布函數(shù)的 反函數(shù)：sample = (1-rand(1, npoints).A(-1/alpha)-1;3、標準 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個以概率p為1,概率1-p為0的隨機數(shù)(m維行向量)。三、特殊分布1、二項分布 simbinom.m將每次成功的概率為p的試驗獨立做n次,設(shè)x是成功的個數(shù)x=sum(rand(n,m)=p);% x是服從 Bin(n,p)的唯！隨機數(shù)向量。2、幾何分布 simgeom.m實驗每次成功的概率為p,設(shè)x為第一次成功前失敗的次數(shù)。x=floor(-log(rand(1,m)./(-log(1-p); %服從參數(shù)為 p的幾何分布 Ge(p)的 m維隨機數(shù)行向量。floor 是取

23、小于它的最小整數(shù)的函數(shù)。3、泊松分布Matlab的統(tǒng)計工具箱含有產(chǎn)生泊松分布隨機數(shù)的命令，為 poissrnd。poissrnd(lambda);poissrnd(lambda, n, m); %產(chǎn)生參數(shù)為lambda的泊松分布 Po(lambda)隨機數(shù)的nxnm巨陣。如果沒有上面的命令，也可以用如下的命令替代之。arrival=cumsum(-log(rand(1,5)./lambda);n=length(find(arrival=lambda); %find是找出非 0值所在的位置，length 是它的維數(shù)。4、正態(tài)分布高斯分布，或正態(tài)分布的隨機數(shù)用Matlab生成的命令是randn。r

24、andn(1,m);%服從標準正態(tài)分布 N(0,1)的m維隨機數(shù)行向量。randn(n,m);%每個分量是服從N(0,1)的nxm巨陣。mu+sigma.*randn(1,m); % m 個服從 N(mu,sigmaA2)分布的隨機數(shù)四、離散試驗的模擬1、從1, , ,n中任取一個。int(n*rand(1,m)+1);從1, , ,n中任取不可重復(fù)兩個。a=int(n*rand(1)+1);b=int(n*rand(1)+1);while(a=b)b=int(n*rand(1)+1);end2、隨機子集模擬集合1, , ,n的隨機子集，我們是定義序列S(j)=0,1, S(j)=1即表示將j

25、在S中。每個S(j) , j=1, , ,n,以1/2的概率獨立選擇0或1。for j=1:ns(j)=int(rand(1)+1);j=j+1；end3、隨機排列假如我們向隨機地排列 a(1), , ,a(n), 一個快速的方法是每一次互換兩個數(shù)的位置，共 n-1 次。for j=n:2N=int(j*rand(1)+1);y=a(N);a(N)=a(j); a(j)=y; end五、外部參數(shù)的隨機數(shù)產(chǎn)生器distrmu.m1、在一些模才程序中(比如更新過程),把概率分布作為一個外部參數(shù)來傳遞是很方便的。這是通過創(chuàng)建一個 MATLA函數(shù)來實現(xiàn)。例如， rand, simpareto。分布的參

26、數(shù)是作為數(shù)組來傳遞的 (在rand中為空數(shù)組 , simpareto中為參量alpha1.4。distrmu.m是一個表-查找函數(shù)，從它的參數(shù)列表中提取期望參數(shù)的外部隨機數(shù)發(fā)生器，如:mu = distrmu(simparetonrm, 1.4, 2.5);2、平穩(wěn)分布 distrstat.m假設(shè)有一個分布函數(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)更新過程或者排隊系統(tǒng)的平穩(wěn)版本的 例子中。distrstat.m 是一個表-查找函數(shù)，給定一個外部隨機數(shù)生成器，返回它的平穩(wěn)分布隨機數(shù)生成器。兩個參數(shù)都以數(shù)組的形式給出。至于應(yīng)用，可參見平穩(wěn)更新計數(shù)過程。例如：statdist, statpar = distrstat(rand, );第三章基本隨機過程兩個基本機制是離散時間的隨機游動與連續(xù)時間的泊松過程。這些過程是基于獨立的簡單 模擬算法的原型。擴展到二維和三維中的模擬是直接的。在這個部分討論的 Matlab 文件有：ranwalk.m, brownian.m,

28、 poissonti.m, poissonjp.m, ranwalk2d.m, ranwalk3d.m, bm3plot.m, poisson2d.m, poisson3d.m一、一維情形1、隨機游動.簡單隨機游動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出計數(shù)過程stairs(arrt, 0:npoints);.固定時間區(qū)間，NHi程poissonti.m它被復(fù)雜化為前面算法的向

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

30、不同顏色畫了同一隨機游動的四條軌道。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、三維隨機游動 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 ,遞歸定義一個 Markov鏈n(k) , k=1.試試參數(shù)p的不同取值。從

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

32、mU勺指數(shù)分布 時間(如果n_0=0,強度為lambda),然后以可能性lambda/(lambda+mu)向前跳躍和以可能性 mu/(lambda+mu)向后跳躍。如此循環(huán)。在N=0時動力學的改變是在開始用一個短循環(huán)來實現(xiàn)的：if i=0mutemp=0; elsemutemp=mu; end主循環(huán)僅僅檢查向前跳或者向后跳：if rand=0, 在t_k時向前跳躍，在u_k時向后跳躍。假設(shè)我們有長度為n到達時向量t_k。這里是得到離開時間的一種方法：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ù)時間分布S,其均值為1/muo 一個相似的計算離開時間的技術(shù)在此情形也是可以的。特別地，取Exp(mu)分布的服務(wù)時間就是 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); %泊松隨機變量%平穩(wěn)分布句柄隨機數(shù)生成器參數(shù)平穩(wěn)服務(wù)時間在t = 0時刻前到達的顧客按到達時間為 0來看待。9、end 10、一旦創(chuàng)建了計數(shù)過程，就刪除開始時負到達額外的0點。增加maxtime到跳躍，以得到正確 的圖形。七、生滅過程1、一般的生滅強度 birthdeath.m作為例子，我們選擇在水平i上出生

35、強度為lambda_i=lambda/(1+i),死亡強度為mu_i=mu*i 的模型，其中l(wèi)ambda和mu固定的常數(shù)。要求循環(huán)滿足直到下次跳躍，跳躍強度和等待時間才被更新，即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另一個起源于遺傳學的生滅過程。有限狀態(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ù)分布率。八、分支過程Galton-Watson 過程 galtonwatson.m離散時間，生命長度為1。死亡的每個個體產(chǎn)生隨機的后代個數(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); %子孫的平均個數(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);九、計數(shù)過程計數(shù)過程N(t)記錄

37、的是實值隨機點過程T_k在區(qū)間0, t) 內(nèi)的點數(shù)。泊松過程及排隊系 統(tǒng)中遇到的系統(tǒng)規(guī)模過程都是計數(shù)過程。十、更新過程更新過程，是一列獨立同分布的正隨機變量的部分和序列。這個過程可以被想象為：當同 種生物的一些個體生命期結(jié)束，同時他們也被新的生命所代替的時間點序列。更新計數(shù)過程記 錄的是在時間區(qū)間0, t)內(nèi)更新的次數(shù)。它是一個隨機階梯函數(shù)。15第五章模擬的應(yīng)用和例子大數(shù)定律表明：1、經(jīng)驗均值收斂到它的期望值；2、統(tǒng)計物理中，軌道平均與總體平均是漸進相同的；3、為隨機模擬提供了理論基礎(chǔ)，并建立了事件頻率和概率的聯(lián)系。一、計算積分 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ù)定律說明，可以用獨立試驗的頻率來近似期望值。2、選取一個包含函數(shù)f圖形的矩形，比如a,bminf,maxf。生成該矩形上n對均勻分布隨機數(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服從標準正態(tài)。因此， P(|1/nsum_j=1AnX_j-EX|t)=2phi(fractVar(X)/n)-1=q。即以 q 的概率誤差小于fractVar(X)/n 。當然，Var(X)也如上面一樣來近似。中心極限定理的另一個用處是模擬某種連續(xù)過程的軌道，即X_n(t)=frac1sqrtnsum_k=ntxi_k,其中 xi_k 為獨立的 0均值隨機變量。當然，它還可以用來近似模擬標準正態(tài)隨機數(shù)。三、減小方差技術(shù)(對立變量、條件期望、控制變量、重要樣本)上面提到，方差往往是不知道，它也是通過產(chǎn)生的隨機數(shù)來估計。由估計的誤差分析知，該方差越小越好，因此給出幾種減小方差的一般方法。

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

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

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

43、 同的。如果打開門的不是主持人，而是另一個參與者，并且當他打開門時發(fā)現(xiàn)什么也沒有，那 么，剩下的兩個門后是轎車的概率才是相等的。計算機模擬為了驗證這一結(jié)果，我們就要比較不改變選擇中獎的幾率和改變選擇中獎的幾率。模擬方法是：我們從0, 1, 2這3個數(shù)中隨機一個為轎車（即中獎號碼），另隨機一個數(shù)為你 的選擇。如果你的選擇與中獎號相同，則計這次為不改變選擇中獎；如果你的選擇不對，則是 改變選擇中獎。分別累積出不改變選擇中獎和改變選擇中獎的次數(shù)，就可以得到不改變選擇中 獎的幾率和改變選擇中獎的幾率了。為了將結(jié)果表示的明顯，我們可以假設(shè)有100扇門，參與的嘉賓選擇了其中的一扇，而主持人隨后把剩下的99

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

45、 山羊這，我們可以通過同樣的隨機模擬過程來模擬這種情況。并可以驗證改變選擇與否對自己 中獎的影響是相同的。當模擬的次數(shù)逐漸的增多時，其模擬值越接近理論值，這說明模擬的效果越好。在隨機事 件的大量重復(fù)出現(xiàn)中，往往呈現(xiàn)幾乎必然的規(guī)律，這個規(guī)律就是大數(shù)定律。通俗地說，這個定 理就是，在試驗不變的條件下，重復(fù)試驗多次，隨機事件的頻率近似于它的概率。偶然必然中 包含著必然。此次模擬試驗也正好用實際的模擬例子說明了大數(shù)定理的正確性和應(yīng)用性。Matlab程序1、編寫函數(shù)n=10000; %實驗次數(shù)stick=0; %堅持選擇的獲獎次數(shù)18alter=0; %改變選擇的獲獎次數(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行為固定坐標軸和柵格的)plot(N,bubian)3、bubianmoni.m這是模擬并返回不改變中獎的概率的程序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; %不變中獎的概率 其中對改變選擇的累計在此處是可以加以省略的。問題二蒲豐投針問題1777年法國科學家蒲豐提出的一種

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

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

50、176620n =100000所模擬出的圓周率為：3.136664可見，模擬的次數(shù)越多，所得的結(jié)果就越準確。這與理論所知的結(jié)果是一樣的。問題三擲骰子問題21 擲三只相同的色子，求有兩個點數(shù)相同的概率。問題分析構(gòu)造3個整數(shù)隨機變量，使之服從1到6之間的均勻分布，比較每次模擬出的 3個值, 如果至少有兩個相同，則將成功的次數(shù)加1,分析可知，成功的頻率即為所求的概率。Matlab程序n=input(,請輸入模擬的次數(shù):)1=0;for i=1:nx=unidrnd(6,1,3); %x 中含有3個元素，相互獨立，且服從 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;%計算成功的頻率，即為概率fprintf(至少兩個相同的概率是%fn,p)運行結(jié)果與分析以下是取4個不同模擬次數(shù)所得的結(jié)果:n =100至少兩個相同的概率是0.470000n =1000至少兩個相同的概率是0.456000n =10000至少兩個相同的概率是 0.437300n =100000至少兩個相同的概率是0.444390由概率知識可知：P=1-(4*5*6)/(6*6*6)=5/9=0.4445。這與以上模擬出的值基本上是一樣的。問題四 無記憶性的例子考慮有兩個營業(yè)員白郵局。假設(shè)當 Smith先生進入郵局的時候，

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

53、X;Y+ S=1/(2* 入)+1/ 入,其中 S 為 Smith 服務(wù)花費的時間,X;Y 分別為Jones和Brown服務(wù)所花費的時間。Matlab程序編寫函數(shù) function mytwo(lambda) n=100000; % 實驗次數(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三個變量記錄各人最后離開次數(shù)統(tǒng)計n次實驗smith花費時間總和分別為三人服務(wù)時間最后離開最后離開最后離開%t=t3+min(t1,t2)轉(zhuǎn)化為概率平均花費時間ETreal=3/(2*lambda) %分析出的ET,用來與模擬結(jié)果對比保存為 mytwo.m,在 command window輸入 my

55、two(3),結(jié)果為 smith =0.5019 jones = 0.2495 brown =0.248623ET =0.5009ETreal =0.5000問題五生日問題有n個人的班級中至少有兩個人生日是同一天的概率，模擬之。問題分析假設(shè)：所有人生日是相互獨立的。設(shè)一年有y天，為求至少兩個人生日在同一天的概率，我們可以先求n個人生日全不同的概率：nn個人生日所有可能的組合(不考慮次序)種數(shù)有：工n!而n個人生日都不一樣的可能的組合方式種數(shù)有：cy = y!(y -n)!* n!故任意兩人生日不同的理論概論是：(這樣寫可以方便編程計算)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注：顯然當學生人數(shù)n比一年的天數(shù)y大時，有抽屜原理易知：p=1.編程思想主程序：在程序初始位置可以定義學生人數(shù)n(默認為50),隨機模擬次數(shù)tn (默認為100),一年的天數(shù)y (默認為365)。通過調(diào)用子函數(shù)計算隨機模擬的概率以及理論概率，將二者畫在 同一個圖中，便于對比分析。子程序1 ： func1.m計算隨機試驗頻率。在每一次隨機試驗中，首先隨機生成1到y(tǒng)的一個隨機數(shù)，對應(yīng)于一個人的生日，根據(jù)人數(shù)存放于數(shù)組a中。使用

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

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

59、機模擬的試驗次數(shù)在增加，當比較小的時候，如圖一tn=30,頻率線一直在理論概率線附近上下波動，而且波動比較大，十分顯著。當試驗次數(shù)逐漸增加的時候，頻率線也逐漸趨近于概率線，當 tn=1000時，兩者幾乎 重合了。當試驗次數(shù)tn越大時，所得的頻率曲線就越精確接近于概率線，但是同時付出的時間代價也是巨大的。當 tn=1000時，運行時間需要約三分鐘。(2)觀察每一幅圖，都有一個共同的特征，隨著人數(shù)的增加，至少兩個人同一天生日的概率也逐漸增大，而且通過概率線可以看出。當人數(shù)為0和1時，p=0,當人數(shù)達到40的時候，至少兩人同一天生日的概率就可以達到90%當人數(shù)達到60或者70時候，概率已經(jīng)達到99%

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