1、實驗目的實驗目的實驗內容實驗內容學習學習主要的隨機變量抽樣方法主要的隨機變量抽樣方法1 1、均勻分布均勻分布U(0,1)的隨機數的產生的隨機數的產生2 2、其他各種分布的隨機數的產生方法其他各種分布的隨機數的產生方法3 3、隨機數生成實例、隨機數生成實例4 4、實驗作業、實驗作業隨機數的生成及隨機變量抽樣隨機數的生成隨機數的生成 隨機數的產生是實現MC計算的先決條件。而大多數概率分布的隨機數的產生都是基于均勻分布U(0,1)的隨機數。 首先，介紹服從均勻分布U(0,1)的隨機數的產生方法。 其次，介紹服從其他各種分布的隨機數的產生方法。以及服從正態分布的隨機數的產生方法。 最后，關于隨機數的幾

2、點注。一、一、均勻分布均勻分布U(0,1)的隨機數的產生的隨機數的產生 產生均勻分布的標準算法在很多高級計算機語言的書都可以看到。算法簡單，容易實現。使用者可以自己手動編程實現。Matlab 中也提供給我們用于產生均勻分布的各種函數。我們的重點是怎樣通過均勻分布產生服從其他分布的隨機數。因此，直接使用Matlab提供的可靠安全的標準函數，當然不用費事了。 IMSL庫中的函數使用庫中的函數使用 RNSET: 種子的設定 CALL RNSET (ISEED) RNOPT: 產生器的類型的設定 CALL RNOPT (IOPT) RNUN/DRNUN: 產生均勻分布的隨機數 CALL RNUN (N

3、R, R) 例例1 1生成生成1 1行行10001000列的列的1 11010上離散均勻分布的隨機上離散均勻分布的隨機數；數；生成生成1 1行行10001000列列21213030上離散均勻分布的隨機數；上離散均勻分布的隨機數；生成生成1 1行行10001000列列50150110001000上離散均勻分布的隨機上離散均勻分布的隨機數。數。 并畫經驗分布函數曲線。并畫經驗分布函數曲線。Randnum=unidrnd(10,1,10000)；cdfplot(Randnum)；pauseRandnum=unidrnd(10,1,10000)+10；cdfplot(Randnum)；pauseRan

4、dnum=unidrnd(500,1,10000)+500；cdfplot(Randnum)cdfplot(x)解解:由密度函數知由密度函數知 例例2設總體設總體X的密度函數為的密度函數為為未知參數其它 , 0,1)()(xexfXx其中其中 0, 生成生成 1行行10000列的隨機數列的隨機數. X具有均值為具有均值為 的指數分布的指數分布 Randnum=exprnd(2,1,10000)+55, 2并畫經驗分布函數曲線。并畫經驗分布函數曲線。cdfplot(Randnum)二、其他各種分布的隨機數的產生二、其他各種分布的隨機數的產生 基本方法有如下三種： 逆變換法 合成法 篩選法 逆變換

5、法逆變換法 設隨機變量 的分布函數為 ，定義 定理定理 設隨機變量 服從 上的均勻分布，則 的分布函數為 。 因此，要產生來自 的隨機數，只要先產生來自 的隨機數，然后計算 即可。 其步驟為 X xF 10 ,:inf1yyxFxyFU) 1 , 0( UFX1 xF xF1 , 0U uF1 uFxuU11 , 0計算，抽取由的隨機數生成,min21nXXXY00, 01);(xxexfx0為常數例例3 設密度函數為并畫經驗分布函數曲線。并畫經驗分布函數曲線。0, 00,1)(0, 00,1)(yyeyFxxexFnyYxX，的隨機數生成,min21nXXXY例例4 設X分布函數為F(X)1

6、 (1 (,)(11)(11nXYnXYUFRandyFyF，0, 001,)(xxxxFXnnXYnYUUFRandyyF111)1 (1)1 (1 (11)(生成生成n=20的的1行行10000列隨機數，并畫經驗分布列隨機數，并畫經驗分布函數曲線。函數曲線。n=20Randnum=1-(1-unifrnd(0,1,1,10000).(1/n);cdfplot(Randnum)的隨機數生成X1n為常數例例5 設密度函數為并畫經驗分布函數曲線。并畫經驗分布函數曲線。1,.,00, 0,)(1nixxxxcxfiiiX，其它合成法合成法 合成法的應用最早見于Butlter 的書中。 構思如下：

7、如果 的密度函數 難于抽樣，而 關于 的條件密度函數 以及 的密度函數 均易于抽樣，則 的隨機數可如下產生： 可以證明由此得到 的服從 。X xpXYyxpY ygX xyxpyygY抽取由條件分布，抽取的分布由X xp篩選抽樣篩選抽樣 假設我們要從 抽樣，如果可以將 表示成 ，其中 是一個密度函數且易于抽樣，而 ， 是常數，則 的抽樣可如下進行： 定理定理 設 的密度函數 ，且 ，其中 ， ， 是一個密度函數。令 和 分別服從 和 ，則在 的條件 下， 的條件密度為 xp xgxhcxp xph 10 xg1cX 。，回到如果，停止，則如果，抽取，由抽取由1321 , 01yguyxyguy

8、yhuUX xp xgxhcxph 10 xg1cUY1 , 0U yh YgU xpYgUxpYY三、生成標準正態分布的隨機數三、生成標準正態分布的隨機數 的隨機數產生方法很多。簡要介紹三種。 法法1、 變換法（Box 和Muller 1958） 設 ， 是獨立同分布的 變量，令 則 與 獨立，均服從標準正態分布。 法法2、 結合合成法與篩選法。（略） 法法3 3、 近似方法（利用中心極限定理） 即用 個 變量產生一個 變量。 其中 是抽自 的隨機數， 可近似為一 個 變量。1 , 0N1U2U1 , 0U2122112sinln22cosln22121UUXUUX1X2Xn1 , 0U1

9、, 0N2112unxiu1 , 0Uniiunu111 , 0N例例6生成單位圓上均勻分布的生成單位圓上均勻分布的1行行10000列隨機數，列隨機數，并畫經驗分布函數曲線。并畫經驗分布函數曲線。Randnum=unifrnd(0,2*pi,1,10000);xRandnum=cos(Randnum)Y,II =sort(xRandnum)yRandnum=sin(Randnum)plot(xRandnum(II),yRandnum(II),.)例例7 生成單位正方形上均勻分布的生成單位正方形上均勻分布的1行行10000列隨列隨機數，并畫散點圖。機數，并畫散點圖。mm=10000;Randnu

10、m=unifrnd(0,4,1,mm);xRandnum=zeros(1,mm);yRandnum=zeros(1,mm);for ii=1:mm if Randnum(1,ii)=1 xRandnum(1,ii)=0; yRandnum(1,ii)=Randnum(1,ii); else if Randnum(1,ii)=2 xRandnum(1,ii)=Randnum(1,ii)-1; yRandnum(1,ii)=1; else if Randnum(1,ii)=3 xRandnum(1,ii)=1; yRandnum(1,ii)=1-(Randnum(1,ii)-2); else xR

11、andnum(1,ii)=1-(Randnum(1,ii)-3); yRandnum(1,ii)=0; end end end end Y,JJ =sort(xRandnum);plot(xRandnum(JJ),yRandnum(JJ),.)離散型隨機變量的生成離散型隨機變量的生成離散型隨機變量X,它的取值是非光滑連續的值，它只能間斷地即離散地取值x1,x2,x3，xn，且規定x1x2x3xn。其概率密度函數為 p(xi)=pX= xi概率分布函數為xixiixpxXpXF)()(例例1010 對某車間每天需求某種零件的數量歷史數據中統計獲得表1的結果。生成1行1000列零件需求的隨機數。并

12、畫經驗分布函數曲線。并畫經驗分布函數曲線。 表1 某零件每天需求量 X需求量x(件)概率P(x)累積概率F(x)可分配的隨機數范圍X1=10 0.10 F(X1)=0.10（ .00 - .10X2=20 0.20 F(X2)=0.30（ .10 - .30X3=30 0.40 F(X3)=0.70（ .30 - .70X4=40 0.25 F(X4)=0.95 （ .70 - .95X5=50 0.05 F(X5)=1.00（ .95 1)隨機變量生成的算法為產生一個u(0，1)，并令i=0；令i=i+1；若uF(xi)，轉回到第步，否則轉至；輸出得 Xxi。mm=10000;Randnum

13、=unifrnd(0,1,1,mm);xRandnum=zeros(1,mm);for ii=1:mm if Randnum(1,ii)=0.1 xRandnum(1,ii)=10; else if Randnum(1,ii)=0.3 xRandnum(1,ii)=20; else if Randnum(1,ii)=0.7 xRandnum(1,ii)=30; else if Randnum(1,ii)=0.95 xRandnum(1,ii)=40; else xRandnum(1,ii)=50; end end end end end cdfplot(xRandnum) 三角分布(a,m,b

14、)的隨機變量其密度函數為 bx 1 a)-a)(b-a)/(m-2(x 0)(時當時當時或當mmxabxxaxf bx 1bx )(/()(1 a)-a)(b-/(ma)-(x 0)(22時當時當時當時當mabmbxbmxaaxxF其分布函數為 在用Monte Carlo等方法解應用問題時,隨機向量的抽樣也是經常用到的. 若隨機向量各分量相互獨立,則它等價于多個一元隨機變量的抽樣。隨機向量的抽樣方法例例8生成單位正方形內均勻分布的生成單位正方形內均勻分布的1行行10000列隨機列隨機數，并畫散點圖。數，并畫散點圖。mm=10000 xRandnum=unifrnd(0,1,1,mm);yRan

15、dnum=unifrnd(0,1,1,mm);plot(xRandnum,yRandnum,.)mm=100000 xRandnum=unifrnd(0,1,1,mm);yRandnum=unifrnd(0,1,1,mm);Y,JJ =sort(xRandnum)plot(xRandnum(JJ),yRandnum(JJ),.)例例9 生成單位圓內均勻分布的生成單位圓內均勻分布的1行行10000列隨機數，列隨機數，并并畫散點圖。畫散點圖。mm=10000；Randnum1=unifrnd(-1,1,1,2*mm);Randnum2=unifrnd(-1,1,1,2*mm);xRandnum=z

16、eros(1,mm);yRandnum=zeros(1,mm);s=Randnum1.2+Randnum2.2；ii=1;jj=1;while iimm if s(1,jj)=1; xRandnum(1,ii)=Randnum1(1,jj); yRandnum(1,ii)=Randnum2(1,jj); ii=ii+1; end jj=jj+1;endplot(xRandnum,yRandnum,.)關于隨機數的幾點注關于隨機數的幾點注 注注1 1 由于均勻分布的隨機數的產生總是采用某個確定的模型進行的，從理論上講，總會有周期現象出現的。初值確定后，所有隨機數也隨之確定，并不滿足真正隨機數的要

17、求。因此通常把由數學方法產生的隨機數成為偽隨機數。 注注2 2 應對所產生的偽隨機數作各種統計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。 但其周期又相當長，在實際應用中幾乎不可能出現。因此，這種由計算機產生的偽隨機數可以當作真正的隨機數來處理。2.設設密度函數為密度函數為1.1.生成單位球內均勻分布的生成單位球內均勻分布的1 1行行1000010000列隨機數，列隨機數，并畫散點圖。并畫散點圖。作業作業:00, 01);(xxexfx0為常數的隨機數生成,max21nXXXY并畫經驗分布函數曲線。并畫經驗分布函數曲線。3.3.生成三角分布生成三角分布(0,1,2)(0,1,2)的的1 1行

18、行1000010000列隨機數，并列隨機數，并畫散點圖。畫散點圖。作業作業:并畫經驗分布函數曲線。并畫經驗分布函數曲線。 5.2隨機數與隨機變量的生成隨機數與隨機變量的生成 5.2.1隨機數的生成隨機數的生成 在系統模擬中只要有隨機變量，則在模擬運行的每一步中都要對隨機變量確定一個具體的值。我們將會遇到各種概率分布的隨機變量，但其中最簡單或最基本的隨機變量是在(0， 1)區間上均勻分布的隨機變量。服從某一分布的隨機變量都可以通過對(0， 1)均勻分布的隨機變量進行適當轉換而得到。(0， 1)均勻分布的隨機變量的取值也是在(0，1)區間上均勻分布的隨機數ui序列(流)的獨立采樣，其密度函數是(5

19、.1) 010 1)(其它時當xxf(5.2) 1 x 1 10 )(0)(時當時當xxdxxfxxFui的數學期望和方差分別為(5.3) 21022)(0)(2xdxxfxuE(5.4) 121)()(0)(2xEdxxfxxV 因此，若能獲得(0，1)均勻分布的隨機數，也就能通過對其適當的轉換而獲得某一規定分布的隨機變量的取值，這就是隨機變量的生成。為此，首先要掌握(0，1)區間上均勻分布隨機數的生成方法。 均勻分布隨機數必須具備均勻性和獨立性的要求；要生成符合上述要求的隨機數流，現在多用數學算法來產生，一般是采用遞推算法，確定一個初始值(種子數)以后，逐次遞推算得隨機數流。數學算法獲得的

20、隨機數、常稱之為偽隨機數(Pseudo Random Number)序列。 數學方法計算產生的隨機數流必須滿足下列要求： (1)盡可能在(0，1)區間均勻分布； (2)具有統計上的獨立性； (3)產生的隨機數流能夠重復出現，即給以相同的初值(種子數)能獲得相同的隨機數流； (4)有足夠長的周期，即在出現周期性重復之前，能生成足夠多個的隨機數； (5)算法占用計算機內存較少而計算生成速度較快。 目前廣泛應用的算法是線性同余法( Linear congruential Method)，其中又分為： 1混合線性同余法。它是由Lehmer于1951年提出的，其算式為xi+1=(axi+c) mod m

21、 ui+1=xi+1/m式中 a乘數(常數)； C增量(常數)； x0種子數； m模數。 a,c,m和x0的選取對隨機數流的統計特性和周期長度有極大影響。 上述第一式的含義是(5.6) m )(1mcaxcaxxiii式中表示取整數， a，c，m皆為整常數。 2、 乘法線性同余法。若混合線性同余法中c=0，則為乘法線性同余法，其算式為 x i+1=ax imod m u i+1=xi+1/m (5.7)可參考選用的數據有：（1） a=16807, m=2147483647, x0=123457;（2） a=655393, m=33554432。 5.2.2隨機數流的檢驗隨機數流的檢驗 一、均勻

22、分布性檢驗一、均勻分布性檢驗 1參數檢驗。檢驗ui的數字特征，如均值、方差的估計值和其理論值的差異是否顯著。 設有u1, u2,，un隨機數流，則它們的 (5.8) 1 1niiunu均值(5.9) )(11 22uuni方差若 ui序列在(0，1)上均勻分布，可假設：u的期望和方差分別為 (5.10) 121)( 21)(nuVuE2的期望和方差分別為則上列假設(810)與(811)應該成立。據此，可對n個ui計算下列統計量若取顯著性水平a=005， (5.11) 1801) ( 121)(22nVE(5.12) )21(12)()(1unuVuEuV(5.13) )21 (180) ()

23、( 22222nVEV 當|V1|196時。則可認為假設(810)式成立； 當|V2|196時，則可認為假設(811)式成立。因而可以接受此假設，檢驗通過;否則拒絕接受。 2均勻性檢驗。它是檢驗所生成的隨機數落在(0， 1)各子區間的頻率的均勻程度，是否與理論上的均勻分布頻率有顯著性差異。此處介紹常用方法之一，x2檢驗方法如下： 將(0，1)區間劃分為相等的k個子區間，假如落在第i個( il，2，3，k)子區間的隨機數有ni個；而在理論上第i個子區間的隨機數個數為mi = Nk，其中 N為隨機數流總個數(擬檢驗的)。由此，可計算 x2統計量(5.14) )()(21122KiiKiiiiKNn

24、NKmmnx再按k-1為自由度、顯著性水平取005，查得2(a)表值。當算得統計量x22(a)時，可認為在顯著性水平a下能接受為均勻分布假設。 二、獨立性檢驗二、獨立性檢驗獨立性檢驗是檢驗隨機數流中前后各數之間是否存在相關性。常用的方法是進行自相關檢驗。此外，還有Poker Test和Run檢驗，一般應用較少。 5.2.3隨機變量的生成隨機變量的生成 一、離散型隨機變量的生成一、離散型隨機變量的生成離散型隨機變量 X,它的取值是非光滑連續的值，它只能間斷地即離散地取值x1,x2,x3，xn，且規定x1x2x30.0025 1 0.172 0.076 P0.0025 2 0.423 0.032 P0.0025 3 0.819 0.026 P0.0025 4 0.255 0.0066 P0.0025 5 0.749 0.0049 P0.0025 6 0.225 0.0011 P0.0025,置 X=6 3近似計算法近似計算法。