第四章隨機(jī)數(shù)產(chǎn)生原理

1、第四章 隨機(jī)數(shù)產(chǎn)生原理,一、引言二、偽隨機(jī)數(shù)產(chǎn)生原理三、[0，1]均勻分布隨機(jī)數(shù)的算法四、其他分布隨機(jī)數(shù)的產(chǎn)生五、正態(tài)分布隨機(jī)數(shù)的產(chǎn)生六、MATLAB統(tǒng)計庫中的隨機(jī)數(shù)發(fā)生器七、隨機(jī)數(shù)的檢驗八、案例3九、習(xí)題,一、引言,以隨機(jī)數(shù)產(chǎn)生為基礎(chǔ)的Monte-Carlo方法已成為現(xiàn)代科研的重要手段之一。其意義早以超出了概率論與數(shù)理統(tǒng)計的范疇。廣泛應(yīng)用于計算方法、隨機(jī)數(shù)規(guī)劃、管理科學(xué)、物理化學(xué)中的高分子結(jié)構(gòu)的研究等領(lǐng)域。我們來看一些

2、例子。,1、數(shù)值計算的研究,數(shù)值計算的研究可以說是一切Monte-Carlo應(yīng)用的基礎(chǔ)，在計算數(shù)學(xué)領(lǐng)域我們遇到了很多的復(fù)雜計算，一個典型的例子是計算積分。對于一維、二維的問題早已獲得解決。但是當(dāng)遇到高維積分問題時，我們傳統(tǒng)的數(shù)值方法都由于計算量太大而陷于了困境。但是高維積分問題又偏偏是物理、高分子化學(xué)、運(yùn)籌學(xué)和最優(yōu)化問題迫切而必須解決的問題。我們來看一個例子。,,這里,這是一個眾所周知的積分公式，我們當(dāng)然也可以把它一般的看為是一個高維積

3、分，如果從傳統(tǒng)的數(shù)值計算方法來看待，則高維問題是隨著維數(shù)的增加計算量成指數(shù)增加的，計算很快就失去控制。但是如果我們換一個角度來看待這個問題，從概率論的角度，它實際是：,,即是f（x）的均值，對于均值我們有一個很好的估計，即,,【例4.1.1】 用Monte-Carlo 對 積分,解：將積分區(qū)域和值域看成是一個邊長為一的正方形。利用均勻分布隨機(jī)數(shù)將點撒在正方形中，計算小于函數(shù)的個數(shù)并除全部點數(shù)。這就

4、是積分的近似值。,% 利用Monte-Carlo方法計算定積分x=rand(1,1000);x_2=x.^2;JF=mean(x_2)% 作Monte-Carlo積分示意圖for i=1:1000 xx=rand(1,100); yy=rand(1,100);endx1=linspace(0,1,50);y1=x1.^2;plot(xx,yy,'.',x1,y1,'linewid

5、th',2)axis equalh = legend('x-y','x^2');JF = 0.3346,面積計算結(jié)果為：s = 0.3482,【例4.1.2】 利用Monte-Carlo方法計算定積分。,,解：抽兩組隨機(jī)數(shù)，求每組元素的平方代入給定的函數(shù)，然后求平均值即得積分的近似值。% Monte-Carlo方法積分二重積分并與數(shù)值方法的結(jié)果進(jìn)行比較Q = dblquad

6、('sin(x.^2+y.^2)',0,1,0,1) % 數(shù)值求積分命令x=rand(2,100000); % 產(chǎn)生兩組隨機(jī)數(shù)Sin_xy= sin(x(1,:).^2+x(2,:).^2); % 代入函數(shù)JF_Sin_xy=mean(Sin_xy) % 用Monte-Carlo方法求積分值,計算結(jié)果為：Q = 0.5613JF_Sin_xy =

7、 0.5612當(dāng)抽樣數(shù)很大時結(jié)果很接近。我們可以從Monte-Carlo方法中看出，如果維數(shù)增加實際計算難度并不增加，因此是處理高維問題的有效方法。,這里 x 是積分定義域中的均勻分布隨機(jī)數(shù)，這是革命性的一個視角。從這個視角，我們把繁難的積分計算變成了簡單的算術(shù)平均，而大量的抽樣計算又是計算機(jī)的拿手好戲，更重要的是當(dāng)維數(shù)增加時并不增加計算難度，從而用 Monte-Carlo 方法研究高維積分問題已是當(dāng)今計算數(shù)學(xué)界的熱門課題。,2、管

8、理科學(xué)的系統(tǒng)仿真研究,管理科學(xué)中的系統(tǒng)仿真研究，如服務(wù)系統(tǒng)、庫存系統(tǒng)等。其共性就是研究的對象是隨機(jī)數(shù)，如顧客到達(dá)時間一般是一個服從指數(shù)分布的隨機(jī)數(shù)，而服務(wù)時間也可以看成是服從某種分布的隨機(jī)數(shù)，當(dāng)一個系統(tǒng)是多隊多服務(wù)體系時，問題就變的相當(dāng)復(fù)雜了。我們很難用數(shù)學(xué)的解析式來表達(dá)。這時Monte-Carlo方法也是有利的武器。對于這個領(lǐng)域的已有各種比較成熟的專用軟件如GPSS、SIMULATION等可以使用。,,3. 物理化學(xué)中的分子領(lǐng)域,50

9、年代科學(xué)家已經(jīng)在高分子領(lǐng)域使用Monte-Carlo方法了。這一領(lǐng)域所研究的問題更加復(fù)雜，計算量非常之大。高分子材料是由幾乎是無窮的高分子鏈組成，而每一個鏈的長度又是10的好幾次方。而鏈的構(gòu)象又是千差萬別，而且是隨機(jī)游動的。如何在中微觀上幾乎是無規(guī)律的現(xiàn)象中去判斷其宏觀的性質(zhì)？用數(shù)學(xué)的解析式來說明這樣的現(xiàn)象是蒼白無力的。Monte-Carlo方法是一個很好的工具，它使得科學(xué)家用Monte-Carlo方法去探索高分子運(yùn)動的規(guī)律。一個典型的

10、例子是：對于高分子多鏈體的研究這是一個很復(fù)雜的問題，直到標(biāo)度理論和重正化群理論方法的引入，才使得單鏈構(gòu)象統(tǒng)計問題得到了較好的解決。,例：用計算機(jī)模擬高分子鏈,鏈的末端距,末端距：空間一鏈的末端與始端的距離為末端距，由于我們將始端放在坐標(biāo)原點，所以末端距的 計算公式為：末端距=（X2+Y2+Z2)1/2這里X，Y，Z為鏈的末端點的坐標(biāo)。顯然末端距隨鏈的不同而不同，即為隨機(jī)變量。,三根鏈的起點(0,0,0),藍(lán)鏈的末端距,綠鏈的

11、末端距,紅鏈的末端距,二、偽隨機(jī)數(shù)產(chǎn)生原理,前面Monte-Carlo方法的例子是以高質(zhì)量的隨機(jī)數(shù)為基礎(chǔ)的。通過完全的隨機(jī)抽樣或調(diào)查可以產(chǎn)生隨機(jī)序列。當(dāng)我們用Monte-Carlo方法研究一個實際問題時，我們需要快速地獲得大量的隨機(jī)數(shù)。用計算機(jī)產(chǎn)生這樣的隨機(jī)數(shù)是非常方便的，用數(shù)學(xué)方法在計算機(jī)上產(chǎn)生的隨機(jī)數(shù)稱為偽隨機(jī)數(shù)。,,基本定理：如果隨機(jī)變量X的分布函數(shù)F（x）連續(xù)，則R = F（x）是[0，1]上的均勻分布的隨機(jī)變量。證

12、：因為分布函數(shù)F（x）是在（0，1）上取值單調(diào)遞增的連續(xù)函數(shù)，所以當(dāng)X在（－∞，x）內(nèi)取值時，隨機(jī)變量R則在（0，F(xiàn)(x)）上取值，且對應(yīng)于（0，1）上的一個R值，至少有一個x滿足，見圖4.2.1,,以 表示隨機(jī)變量 R 的分布函數(shù)，則有,,,(4.2.1),,,=,證畢,,圖4.2.1,(4.2.2),基本定理給出了任一隨機(jī)變量和均勻分布R之間的關(guān)系。而有些隨機(jī)變量可以通過分布函數(shù)的逆變換來獲得，因此如果

13、我們可以產(chǎn)生高質(zhì)量的均勻分布，我們就可以通過變換獲得高質(zhì)量的其他分布。見公式 (4.2.3),,(4.2.3),例4.2.1 求指數(shù)分布的隨機(jī)數(shù)。令,,,,從而我們用服從[0，1]上的隨機(jī)數(shù)R，通過上面的公式就可以得到指數(shù)分布的隨機(jī)數(shù)了。,例4.2.1 產(chǎn)生1000個均勻分布隨機(jī)數(shù)，利用變換產(chǎn)生λ=6的指數(shù)分布并進(jìn)行擬合優(yōu)度檢驗。,clc,clearx=linspace(0,20,100);R=rand(1,1000);

14、 % 產(chǎn)生1000個（0，1）均勻隨機(jī)數(shù)ex=-6*log(1-R); % 變換為指數(shù)分布隨機(jī)數(shù)ex=ex';m=mean(ex)v=var(ex)subplot(1,2,1),cdfplot(ex)subplot(1,2,2),hist(ex)kstest(ex, [ex expcdf(ex, 6)]) % 擬合優(yōu)度檢驗,結(jié)果為：H=0, 接受原假設(shè)，變換后的確為

15、λ=6的指數(shù)分布,三、 （0，1）均勻分布隨機(jī)數(shù)的產(chǎn)生,1、算法要求（1）       產(chǎn)生的數(shù)值序列要具有均勻總體簡單子樣的一些概率統(tǒng)計特性，通常包括分布的均勻性，抽樣的隨機(jī)性，試驗的獨立性和前后的一致性。（2） 產(chǎn)生的隨機(jī)數(shù)要有足夠長的周期。（3） 產(chǎn)生隨機(jī)數(shù)速度快，占用內(nèi)存小。,為了達(dá)到快速的要求，一般采用遞推公式,,（4.3.1),目前最常用的方法是上述方法的

16、一個特例：,混合同余法,,(4.3.2),其中a，b，M 以及初值 y 都是正整數(shù)，容易看出 x 滿足0≤x≤1。其中 mod M 運(yùn)算定義為：任一整數(shù) y 可唯一表示為公式,,,則,乘同余法,當(dāng) b = 0 時，有,,(4.3.4),加同余法,以下形式的同余法稱為加同余法,,(3.4.5),例4.3.1 歷史上比較有名的稱為“菲波那西”數(shù)列為加同余法 的特例。,,(4.3.6),當(dāng) M=8 時，取初值

17、得“菲波那西”數(shù)列。 0，1，1，2，3，5，8，13，21，34，55，89，144，253 ……對上述數(shù)列取模得： 0，1，1，2，3，5，0，5，5，7，1，1…… (4.3.7) 再除以模 M 我們可得到 (0,1) 之間的序列 。,我們知道對于一個來自均勻分布的隨機(jī)序列，應(yīng)該滿足獨立性、均勻性等統(tǒng)計特性，但偽隨機(jī)數(shù)往往有一些缺陷，例如 (4.3.7) 序列到一定長度后，又開始重復(fù)下面

18、的序列，這稱為周期性是一種明顯的規(guī)律，與隨機(jī)性矛盾。通常我們只能選用一個周期內(nèi)的序列作為我們的偽隨機(jī)數(shù)。因此研究一種算法，使得其產(chǎn)生的隨機(jī)序列的周期盡可能長，我們可以通過調(diào)節(jié)（4.3.1）的參數(shù)來實現(xiàn)。因此如何來獲得一個周期比較長的序列，就成了我們研究的一個內(nèi)容：有關(guān)偽隨機(jī)數(shù)序列的周期有如下的一些結(jié)論：,定理 4.3.1 混合同余法產(chǎn)生序列達(dá)最大周期 M 的充要條件： （1） b 與 M 互素 （2）

19、 對于 M 的任意素因子 p，有 （3） 如果 4 是 M 的因子，則 顯然乘同余法產(chǎn)生的序列達(dá)不到周期 M（不滿足（1））。當(dāng)取 （k為任意整數(shù)）時，因為 M 只有一個素因子2， 且4是 M 的因子，則由條件（2）、（3）有 ，從而混合同余發(fā)生器達(dá)到最大周期的算法為：,,,,,（3.4.8),其中c，d為任

20、意整數(shù)?；旌贤喟l(fā)生器是否達(dá)到最大周期M與初始值無關(guān)。,,對于乘同余發(fā)生器，由同余運(yùn)算的定義，知其由如下性質(zhì),(1) 如果,則有：,,,（2）如果,,則,其中（c，M）是c，M的最大公約數(shù)。,利用這些性質(zhì)可得到以下定理。定理4.3.2 對乘同余發(fā)生器，若 ，則使成立得最小正整數(shù) V 就是此發(fā)生器得周期。,,,在數(shù)論中稱 V 為 a 關(guān)于 M 的階數(shù)，對于乘同余發(fā)生器，若種子與 M

21、互素，則其周期就是關(guān)于 M 的階數(shù)。這樣一來，尋找達(dá)到最大周期的同余發(fā)生器的問題就轉(zhuǎn)化為數(shù)論方面尋求 M達(dá)到最大階數(shù) a 的問題了。Knuth 對這一問題的研究作了總結(jié)。,從算法上我們知道公式是遞推的，因此一般的隨機(jī)發(fā)生器程序都要預(yù)先賦初值，這種初值為種（Seed），有些統(tǒng)計軟件如SPSS要求用戶給出Seed.以均勻分布（0，1）隨機(jī)變量R變換成的隨機(jī)變量。以r,ε,u,分別表示 (0,1) 均勻分布,指數(shù)分布，N(0,1) 標(biāo)準(zhǔn)正態(tài)

22、分布。其他常見的分布如卡方分布、F分布等的抽樣方法見表4.3.1。,四、其他分布隨機(jī)數(shù)的產(chǎn)生,1、 直接抽樣法 由基本定理我們知道，對于有些隨機(jī)變量可以建立與R的函數(shù)關(guān)系，因此只需對R進(jìn)行抽樣，利用函數(shù)的映射關(guān)系我們就可以方便地得到該隨機(jī)變量的抽樣了。如前面的指數(shù)分布隨機(jī)數(shù)。2、變換抽樣 產(chǎn)生隨機(jī)變量的變換抽樣方法，是討論均勻分布的不同函數(shù)分布，為隨機(jī)變量抽樣提供一些簡單可行的算法。在概率論中，從不同的角

23、度出發(fā)，對隨機(jī)變量函數(shù)進(jìn)行了討論，以下列出一些結(jié)果。,,設(shè)隨機(jī)變量 X 具有密度函數(shù) ，是對隨機(jī)變量 X 的變換，且 的逆函數(shù)存在，記為 有一階連續(xù)導(dǎo)數(shù)，則隨機(jī)變量的密度函數(shù),,,,，,,,,（4.4.1),例4.4.1 R，1-R均為（0，1）上的均勻分布隨機(jī)變量，設(shè)隨機(jī)向量（X，Y）具有二維聯(lián)合密度。對于隨機(jī)變量 X，Y

24、進(jìn)行函數(shù)變換,,其中，函數(shù) ， 的逆變換存在，記為,,,,且存在一階偏導(dǎo)數(shù)，設(shè) J 為 Jacobi 變換,,則隨機(jī)變量的二維聯(lián)合密度為,,(4.4.2),例4.4.2 用變換抽樣產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量 U 隨機(jī)變量 U 的密度函數(shù)為:,,取隨機(jī)數(shù),，則,,相互獨立、服從N（0，1）分布。解上面的兩個方程，得逆變換公式,,,由(4.4.2)，隨機(jī)變量 的密度函數(shù),,,從而我

25、們知道是兩個獨立的標(biāo)準(zhǔn)正態(tài)分布。下面還有幾個常用的二維函數(shù)的變換結(jié)果。 1）  隨機(jī)變量和 的密度函數(shù),,,（4.4.3),2) 隨機(jī)變量的積 的密度函數(shù),,(4.4.4),3) 隨機(jī)變量的商 X/Y,,(4.4.5),3、值序抽樣 值序統(tǒng)計量 ，通常稱為值序量。是隨機(jī)向量

26、 的一個函數(shù)，取其一個實現(xiàn)并排序得其中的第 l 個值為函數(shù)值。顯然這是一個統(tǒng)計量,,,,若隨機(jī)變量的各個分量獨立且同分布，則值序統(tǒng)計量 的密度函數(shù)和分布函數(shù)分別為：,,,,,這里，F(xiàn)(x), f(x) 分別為隨機(jī)變量在分布函數(shù)和密度函數(shù)。,特別當(dāng)隨機(jī)變量 為[0，1]上的均勻分布時，得密度函數(shù)

27、為：,,,這是 分布的密度函數(shù),,因此我們可以很容易地產(chǎn)生特殊的 分布的隨機(jī)數(shù)。,,,例4.4.3 產(chǎn)生服從 分布的隨機(jī)數(shù)若隨機(jī)變量 的密度函數(shù)為：,,,,,,,利用（0，1）均勻分布隨機(jī)數(shù)我們可以如下產(chǎn)生特殊 隨機(jī)數(shù),（1）,（2）,（3）,（4）,我們來驗證（2），即是否服從a=1，b=5的貝塔分布，按公式抽5個均勻分布的隨機(jī)數(shù)，取其中最小的為一個樣本，共取10

28、00個，然后用分布的Kolmogorov-Smirnov擬合優(yōu)度檢驗，程序如下：,% 例4.4.3驗證(2)% 抽5個均勻分布隨機(jī)數(shù)，產(chǎn)生一個服從bata分布的隨機(jī)數(shù)for i=1:1000 B(i)=min(rand(1,5));endB=B'; % 轉(zhuǎn)置為列向量h=kstest(B, [B betacdf(B,1,5)]) % 檢驗服從bata分布嗎？,計算結(jié)果為：h=0接受原假設(shè)，服從a=1，b=5的

29、貝塔分布。學(xué)習(xí)者可以驗證其余的抽樣公式。,五、正態(tài)分布隨機(jī)數(shù)的產(chǎn)生,正態(tài)分布在數(shù)理統(tǒng)計中具有基礎(chǔ)性的作用，因此產(chǎn)生高質(zhì)量的正態(tài)分布有重要的 意義，這一節(jié)我們將介紹幾種數(shù)值方法求正態(tài)分布，利用中心極限定理產(chǎn)生正態(tài)分布的隨機(jī)數(shù)。,1）    利用中心極限定理產(chǎn)生隨機(jī)數(shù) 中心極限定理：,,設(shè),服從均值為μ，方差為σ2 的某一分布，令,,（4.5.1),則當(dāng) n 充分大時， 漸近地服從標(biāo)準(zhǔn)

30、正態(tài)分布N（0，1）注意，這個定理沒有指出隨機(jī)變量x 是服從什么分布的。這正是該定理的神奇之處。我們現(xiàn)在已經(jīng)能產(chǎn)生[0，1]均勻分布的隨機(jī)數(shù)了，那么我們可以利用這個定理來產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)。,,現(xiàn)在我們產(chǎn)生n個[0，1]均勻分布隨機(jī)數(shù)，由公式（4.5.1）我們有：,,,為編程方便，我們特別選 n = 12，則,,這樣我們很方便地就把標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)計算出來了。,例4.5.1 利用中心極限定理產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)并檢

31、驗,% example 4.5.1clc,clearfor i=1:1000R=rand(1,12);X(i)=sum(R)-6;endX=X';m=mean(X)v=var(X)subplot(1,2,1),cdfplot(X)subplot(1,2,2),histfit(X)h=kstest(X, [X normcdf(X, 0,1)]),結(jié)果為：H=0, 接受原假設(shè)，變換后的確為標(biāo)準(zhǔn)正態(tài)分布

32、。,2) Hasiting 有理逼近法,這是一種計算速度快，也能滿足一定精度的算法。我們可以構(gòu)造分布函數(shù)反函數(shù)的近似逼近公式，來產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)。其計算公式為：,,（4.5.2),這里 ，系數(shù)為：,,a0 = 2.515517 b1 = 1.432788 a1 = 0.802853

33、 b2 = 0.189269 a2 = 0.010328 b3 = 0.001308,六、利用統(tǒng)計工具箱產(chǎn)生隨機(jī)數(shù),在MATLAB統(tǒng)計工具箱中為我們提供了大量的產(chǎn)生各種隨機(jī)數(shù)發(fā)生器程序，我們只需要調(diào)用就可以產(chǎn)生我們想要的隨機(jī)數(shù)。,隨機(jī)數(shù)發(fā)生器,例4.5.2,例4.6.1 設(shè)正態(tài)二維分布的密度函數(shù)為：,作二維分布的散點圖和直方圖,% 例4.6.1 二為獨立正態(tài)分布的散點圖和直方圖mu

34、 = [0 0];sigma = [1 0; 0 1];r = mvnrnd(mu,sigma,1000);subplot(1,2,1),plot(r(:,1),r(:,2),'+')subplot(1,2,2),hist3(r,[10 10]),七、隨機(jī)數(shù)的檢驗,我們已經(jīng)基本搞清偽隨機(jī)數(shù)的產(chǎn)生原理，由于并不是真正的隨機(jī)數(shù)，很自然的問題是，它們是否具有真正隨機(jī)數(shù)的那些統(tǒng)計性質(zhì)如參數(shù)大小、獨立性，均勻性等等。設(shè)：

35、隨機(jī)數(shù)具有連續(xù)的分布函數(shù)F（X），則隨機(jī)變量R=（X）是均勻分布（0，1）的隨機(jī)變量，因此如果R通過統(tǒng)計檢驗隨機(jī)變量 X 也可以通過。因此我們以下著重討論均勻分布R的檢驗問題，再簡單地討論正態(tài)隨機(jī)數(shù)檢驗問題。,統(tǒng)計推斷原理統(tǒng)計量的定義：設(shè) 為 隨機(jī)變量序列，則隨機(jī)序列的函數(shù)稱為統(tǒng)計量。記為：,顯然統(tǒng)計量 S 也是隨機(jī)變量。既然是隨機(jī)變量，它們就應(yīng)該有其分布或稱總體的規(guī)律，當(dāng)然也有各種

36、數(shù)字特征。例如均值、標(biāo)準(zhǔn)差、方差等等各階矩。我們的統(tǒng)計推斷方式是：（1）H0：某假定成立（2）在假定成立的條件下構(gòu)造統(tǒng)計量S,（3）統(tǒng)計量構(gòu)造完畢，我們也就知道了該統(tǒng)計量的全部統(tǒng)計 規(guī)律。如它的分布函數(shù)，或密度函數(shù)各階矩等。（4）根據(jù)統(tǒng)計量的分布，在給定的顯著性水平α， 對統(tǒng)計量 S 的一次抽樣確定以 1-α為概率的區(qū)域，該區(qū)域稱為接 受域 。如果該次抽樣計算出統(tǒng)計量 S 的

37、值 s 落入該領(lǐng)域， 我們就接受原假，否則推翻原假設(shè)。這個就是小概率事件在一次實驗實際不可能發(fā)生原理。落入由α 確定的區(qū)域是一個小概率事件，在一次實驗中我們認(rèn)為是不可能發(fā)生的。,,,,,,接受域,1、 統(tǒng)計檢驗中兩類常用統(tǒng)計量的構(gòu)造檢驗方法（1）設(shè)隨機(jī)變量 X 具有數(shù)學(xué)期望E（X）= ，和有限方差D（X）= ，我們抽 N 次樣本。X1，X2……，XN，當(dāng)N充分大時，統(tǒng)計量,,,

38、(4.6.1),以 N（0，1）為極限分布，取顯著水平 α = 0.05, 則拒絕為（—∞，—1.96），（1.96，∞）。當(dāng) ＞1.96，則認(rèn)為差異顯著，拒絕假設(shè) E（X）=μ,,（2）將樣本X1，X2，……，NN，按一定規(guī)則分為不相交的 K 組，記 i 組的觀測頻數(shù)為 ni（i=1,……，k），若隨機(jī)變量 X 落于弟 i 組的概率為 Pi，則得理論頻數(shù)m i= N × Pi，由ni，mi構(gòu)造統(tǒng)

39、計量。,,=,,(4.6.2),漸近服從自由度為 的分布，簡記 這里，l 是確定概率 P 時，由子樣 X 中給出的約束條件數(shù)。為有效地進(jìn)行統(tǒng)計檢驗，一般要求（4.6.2）的樣本數(shù) N＞30；在（4.6.2）中k≥5，mi≥5。,,當(dāng)統(tǒng)計量的自由度 時，U =漸近服從N（0，1）分布。&

40、#160;服從柯爾莫格洛夫—斯米爾諾夫的統(tǒng)計量,,,,取顯著性水平α= 0.05, 可拒絕原假設(shè)。,我們可以用這種方法進(jìn)行獨立性，及擬合優(yōu)度檢驗。,2、   隨機(jī)數(shù)的統(tǒng)計性質(zhì)檢驗,設(shè)均勻分布的隨機(jī)樣本為： （4.6.3）是一組需要檢驗的[0，1]均

41、勻分布隨機(jī)數(shù)，利用(4.6.1)，(4.6.2) 我們構(gòu)造統(tǒng)計量來對常用的參數(shù)進(jìn)行統(tǒng)計檢驗。,1）參數(shù)檢驗參數(shù)檢驗是檢驗參數(shù)估計值與理論值的差異是否顯著的方法，設(shè)我們有[0，1] 均勻分布樣本 R1，R2，……RN （4.6.4）我們可以構(gòu)造以下統(tǒng)計量（參數(shù)估計量）,,,,即隨機(jī)變量R的一階距，二階距和方差這些參數(shù)的估

42、計量。根據(jù)隨機(jī)變量R的理論分布，不難計算,,,,,,,現(xiàn)在我們可以應(yīng)用中心極限定理及（4.5.1）來構(gòu)造統(tǒng)計量了,,,,漸近服從N（0，1）分布，從而可以推斷樣本估計值與理論參數(shù)的差異。我們還可以利用同樣的方法去構(gòu)造三階矩、四階矩，對隨機(jī)變量的偏度和峰度進(jìn)行參數(shù)檢驗。,例4.7.1 對隨機(jī)數(shù)發(fā)生器unifrnd產(chǎn)生的1000個隨機(jī)數(shù)進(jìn)行均值、方差、偏度和峰度等的參數(shù)的檢驗。偏度計算方法為：,u3=mean(((R-R_mean)/R

43、_std).^3)*0.408248*sqrt(n);峰度計算方法為：uu=mean(((R-0.5)/sqrt(1/12)).^4)-1.75u4=uu*0.204124*sqrt(n);,function [s1,s2,s3,s4]=moment_test(R)% 對（0，1)均勻分布隨機(jī)數(shù)進(jìn)行矩檢驗n=length(R);R_mean=mean(R);R_var=var(R);R_std=std(R);u1=sqrt

44、(12*n)*(R_mean-0.5);if abs(u1)<1.96 s1='pass';else s1='*';end,% 對方差進(jìn)行檢驗var(R)u2=sqrt(180*n)*(R_var-1/12)if abs(u2)<1.96 s2='pass';else s2='*';end% 對偏度進(jìn)行檢驗u3

45、=mean(((R-R_mean)/R_std).^3)*0.408248*sqrt(n);if abs(u3)<1.96 s3='pass';else s3='*';end,% 對偏度進(jìn)行檢驗uu=mean(((R-0.5)/sqrt(1/12)).^4)-1.75u4=uu*0.204124*sqrt(n);if abs(u4)<1.96 s4='

46、;pass';else s4='*';end,調(diào)用該函數(shù)的主程序為：R=rand(1,1000);[s1,s2,s3,s4]=moment_test(R);S1 % 顯示對均值的推斷，pass為接受，*拒絕S2 % 顯示對方差的推斷，pass為接受，*拒絕S3 % 顯示對偏度的推斷，pass為接受，*拒絕S4

47、 % 顯示對偏度的推斷，pass為接受，*拒絕,計算結(jié)果為：s1 = passs2 = passs3 = passs4 = pass,2) 均勻性檢驗或稱分布的擬合優(yōu)度檢驗均勻性檢驗，又稱頻率檢驗，檢驗它的經(jīng)驗分布頻率和理論頻率的差異是否顯著。把[0，1]區(qū)間分為k個等區(qū)間，按 ri 取值的大小把（4.6.3）分為 k 組，設(shè)有ni個隨機(jī)數(shù)屬于i組，即共有個滿足不等式（i－1）/k＜r＜i/k。根據(jù)均勻

48、性假設(shè)，滿足在每個小區(qū)間的概率為Pi=1/k，mi=N/k由（4.6.2）得統(tǒng)計量,,,漸近服從 分布，而統(tǒng)計量,,(4.6.8),理論頻數(shù)ni,mi,,為累積頻率檢驗，漸近服從柯莫洛戈洛夫—斯米爾諾夫分布取顯著水平 α = 0.05，當(dāng)DN＞1.35 時拒絕均勻性假設(shè)。,（4.6.9),例4.7.2 對分布的擬合優(yōu)度進(jìn)行檢驗。（1）漸近服從柯莫洛戈洛夫—斯米爾諾夫檢驗,R=unifr

49、nd(0,1,1000,1);h=kstest(R, [R unifcdf(R, 0, 1)]),結(jié)果為：h=0， 接受原假設(shè)，是來自（0，1）均勻分布隨機(jī)數(shù),R=unifrnd(0,1,1,1000);% 構(gòu)造卡方統(tǒng)計量k=12;n=length(R);n1=hist(R,k); % 計算每個區(qū)間的頻數(shù)kf_7 = k/n*(sum((n1-n/k).^2)); % 計算分位點即統(tǒng)計量chi2_p=chi2cdf

50、(k-1,kf_7); % 計算下側(cè)概率 if chi2_p<0.95 chi2_str='pass';else chi2_str='*';endchi2_str,（2）利用卡方分布進(jìn)行檢驗,3）獨立性檢驗隨機(jī)數(shù)獨立性檢驗，其重點是檢驗（4.6.3）中前后各隨機(jī)數(shù)的統(tǒng)計相關(guān)是否異常。我們知道統(tǒng)計性能良好的隨機(jī)數(shù)數(shù)，前后之間應(yīng)該是統(tǒng)計獨立的，若存在明顯的相關(guān)性則

51、視這批隨機(jī)數(shù)為不合格。,這里： 表示i 階滯后相關(guān)系數(shù)，為0時表示沒有相關(guān)性。,（1）相關(guān)系數(shù)檢驗相關(guān)系數(shù)取值為零是兩個隨機(jī)變量獨立的必要條件，其值的大小給出它們之間線性相關(guān)強(qiáng)弱的測度，可以用來檢驗隨機(jī)數(shù)的獨立性。例如隨機(jī)序列與其自身持滯后序列的相關(guān)性。例如： 一階滯后相關(guān)性,,j 階滯后相關(guān)性,利用(4.6.3)數(shù)據(jù)構(gòu)造統(tǒng)計量,,（4.6.10),對充分大的N（如N—j>50），取零假設(shè)則統(tǒng)計量,

52、,（4.6.11),漸近服從N（0，1）分布。,例4.7.3對1~7階滯后相關(guān)系數(shù)進(jìn)行檢驗。,function [sacf1,sacf2,sacf3,sacf4,sacf5,sacf6,sacf7]=acf_test(R)% 獨立性的自相關(guān)AFC檢驗R_mean=mean(R);R_var=var(R);n=length(R);for i=1:7rou(i)=sum(((R(1:n-i).*R(i+1:n)-R_mean^2))

53、/R_var)*sqrt(1/(n-i));endrouif abs(rou(1))<1.96 sacf1='pass';else sacf1='*';endif abs(rou(2))<1.96 sacf2='pass';else sacf2='*';end,（2）聯(lián)列表檢驗在X—Y平面上，將單位正方形分為個相

54、等的小正方形，把隨機(jī)數(shù)（4.6.3）按出現(xiàn)的先后順序兩兩分組，如取,,記落入小正方形 ( i, j) 內(nèi)的觀測數(shù)為 nij ( i, j = 1, 2, ……, k)，令,,,H0：數(shù)據(jù)（4.6.3)來自均勻分布隨機(jī)數(shù),,則可以構(gòu)造統(tǒng)計量,其極限分布是服從自由度為（k-1)2 的卡方分布，記為,,(4.6.12),例4.7.4 對隨機(jī)數(shù)進(jìn)行聯(lián)列表檢驗。,R=unifrnd(0,1,1000,2);N

55、=2000k=6n=hist3(R,[k k]) % 產(chǎn)生每個小正方形落入的個數(shù)ni=sum(n');nj=sum(n);nij=ni'*nj;n_sum=sum(sum(n.^2./nij))-1chi2_2=N*n_sumchi2_p=chi2cdf((k-1)^2,chi2_2);if chi2_p<0.95 chi2_str='pass';else c

56、hi2_str='*';endchi2_str,計算結(jié)果：pass,4) 組合規(guī)律檢驗該檢驗用于對數(shù)據(jù)（4.6.3）的組合規(guī)律進(jìn)行檢驗，按隨機(jī)數(shù)先后出現(xiàn)的順序，根據(jù)一定的規(guī)則把這些隨機(jī)數(shù)組合起來，檢驗他們的組合規(guī)律的樣本性質(zhì)與理論性質(zhì)的差異，并進(jìn)行判斷。  （1） 撲克（顏色）檢驗將（4.6.3）的數(shù)據(jù)按順序分為8個隨機(jī)數(shù)一組，記組數(shù)為L，對每組的隨機(jī)數(shù)取其小數(shù)點后第一位，并以八進(jìn)制來表示。,

57、,例如：如下一組數(shù) 1 2 3 4 5 6 7 8 (0.001, 0.615, 0.516, 0.316, 0.815, 0.95, 0.11, 0.216)小數(shù)點后面第一位用 8 進(jìn)制表示，我們有 （ 0, 6, 5,

58、 3, 0, 1, 1, 2 ）一般地，我們們記為 顯然有,,,,當(dāng)中有 k 種相同數(shù)字時，稱該序列為 k 色向量，例如:（1，1，3，4，4，6，7，0） 為6色（5，5，5，5，5，5，5，5） 為1色（0，2，5，4，7，6，1，3） 為8色將向量按顏色的不同分為 8 類 。從理論上我們

59、有以下頻率。,,,,,記 為落入第 k 組向量的個數(shù)， 為理論頻數(shù)，則我們可以構(gòu)造統(tǒng)計量如下：,(4.6.13),例4.7.5 對1000個隨機(jī)數(shù)進(jìn)行撲克牌檢驗。解題思路：1）對1000個隨機(jī)數(shù)按順序每8個一組，共分125組。2）分別對每組元素小數(shù)點后第一位進(jìn)行模為8的運(yùn)算。3）計算每組的顏色數(shù)，并放入數(shù)組cc4）構(gòu)造自由度為4的卡方分布統(tǒng)

60、計量進(jìn)行檢驗。,% 撲克牌檢驗，2005.10.24clear,clcR=rand(1,1000);% poke testn=length(R);kk=8;jj=1;% 求125組數(shù)據(jù)的顏色數(shù)，并放入矩陣cc中。for ii=1:kk:n% 對每組數(shù)小數(shù)點后第一位取模為8運(yùn)算rr=10*R(ii:ii+7);pk=mod(fix(rr),8);,,% 計算每組的顏色數(shù)pk=sort(pk);j=1;for i=

61、1:7 if pk(i)~=pk(i+1); j=j+1; endendcc(jj)= j+1;jj=jj+1; % 將顏色數(shù)輸入到數(shù)組cc中end,% 構(gòu)造自由度為4的卡方分布，并進(jìn)行檢驗。nn(1)=sum(n_pk(1:3));nn(2:4)=n_pk(4:6);nn(5)=sum(n_pk(7:8));m=[0.02 0.17 0.4205 0.3195 0.0697]*n;chi

62、_4=sum((m-nn).^2./m); % 構(gòu)造統(tǒng)計量p=chi2cdf(4,chi_4);if p<0.95 str_pk='pass';else str_pk='*';endstr_pk,(2) 連（run）檢驗把隨機(jī)數(shù) (4.6.3)按一定規(guī)則進(jìn)行分類，如分為兩類，分別記為a, b, 得到形如 aa bbb aaaaaaaa bb

63、 ……由兩類元素a,b組成的序列，我們把位于異類元素之間的同類元素，如abbbba中的bbb類元素稱為一個連，連中包含同類元素的個數(shù)稱為連長，顯然連長是一個隨機(jī)變量。在隨機(jī)數(shù)序列中，出現(xiàn)連長為 i 的連數(shù)記為 ?？傔B數(shù)記為，構(gòu)成進(jìn)行檢驗的統(tǒng)計量。因此，隨機(jī)數(shù)的連檢驗是按照隨機(jī)數(shù)出現(xiàn)的先后順序，重點檢驗它的連貫現(xiàn)象是否異常的一種方法。,,,l      

64、;   正負(fù)連檢驗把隨機(jī)數(shù)序列 (4.5.3) 變換為 {ri － 1/2} 按正負(fù)分為兩類，這時a = -1，b= 1 組成正負(fù)兩類連，根據(jù)均勻性，獨立性假設(shè)，出現(xiàn) a, b 的概率都是0.5且有E（l）=N/2+1， D（l ）=（N-1）/ 4，P{1=k}=2-k（k=1，……），則可按（4.5.1），（4.5.2）構(gòu)造統(tǒng)計量U，和卡方統(tǒng)計量進(jìn)行檢驗。l   &#

65、160;     升降連檢驗把隨機(jī)序列（4.5.3）按生序或降序規(guī)律分為兩類，表示隨機(jī)數(shù)的增減及其長度的變化規(guī)律，組成升降兩類連，這時有E(l)=(2N—4)/3，D(1) = (16N-29)/90,從而我們又可以按(4.5.1)，(4.5.2)構(gòu)造 U，和卡方統(tǒng)計量進(jìn)行檢驗了。,,例4.7.6 對1000個隨機(jī)數(shù)進(jìn)行升降連檢驗。解題思路：1）對1000個隨機(jī)數(shù)按升降連產(chǎn)生連長隨機(jī)數(shù)。

66、2）分別構(gòu)造標(biāo)準(zhǔn)正態(tài)分布統(tǒng)計量和卡方統(tǒng)計量。3）對兩個統(tǒng)計量進(jìn)行統(tǒng)計檢驗。,function [striud_u,striud_chi2]=run_ud_test(R);% 升降連檢驗，包括正態(tài)和卡方檢驗m1=[0.625 0.275 0.079167 0.020833];n=length(R);R1=diff(R);AN=[0 0 0 0];% 搜索并計算總連長k=1;j=1;i=1;mk(1)=1;while i&

67、lt;=n-2 if R1(i)*R1(i+1)<0 mk(j+1)=i+1;j=j+1; end i=i+1; end,[mkmax,mki]=max(mk);mk=mk(1:mki);mk=diff(mk);nn=hist(mk,5);n1=[nn(1:3) sum(nn(4:5))];Run_tol=sum(nn);m1=m1*Run_tol;% 計算正態(tài)和卡方檢驗

68、檢驗u=(3*Run_tol-2*n+4)/sqrt(1.6*n-2.9);chi_3=sum((n1-m1).^2./m1);% 檢驗if abs(u)<1.96 striud_u='pass';else striud_u='*';end,5、檢驗隨機(jī)序列性質(zhì)的其它一些統(tǒng)計量除了上面介紹的方法外，事實上還有其他很多常用的統(tǒng)計量如，最小值MINMUN、即MINIMUN=

69、MIN{X1, X2, ……XN}。最大值MAXIMUN，即MAXIMUN= MAX{X1, X2, ……XN}。極差RANG=MAXIMUN - MINIMUN等。我們可以用上面的方法構(gòu)造統(tǒng)計量對它們進(jìn)行檢驗。,if chi2cdf(3,chi_3) < 0.95 striud_chi2='pass';else striud_chi2='*';end,l  

70、       模（MODE）這是一種中心趨勢的度量，類似與均值。如對一個隨機(jī)變量X抽了六個樣 6，10，10，4，4，10。則這個樣本的模是10，即出現(xiàn)次數(shù)最多的那個數(shù)，對一個大樣本來講，如果有N個數(shù)都出現(xiàn)同樣的最多數(shù)，則取其中值為最小的那個。l        

71、; 中位數(shù)（MEDIAN）這也是中心趨勢的一種度量，將樣本X1，X2，……XN從小到大排列。記為X（1），X（2），……，X（N），其最中間的那個數(shù)為中位數(shù)。,l         偏度SKEWNESS前面的統(tǒng)計量是反映序列的中心離差趨，而偏度則是衡量X的密度是否偏向一邊，即不對稱的一種度量，在正態(tài)情況下，模，均值，中位數(shù)應(yīng)近乎于重疊于一點，當(dāng)這三點不