第4講隨機(jī)數(shù)的生成及隨機(jī)變量抽樣省公共課一等獎(jiǎng)全國賽課獲獎(jiǎng)?wù)n件

實(shí)驗(yàn)?zāi)吭囼?yàn)內(nèi)容學(xué)習(xí)主要隨機(jī)變量抽樣方法1、均勻分布U(0,1)隨機(jī)數(shù)產(chǎn)生2、其它各種分布隨機(jī)數(shù)產(chǎn)生方法3、隨機(jī)數(shù)生成實(shí)例4、試驗(yàn)作業(yè)隨機(jī)數(shù)生成及隨機(jī)變量抽樣第1頁隨機(jī)數(shù)生成隨機(jī)數(shù)產(chǎn)生是實(shí)現(xiàn)MC計(jì)算先決條件。而大多數(shù)概率分布隨機(jī)數(shù)產(chǎn)生都是基于均勻分布U(0,1)隨機(jī)數(shù)。

首先，介紹服從均勻分布U(0,1)隨機(jī)數(shù)產(chǎn)生方法。其次，介紹服從其它各種分布隨機(jī)數(shù)產(chǎn)生方法。以及服從正態(tài)分布隨機(jī)數(shù)產(chǎn)生方法。最終，關(guān)于隨機(jī)數(shù)幾點(diǎn)注。第2頁一、均勻分布U(0,1)隨機(jī)數(shù)產(chǎn)生產(chǎn)生均勻分布標(biāo)準(zhǔn)算法在很多高級(jí)計(jì)算機(jī)語言書都能夠看到。算法簡單，輕易實(shí)現(xiàn)。使用者能夠自己手動(dòng)編程實(shí)現(xiàn)。Matlab中也提供給我們用于產(chǎn)生均勻分布各種函數(shù)。我們重點(diǎn)是怎樣經(jīng)過均勻分布產(chǎn)生服從其它分布隨機(jī)數(shù)。所以，直接使用Matlab提供可靠安全標(biāo)準(zhǔn)函數(shù)，當(dāng)然不用費(fèi)事了。

第3頁IMSL庫中函數(shù)使用RNSET:種子設(shè)定CALLRNSET(ISEED)RNOPT:

產(chǎn)生器類型設(shè)定

CALLRNOPT(IOPT)

RNUN/DRNUN:產(chǎn)生均勻分布隨機(jī)數(shù)CALLRNUN(NR,R)

第4頁例1生成1行1000列1—10上離散均勻分布隨機(jī)數(shù)；生成1行1000列21—30上離散均勻分布隨機(jī)數(shù)；生成1行1000列501—1000上離散均勻分布隨機(jī)數(shù)。并畫經(jīng)驗(yàn)分布函數(shù)曲線。Randnum=unidrnd(10,1,10000)；cdfplot(Randnum)；pauseRandnum=unidrnd(10,1,10000)+10；cdfplot(Randnum)；pauseRandnum=unidrnd(500,1,10000)+500；cdfplot(Randnum)cdfplot(x)第5頁第6頁解:由密度函數(shù)知

例2設(shè)總體X密度函數(shù)為其中>0,生成1行10000列隨機(jī)數(shù).含有均值為指數(shù)分布Randnum=exprnd(2,1,10000)+5并畫經(jīng)驗(yàn)分布函數(shù)曲線。cdfplot(Randnum)第7頁第8頁二、其它各種分布隨機(jī)數(shù)產(chǎn)生基本方法有以下三種：逆變換法合成法篩選法第9頁逆變換法設(shè)隨機(jī)變量分布函數(shù)為

，定義

定理

設(shè)隨機(jī)變量服從上均勻分布，則分布函數(shù)為。所以，要產(chǎn)生來自隨機(jī)數(shù)，只要先產(chǎn)生來自

隨機(jī)數(shù)，然后計(jì)算即可。其步驟為

第10頁為常數(shù)例3

設(shè)密度函數(shù)為并畫經(jīng)驗(yàn)分布函數(shù)曲線。第11頁例4

設(shè)X分布函數(shù)為F(X)生成n=201行10000列隨機(jī)數(shù)，并畫經(jīng)驗(yàn)分布函數(shù)曲線。第12頁n=20Randnum=1-(1-unifrnd(0,1,1,10000)).^(1/n);cdfplot(Randnum)第13頁為常數(shù)例5

設(shè)密度函數(shù)為并畫經(jīng)驗(yàn)分布函數(shù)曲線。第14頁合成法

合成法應(yīng)用最早見于Butlter書中。構(gòu)思以下：

假如密度函數(shù)難于抽樣，而關(guān)于條件密度函數(shù)以及密度函數(shù)均易于抽樣，則隨機(jī)數(shù)能夠下產(chǎn)生：能夠證實(shí)由此得到服從。第15頁篩選抽樣

假設(shè)我們要從抽樣，假如能夠?qū)⒈硎境?/p>

，其中是一個(gè)密度函數(shù)且易于抽樣，而，是常數(shù)，則抽樣能夠下進(jìn)行：

定理

設(shè)密度函數(shù)，且，其中，，是一個(gè)密度函數(shù)。令和分別服從和，則在條件下，條件密度為

第16頁三、生成標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)隨機(jī)數(shù)產(chǎn)生方法很多。簡明介紹三種。

法1、

變換法（Box和Muller1958）設(shè)，是獨(dú)立同分布變量，令則與獨(dú)立，均服從標(biāo)準(zhǔn)正態(tài)分布。法2、結(jié)合合成法與篩選法。（略）法3、近似方法（利用中心極限定理）即用個(gè)變量產(chǎn)生一個(gè)變量。其中是抽自隨機(jī)數(shù)，可近似為一個(gè)變量。第17頁例6生成單位圓上均勻分布1行10000列隨機(jī)數(shù)，并畫經(jīng)驗(yàn)分布函數(shù)曲線。Randnum=unifrnd(0,2*pi,1,10000);xRandnum=cos(Randnum)[Y,II]=sort(xRandnum)yRandnum=sin(Randnum)plot(xRandnum(II),yRandnum(II),'.')第18頁第19頁例7生成單位正方形上均勻分布1行10000列隨機(jī)數(shù)，并畫散點(diǎn)圖。mm=10000;Randnum=unifrnd(0,4,1,mm);xRandnum=zeros(1,mm);yRandnum=zeros(1,mm);forii=1:mmifRandnum(1,ii)<=1xRandnum(1,ii)=0;yRandnum(1,ii)=Randnum(1,ii);elseifRandnum(1,ii)<=2xRandnum(1,ii)=Randnum(1,ii)-1;yRandnum(1,ii)=1;elseifRandnum(1,ii)<=3xRandnum(1,ii)=1;yRandnum(1,ii)=1-(Randnum(1,ii)-2);elsexRandnum(1,ii)=1-(Randnum(1,ii)-3);yRandnum(1,ii)=0;endendendend[Y,JJ]=sort(xRandnum);plot(xRandnum(JJ),yRandnum(JJ),'.')第20頁第21頁離散型隨機(jī)變量生成離散型隨機(jī)變量X,它取值是非光滑連續(xù)值，它只能間斷地即離散地取值x1,x2,x3，…，xn，且要求x1＜x2＜x3＜…<xn。其概率密度函數(shù)為p(xi)=p{X=xi}概率分布函數(shù)為例10

對(duì)某車間天天需求某種零件數(shù)量歷史數(shù)據(jù)中統(tǒng)計(jì)取得表1結(jié)果。生成1行1000列零件需求隨機(jī)數(shù)。并畫經(jīng)驗(yàn)分布函數(shù)曲線。表1某零件天天需求量X第22頁需求量x(件)概率P(x)累積概率F(x)可分配隨機(jī)數(shù)范圍X1=10

0.10

F(X1)=0.10（.00--.10]X2=20

0.20

F(X2)=0.30（.10--.30]X3=30

0.40

F(X3)=0.70（.30--.70]X4=40

0.25

F(X4)=0.95

（.70--.95]X5=50

0.05F(X5)=1.00（.95–1)隨機(jī)變量生成算法為①產(chǎn)生一個(gè)u(0，1)，并令i=0；②令i=i+1；③若u＞F(xi)，轉(zhuǎn)回到第②步，不然轉(zhuǎn)至④；④輸出得X＝xi。第23頁mm=10000;Randnum=unifrnd(0,1,1,mm);xRandnum=zeros(1,mm);forii=1:mmifRandnum(1,ii)<=0.1xRandnum(1,ii)=10;elseifRandnum(1,ii)<=0.3xRandnum(1,ii)=20;elseifRandnum(1,ii)<=0.7xRandnum(1,ii)=30;elseifRandnum(1,ii)<=0.95xRandnum(1,ii)=40;elsexRandnum(1,ii)=50;endendendendendcdfplot(xRandnum)第24頁第25頁

三角分布(a,m,b)隨機(jī)變量其密度函數(shù)為其分布函數(shù)為第26頁第27頁在用MonteCarlo等方法解應(yīng)用問題時(shí),隨機(jī)向量抽樣也是經(jīng)慣用到.若隨機(jī)向量各分量相互獨(dú)立,則它等價(jià)于多個(gè)一元隨機(jī)變量抽樣。隨機(jī)向量抽樣方法

第28頁例8生成單位正方形內(nèi)均勻分布1行10000列隨機(jī)數(shù)，并畫散點(diǎn)圖。mm=10000xRandnum=unifrnd(0,1,1,mm);yRandnum=unifrnd(0,1,1,mm);plot(xRandnum,yRandnum,'.')第29頁第30頁第31頁mm=100000xRandnum=unifrnd(0,1,1,mm);yRandnum=unifrnd(0,1,1,mm);[Y,JJ]=sort(xRandnum)plot(xRandnum(JJ),yRandnum(JJ),'.')第32頁例9生成單位圓內(nèi)均勻分布1行10000列隨機(jī)數(shù)，并畫散點(diǎn)圖。mm=10000；Randnum1=unifrnd(-1,1,1,2*mm);Randnum2=unifrnd(-1,1,1,2*mm);xRandnum=zeros(1,mm);yRandnum=zeros(1,mm);s=Randnum1.^2+Randnum2.^2；ii=1;jj=1;whileii<mmifs(1,jj)<=1;xRandnum(1,ii)=Randnum1(1,jj);yRandnum(1,ii)=Randnum2(1,jj);ii=ii+1;endjj=jj+1;endplot(xRandnum,yRandnum,'.')第33頁第34頁關(guān)于隨機(jī)數(shù)幾點(diǎn)注注1

因?yàn)榫鶆蚍植茧S機(jī)數(shù)產(chǎn)生總是采取某個(gè)確定模型進(jìn)行，從理論上講，總會(huì)有周期現(xiàn)象出現(xiàn)。初值確定后，全部隨機(jī)數(shù)也隨之確定，并不滿足真正隨機(jī)數(shù)要求。所以通常把由數(shù)學(xué)方法產(chǎn)生隨機(jī)數(shù)成為偽隨機(jī)數(shù)。注2

應(yīng)對(duì)所產(chǎn)生偽隨機(jī)數(shù)作各種統(tǒng)計(jì)檢驗(yàn)，如獨(dú)立性檢驗(yàn)，分布檢驗(yàn)，功率譜檢驗(yàn)等等。

但其周期又相當(dāng)長，在實(shí)際應(yīng)用中幾乎不可能出現(xiàn)。所以，這種由計(jì)算機(jī)產(chǎn)生偽隨機(jī)數(shù)能夠看成真正隨機(jī)數(shù)來處理。第35頁2.設(shè)密度函數(shù)為1.生成單位球內(nèi)均勻分布1行10000列隨機(jī)數(shù)，并畫散點(diǎn)圖。作業(yè):為常數(shù)并畫經(jīng)驗(yàn)分布函數(shù)曲線。第36頁3.生成三角分布(0,1,2)1行10000列隨機(jī)數(shù)，并畫散點(diǎn)圖。作業(yè):并畫經(jīng)驗(yàn)分布函數(shù)曲線。第37頁

§5.2隨機(jī)數(shù)與隨機(jī)變量生成

§5.2.1隨機(jī)數(shù)生成

在系統(tǒng)模擬中只要有隨機(jī)變量，則在模擬運(yùn)行每一步中都要對(duì)隨機(jī)變量確定一個(gè)詳細(xì)值。我們將會(huì)碰到各種概率分布隨機(jī)變量，但其中最簡單或最基本隨機(jī)變量是在(0，1)區(qū)間上均勻分布隨機(jī)變量。服從某一分布隨機(jī)變量都能夠經(jīng)過對(duì)(0，1)均勻分布隨機(jī)變量進(jìn)行適當(dāng)轉(zhuǎn)換而得到。(0，1)均勻分布隨機(jī)變量取值也是在(0，1)區(qū)間上均勻分布隨機(jī)數(shù)ui序列(流)獨(dú)立采樣，其密度函數(shù)是ui數(shù)學(xué)期望和方差分別為第38頁所以，若能取得(0，1)均勻分布隨機(jī)數(shù)，也就能經(jīng)過對(duì)其適當(dāng)轉(zhuǎn)換而取得某一要求分布隨機(jī)變量取值，這就是隨機(jī)變量生成。為此，首先要掌握(0，1)區(qū)間上均勻分布隨機(jī)數(shù)生成方法。均勻分布隨機(jī)數(shù)必須具備均勻性和獨(dú)立性要求；要生成符合上述要求隨機(jī)數(shù)流，現(xiàn)在多用數(shù)學(xué)算法來產(chǎn)生，普通是采取遞推算法，確定一個(gè)初始值(種子數(shù))以后，逐次遞推算得隨機(jī)數(shù)流。數(shù)學(xué)算法取得隨機(jī)數(shù)、常稱之為偽隨機(jī)數(shù)(PseudoRandomNumber)序列。數(shù)學(xué)方法計(jì)算產(chǎn)生隨機(jī)數(shù)流必須滿足以下要求：(1)盡可能在(0，1)區(qū)間均勻分布；(2)含有統(tǒng)計(jì)上獨(dú)立性；(3)產(chǎn)生隨機(jī)數(shù)流能夠重復(fù)出現(xiàn)，即給以相同初值(種子數(shù))能取得相同隨機(jī)數(shù)流；(4)有足夠長周期，即在出現(xiàn)周期性重復(fù)之前，能生成足夠多個(gè)隨機(jī)數(shù)；(5)算法占用計(jì)算機(jī)內(nèi)存較少而計(jì)算生成速度較快。當(dāng)前廣泛應(yīng)用算法是線性同余法(LinearcongruentialMethod)，其中又分為：1．混合線性同余法。它是由Lehmer于1951年提出，其算式為第39頁xi+1=(axi+c)modmui+1=xi+1/m式中a——乘數(shù)(常數(shù))；C——增量(常數(shù))；x0——種子數(shù)；m——模數(shù)。a,c,m和x0選取對(duì)隨機(jī)數(shù)流統(tǒng)計(jì)特征和周期長度有極大影響。

上述第一式含義是式中[]表示取整數(shù)，a，c，m皆為整常數(shù)。2、

乘法線性同余法。若混合線性同余法中c=0，則為乘法線性同余法，其算式為

xi+1=aximodmui+1=xi+1/m(5.7)第40頁可參考選取數(shù)據(jù)有：（1）

a=16807,m=2147483647,x0=123457;（2）

a=655393,m=33554432。

§5.2.2隨機(jī)數(shù)流檢驗(yàn)

一、均勻分布性檢驗(yàn)1．參數(shù)檢驗(yàn)。檢驗(yàn)ui數(shù)字特征，如均值、方差預(yù)計(jì)值和其理論值差異是否顯著。設(shè)有u1,u2,…，un隨機(jī)數(shù)流，則它們?nèi)魎i序列在(0，1)上均勻分布，可假設(shè)：u期望和方差分別為第41頁2期望和方差分別為則上列假設(shè)(8．10)與(8．11)應(yīng)該成立。據(jù)此，可對(duì)n個(gè)ui計(jì)算以下統(tǒng)計(jì)量若取顯著性水平a=0．05，

當(dāng)|V1|≤1．96時(shí)。則可認(rèn)為假設(shè)(8．10)式成立；當(dāng)|V2|≤1．96時(shí)，則可認(rèn)為假設(shè)(8．11)式成立。因而能夠接收此假設(shè)，檢驗(yàn)經(jīng)過;不然拒絕接收。2．均勻性檢驗(yàn)。它是檢驗(yàn)所生成隨機(jī)數(shù)落在(0，1)各子區(qū)間頻率均勻程度，是否與理論上均勻分布頻率有顯著性差異。此處介紹慣用方法之一，x2檢驗(yàn)方法以下：

將(0，1)區(qū)間劃分為相等k個(gè)子區(qū)間，假如落在第i個(gè)(i＝l，2，3，…，k)子區(qū)間隨機(jī)數(shù)有ni個(gè)；而在理論上第i個(gè)子區(qū)間隨機(jī)數(shù)個(gè)數(shù)為mi=N／k，其中N為隨機(jī)數(shù)流總個(gè)數(shù)(擬檢驗(yàn))。由此，可計(jì)算x2統(tǒng)計(jì)量第42頁再按k-1為自由度、顯著性水平取0．05，查得

2(a)表值。當(dāng)算得統(tǒng)計(jì)量x2≤

2(a)時(shí)，可認(rèn)為在顯著性水平a下能接收為均勻分布假設(shè)。

二、獨(dú)立性檢驗(yàn)獨(dú)立性檢驗(yàn)是檢驗(yàn)隨機(jī)數(shù)流中前后各數(shù)之間是否存在相關(guān)性。慣用方法是進(jìn)行自相關(guān)檢驗(yàn)。另外，還有PokerTest和Run檢驗(yàn)，普通應(yīng)用較少。

§5.2.3隨機(jī)變量生成

一、離散型隨機(jī)變量生成離散型隨機(jī)變量X,它取值是非光滑連續(xù)值，它只能間斷地即離散地取值x1,x2,x3，…，xn，且要求x1＜x2＜x3＜…<xn。其概率密度函數(shù)為p(xi)=p{X=xi}

概率分布函數(shù)為第43頁譬如，對(duì)某車間天天需求某種零件數(shù)量歷史數(shù)據(jù)中統(tǒng)計(jì)取得表5．2結(jié)果。表5．2某零件天天需求量X需求量x(件)概率P(x)累積概率F(x)可分配隨機(jī)數(shù)范圍X1=100.10F(X1)=0.10.00--.09X2=200.20F(X2)=0.30.10--.29X3=300.40F(X3)=0.70.30--.69X4=400.25F(X4)=0.95.70--.94X5=500.05F(X5)=1.00.95--.99隨機(jī)變量生成算法為①產(chǎn)生一個(gè)u(0，1)，并令i=0；②令i=i+1；③若u＞F(xi)，轉(zhuǎn)回到第②步，不然轉(zhuǎn)至④；④輸出得X＝xi。第44頁二、連續(xù)型隨機(jī)變量生成

以下介紹慣用幾個(gè)概率分布隨機(jī)變量生成方法。1．逆變換法。任何概率分布隨機(jī)變量X，其分布函數(shù)F(x)值域是(0,l)；所以，可令

F(x)=u(5．17)假如對(duì)上式能解出顯式逆函數(shù)

x=F-1(u)(5．18)則調(diào)用隨機(jī)數(shù)生成算法產(chǎn)生一個(gè)隨機(jī)數(shù)u，將其代入顯式(8．18)，即可取得隨機(jī)變量X一次取值x?？上芮蟮蔑@式逆函數(shù)隨機(jī)變量，只有以下幾個(gè)分布。(l)在(a，b)區(qū)間均勻分布隨機(jī)變量X。其概率密度函數(shù)為其分布函數(shù)為第45頁均值為(a十b)／2；方差為(b—a)2／12。兩函數(shù)圖形如圖8．1。令u=(x-a)／(b-a)，則可解得x=a+(b–a).u(5．21)(2)三角分布(a，m，b)隨機(jī)變量其密度函數(shù)為其分布函數(shù)為

第46頁第47頁第48頁第49頁kUP結(jié)果00.4430.443P>0.002510.1720.076P>0.002520.4230.032P>0.002530.8190.026P>0.002540.2550.0066P>0.002550.7490.0049P>0.00250.2250.0011P<0.0025,置X=6kUP結(jié)果00.4430.443P>0.002510.1720.076P>0.0025