一元連續(xù)型隨機(jī)數(shù)的生成與檢驗(yàn)

課程設(shè)計(jì)（綜合實(shí)驗(yàn)）報(bào)告（2011--2011年度第1學(xué)期）名稱：數(shù)學(xué)建模課程設(shè)計(jì)題日：一元連續(xù)型隨機(jī)數(shù)生成及檢驗(yàn)院系：數(shù)理系班級(jí)：學(xué)號(hào)：學(xué)生姓名：指導(dǎo)教師：設(shè)計(jì)周數(shù)：2周成績(jī)：日期：2012年1月12日摘要：本文具體闡述了線性混合同余法生成隨機(jī)數(shù)的過程，并對(duì)生成的隨機(jī)數(shù)進(jìn)行了參數(shù)檢驗(yàn)、均勻性檢驗(yàn)和獨(dú)立性檢驗(yàn)，驗(yàn)證了線性混合同余法生成隨機(jī)數(shù)的正確性。關(guān)鍵字：隨機(jī)數(shù)均勻分布偽隨機(jī)數(shù)2檢驗(yàn)1、問題重述：（1）背景人類研究產(chǎn)生隨機(jī)數(shù)已有一段漫長(zhǎng)的歷史，最早的隨機(jī)數(shù)產(chǎn)生基本上是手工實(shí)現(xiàn)的，如抽簽、擲骰子、發(fā)紙牌以及從箱子里摸編號(hào)球等都屬此類；20世紀(jì)初，統(tǒng)計(jì)學(xué)家和賭博者對(duì)隨機(jī)數(shù)發(fā)生興趣，制造出了很多能更快產(chǎn)生隨機(jī)數(shù)的機(jī)械裝置；隨著計(jì)算機(jī)的廣泛應(yīng)用，隨機(jī)數(shù)的產(chǎn)生更加得到重視，40?50年代，人們轉(zhuǎn)向用數(shù)學(xué)方法來產(chǎn)生隨機(jī)數(shù)。（2）問題設(shè)隨機(jī)變量X（總體）服從某種隨機(jī)分布，對(duì)其進(jìn)行了n次獨(dú)立觀察，得到一組簡(jiǎn)單隨機(jī)樣本氣,x2，...,X.，滿足：1）X],X2,...,X相互獨(dú)立；2）每一個(gè)X],X2,...,X都與總體X同分布.利用某種方法得到一串?dāng)?shù)列r,r，…,r，在一定的統(tǒng)計(jì)意義下可作12n為隨機(jī)樣本X,X，…,X的一組樣本值，稱r,r,...,r一組具有與X相同12n12n分布的隨機(jī)數(shù).現(xiàn)在需要尋求一種簡(jiǎn)便、經(jīng)濟(jì)、可靠并能在計(jì)算機(jī)上實(shí)現(xiàn)的產(chǎn)生隨機(jī)數(shù)的方法.而最常用、最基礎(chǔ)的隨機(jī)數(shù)是在(0,1)區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)。在得到一系列隨機(jī)數(shù)后，我們又面臨如何檢驗(yàn)所得的隨機(jī)數(shù)是否滿足使用要求的問題。2、問題分析直接抽樣法時(shí)，利用均勻隨機(jī)數(shù)r通過給定分布函數(shù)的逆F一i(r)來得到給定分布隨機(jī)數(shù)的抽樣法，其主要依據(jù)為下列統(tǒng)計(jì)中的結(jié)果。定理：設(shè)F(x)為一連續(xù)的分布函數(shù)，則若X?F(x)，則R=F(X)?U(0,1)；若R?U(0,1)，則X=F-i(R)?F(x)；其中F-1(y)=inf{x:F(x)>】0<y<1)稱為F(x)的逆。據(jù)此定理，生成均勻隨機(jī)數(shù)r后，由r=F(x)所能解出的x=F-1(r)即為給定分布F(x)的隨機(jī)數(shù)x的直接抽樣公式。例：用直接抽樣法產(chǎn)生指數(shù)分布f(x)二人e-人x，(x>0)的隨機(jī)數(shù)，其中人〉0為參數(shù)。解：指數(shù)分布的分布函數(shù)F(x)=1-e&，(x>0)，則對(duì)均勻隨機(jī)數(shù)r，由r=1-e九解得x=-ln(1-r)。又因r和1-r均服從(0,1)上的均勻分布，故x=-llnr，即為指數(shù)分布的直接抽樣力公式，在產(chǎn)生均勻隨機(jī)數(shù)r后，由該公式即得指數(shù)分布的隨機(jī)數(shù)x。而我們所真正需要的隨機(jī)數(shù)應(yīng)該不是真正的隨機(jī)數(shù)，而是統(tǒng)計(jì)上滿足使用需求的隨機(jī)數(shù)，產(chǎn)生方便，使用方便，易于產(chǎn)生。并滿足一下特性：分布均勻，盡可能接近U(0,1)；統(tǒng)計(jì)上獨(dú)立；產(chǎn)生速度足夠快；可以重現(xiàn)；有足夠長(zhǎng)的周期；用內(nèi)存盡可能的少。嚴(yán)格意義上來說，現(xiàn)在所有的隨機(jī)數(shù)發(fā)生器都不是在概率論意義下的真正隨機(jī)數(shù)，而只能是稱為偽隨機(jī)數(shù)，因?yàn)闊o論哪一種隨機(jī)數(shù)發(fā)生器都是采用遞推算法一一即有一定規(guī)律可言，這已經(jīng)違反了概率論意義下隨機(jī)的含義。但是，如果算法選擇合適，則通過該算法得到的數(shù)據(jù)，可以在一定程度下滿足統(tǒng)計(jì)檢驗(yàn)的要求，具有較好的統(tǒng)計(jì)特性，也就可以用于模擬。由于“偽”隨機(jī)數(shù)應(yīng)該是隨機(jī)數(shù)，又具有均勻、獨(dú)立的特點(diǎn)，因此要進(jìn)行以下三個(gè)方面的檢驗(yàn)：參數(shù)檢驗(yàn)；均勻性檢驗(yàn)；獨(dú)立性檢驗(yàn)。3、模型假設(shè)把偽隨機(jī)數(shù)近似看成隨機(jī)數(shù)。4、模型建立與求解分析歷史上常用的隨機(jī)數(shù)的生成方法有平方取中法、乘法取中法、常數(shù)取中法和菲波那契法。這些方法雖然簡(jiǎn)單、百觀、易于實(shí)現(xiàn)，但同時(shí)也存在數(shù)據(jù)內(nèi)部特性難以控制，容易出現(xiàn)退化等問題。此處我們使用混合線性同余法。通式：X=[aX+c]modmU:=X，m其中:a表示乘子；c表示增量；X0表示初值；m表示模；通常4個(gè)參數(shù)都取整數(shù)。因?yàn)閰?shù)的選擇不同對(duì)產(chǎn)生的隨機(jī)序列和統(tǒng)計(jì)特性有極大的影響，所以問題的關(guān)鍵在于如何選擇參數(shù)以得到好的結(jié)果。顯然，根據(jù)通式X.廣[aX+c]modm，可得：0<X<m—1；由于X為整數(shù)，因此由此所得到的隨機(jī)數(shù)肯定是有限個(gè)數(shù);I通式：X.1=[aX+c]modm=(aX+c)—[(aX+c)/m]*mlI令N=axX.+c貝U：X.1=N—[N/m]xmX.的范圍0<X<m；U]=X"m<1；X.為整數(shù)，所以max(X.])=m—1，0<X.1<m—1而且，其數(shù)據(jù)周期的重復(fù)出現(xiàn);進(jìn)一步討論因?yàn)閤.取值有限，因此肯定會(huì)在m次迭代內(nèi)出現(xiàn)重復(fù)值；又因?yàn)閄.+1取值，只與X.相關(guān)，則一旦乂重復(fù)出現(xiàn)，后續(xù)X.+1都會(huì)完全重復(fù)。X.+1為整數(shù)，所以只能為｛0,1,2,...師-1｝中的一個(gè)，則*也為有限長(zhǎng)，為｛0，1/m，2/m，…，m-1/m｝之一。由此帶來的問題：可能會(huì)不均勻，即最大的周期長(zhǎng)度可能是m，但是并不一定。m越大，X.+1的選擇范圍越大，數(shù)字出現(xiàn)的頻率越均勻氣的可選擇范圍也越大。反之如果m取的小，則U+1很難為均勻分布。通常，取m盡可能的大，如：取m=2b其中b為計(jì)算機(jī)最大位長(zhǎng)，如果取b=31，m=2x109，即該隨機(jī)數(shù)發(fā)生器能夠產(chǎn)生21億以上個(gè)不重復(fù)隨機(jī)數(shù)。Um為周期式的重復(fù)，在每個(gè)m段中，計(jì)算公式是相同的。且X..肯定是不隨機(jī)的，因?yàn)榭梢愿鶕?jù)遞歸調(diào)用發(fā)現(xiàn)，只要m，a，C，X0一旦確定，則X是函數(shù)可得的，本質(zhì)上不存在任何隨機(jī)性。i5、模型檢驗(yàn)分析：由于在實(shí)際應(yīng)用中，我們并不真正需要“真”的隨機(jī)數(shù)，只要滿足“某些條件”的“偽”隨機(jī)數(shù)已經(jīng)足夠了。所以我們現(xiàn)在面臨的問題是如何驗(yàn)證由混合線性同余法所得的隨機(jī)數(shù)是否滿足使用要求。由于“偽”隨機(jī)數(shù)應(yīng)該是隨機(jī)數(shù)，又具有均勻、獨(dú)立的特點(diǎn)，因此要進(jìn)行以下三個(gè)方面的檢驗(yàn)：參數(shù)檢驗(yàn)；均勻性檢驗(yàn)；獨(dú)立性檢驗(yàn)。檢驗(yàn)過程：1、參數(shù)檢驗(yàn)對(duì)于隨機(jī)數(shù)u,u，…,u(N為數(shù)據(jù)長(zhǎng)度)，一些常見的進(jìn)行檢驗(yàn)的參12N數(shù)包括：均值的估計(jì)值：u=￡u/n；方差的估計(jì)值：S2=上Z(u-u)2；n一1均值估計(jì)值的期望：E面=1；2均值估計(jì)值的方差：見u)=上；12n方差估計(jì)值的期望：E(s2)二1；12方差估計(jì)值的方差：D(s2)^1；180n可構(gòu)造如下統(tǒng)計(jì)量：x廣^一^=、：頑(u一2)，X2=S2-ES°=\E(s2-￡)；注意到當(dāng)N足夠大時(shí)，X,X近似服從n(0,1)分布，對(duì)給定的顯著水平a(假定為0.05)，查表可得臨界值為1.96。當(dāng)1yI>1.96時(shí)，拒絕eu)=1和D(訪=1的假設(shè)；1212n當(dāng)|匕|>1.96時(shí)，拒絕E(s2)二￡和D(s2)=-8-假設(shè)，否則接受這兩個(gè)假設(shè)。這樣，我們就把參數(shù)檢驗(yàn)轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)分布的檢驗(yàn)。卜面用混合線性同余法生成了400個(gè)隨機(jī)數(shù)(具體隨機(jī)數(shù)見附錄1)，我們先假設(shè)這些數(shù)服從u(0,1)分布，則有：TOC\o"1-5"\h\zU=-￡u=0.4979，s2=-^￡(u-u)2=0.0875；nin-1^^且有：E(u)=，D(u)=，E(s2)=，D(s2)=；212n12180n我們知道服從u(0,1)分布的均值u=1，方差82=.1。212我們?nèi)〗y(tǒng)計(jì)量X,X，u—E(u)—廠】、X=—_=%.12n(u-—)，1項(xiàng)2X2=S2D；)=、E(s2-￡)；當(dāng)n充分大時(shí)，則應(yīng)有X,X都近似服從n(0,1)。取以=0.05，查止態(tài)分布表Z=1.96；2由于|X]|=0.1441<1.96，叩=1.1117<1.96；所以在顯著性水平a=0.05條件下，可以接受總體均值u=1，方差82=1的假設(shè)。212另外，用matlab擬合出了正態(tài)分布函數(shù)及正態(tài)隨機(jī)數(shù)的經(jīng)驗(yàn)分布函數(shù)和u(0,1)分布函數(shù)及u(0,1)的經(jīng)驗(yàn)分布函數(shù)，分別取產(chǎn)生的隨機(jī)數(shù)個(gè)數(shù)n=20,50,200,2000正態(tài)分布函數(shù)及正態(tài)隨機(jī)數(shù)的經(jīng)驗(yàn)分布函數(shù):n=20n=50n=2000u(0,1)分布函數(shù)及u(0,1)的經(jīng)驗(yàn)分布函數(shù):n=20EmpiricalCDF°-91—°-8-----…-…十——°-7-—:?！灰?--一一--■-X-一一....i-.....一一........；0-3"-一一…/…一十一一…一一…一一……°；i01-tX---…--…十一一…一十…一。iin=2000'：...'2、均勻性檢驗(yàn)對(duì)于給定的隨機(jī)數(shù)，其均勻性檢驗(yàn)主要是采用頻率檢驗(yàn)法，手段是檢驗(yàn)所產(chǎn)生的隨機(jī)數(shù)落在各子區(qū)間的頻率和理論頻率之間的差異是否顯著。常用的頻率檢驗(yàn)法有兩種：K-S檢驗(yàn)(Kolmogorov-Smirnov)和/2檢驗(yàn)；在這里，我們使用Z2檢驗(yàn)。思路:將(0,1)區(qū)間分成k個(gè)相等的子區(qū)間，檢查總數(shù)為N的隨機(jī)數(shù)落在各子區(qū)間的實(shí)際頻數(shù)n與理論頻數(shù)的差異是否顯著。設(shè)落在第/個(gè)子區(qū)間的隨機(jī)數(shù)的實(shí)際頻數(shù)為n，落在第/個(gè)子區(qū)間的理論頻數(shù)為m，m=N/k；建立x2統(tǒng)計(jì)量:木2=8上冬=kE(n—N)2i=1'i=1當(dāng)樣本量足夠大時(shí)，該統(tǒng)計(jì)量近似服從自由度為k-1的X2分布，顯著水平a一般取0.01?0.05。x2當(dāng)統(tǒng)計(jì)量>臨界值(k—1),(1—a)時(shí)，則在該顯著水平上拒絕該批樣本為均勻的假設(shè)。X2檢驗(yàn)步驟：分區(qū)間k；統(tǒng)計(jì)n和計(jì)算m；計(jì)算統(tǒng)計(jì)量x2；根據(jù)顯著水平及自由度k-1從x2表查臨界值；比較x2和臨界值，確定是否接受。下面用x2擬合優(yōu)度檢驗(yàn)法對(duì)上面得到的隨機(jī)數(shù)做分布均勻性檢驗(yàn)，我們將樣本的取值范圍分布在10個(gè)等寬區(qū)間，即分成0?0.1，0,1?0.2，…，0.9?1，10個(gè)區(qū)間。統(tǒng)計(jì)實(shí)際落在每個(gè)區(qū)間內(nèi)的樣本個(gè)數(shù)n.如下：i12345678910ni32603040304052343840如對(duì)均勻分布，則這些數(shù)落在第"區(qū)間內(nèi)的概率1。按X2擬合法，取統(tǒng)計(jì)量V=—￡〃2/p-n=—￡n2-n=13.240i取以=0.05，查自由度k-1=10-1=9的X2分布表X2059=16.92>|V|=13.2故應(yīng)接受顯著水平為的a=0.05均勻性假設(shè)。3、獨(dú)立性檢驗(yàn)檢查隨機(jī)數(shù)樣本間是否存在相關(guān)性，常用的有下列三種方法:(1)相關(guān)性檢驗(yàn)’Zuu—u2'ii+k"N—k)設(shè)給定隨機(jī)數(shù)列U1,u2，..,un，前后相距為k的樣本的相關(guān)系數(shù)的估計(jì)值為pk，則若k階自相關(guān)系數(shù)為0，即p廣0。建立統(tǒng)計(jì)量：Vk=pkJN-k當(dāng)N充分大時(shí)，該統(tǒng)計(jì)量近似服從N(0,1)分布；給定顯著水平a=0.05,計(jì)算出k=1,2,…,k...時(shí)各階自相關(guān)系數(shù)的估計(jì)值；若I匕I>1.96，則拒絕相關(guān)系數(shù)pk=0的假設(shè)。（2）組合規(guī)律性檢驗(yàn)-Poker檢驗(yàn)用于檢查隨機(jī)數(shù)樣本中重復(fù)數(shù)字出現(xiàn)的頻率，從而檢查該組樣本的組合規(guī)律性。如，三位數(shù)：三個(gè)數(shù)字相同出現(xiàn)的概率為：P三個(gè)數(shù)字相同）=0.10.1=0.01三個(gè)數(shù)字不同出現(xiàn)的概率為：P三個(gè)數(shù)字不同）=0.90.8=0.72有兩個(gè)數(shù)字同時(shí)出現(xiàn)的概率為：p兩個(gè)數(shù)字相同）二C（3,2）（0.0.9）=1-0.01-0.72=0.27以此作為理論頻率與實(shí)際頻數(shù)作比較，然后對(duì)其進(jìn)行2檢驗(yàn)。在給定的顯著水平下，找出臨界值，根據(jù)臨界值確定差距是否顯著。（3）連貫性檢驗(yàn)根據(jù)給定樣本中是否存在每隔一些數(shù)字就連續(xù)出現(xiàn)增大或減少的趨勢(shì)，或連續(xù)在均值之上或在均值之下的現(xiàn)象來判斷給定的樣本是否相互獨(dú)立?？梢圆捎靡环N符號(hào)來表示一種事件的出現(xiàn)，連續(xù)出現(xiàn)的相同（或單個(gè)）事件的符號(hào)連就稱為一個(gè)“連”，這樣容易看到相同事件出現(xiàn)的規(guī)律性。如兩種類型的10個(gè)事件的符號(hào)+--+++-就構(gòu)成了6個(gè)連（相同狀態(tài)的事件構(gòu)成一個(gè)連）。常用的連貫性檢驗(yàn)方法有：升降連檢驗(yàn)和正負(fù)連檢驗(yàn)。在此，我們使用正負(fù)連檢驗(yàn)。判定給定樣本中是否每隔一些數(shù)字就出現(xiàn)連續(xù)在均值之上或在均值之下的趨勢(shì)。大于均值的樣本用+號(hào)表示，其總數(shù)為n「小于均值的樣本用-號(hào)表示，其總數(shù)為〃;2用這些+、-號(hào)組成的符號(hào)串中最大可能的連數(shù)為N=氣+七(N為樣本個(gè)數(shù))，最小連數(shù)為1；設(shè)》為一個(gè)樣本序列的連的總數(shù)(連續(xù)發(fā)生的+或-事件構(gòu)成一個(gè)連)則其均值和方差分別為：2nn1廣N22nn(2nn—N)bbN2(N—1)當(dāng)數(shù)字總數(shù)N>20時(shí)，連的總數(shù)a近似服從正態(tài)分布N(u,a2)；建立統(tǒng)計(jì)量：7_b—(2nnJN)—1/2o2nn(2nn—N)1N，(N—y1/2據(jù)顯著性水平確定臨界值(正態(tài))，以確定是否顯著?，F(xiàn)在我們將上面隨機(jī)數(shù)分別與均值u_0.5相減，取其符號(hào)得到一個(gè)正負(fù)號(hào)序列，統(tǒng)計(jì)其中的正號(hào)個(gè)數(shù)n_208，負(fù)號(hào)個(gè)數(shù)n=192，b=188；假設(shè)這些數(shù)式獨(dú)立分布的，取統(tǒng)計(jì)量b則應(yīng)有：u=丑2+-=200.18bN2^2=2nin2(2nin2—N)=99.4296bN2(N—1)b近似于正態(tài)分布。由于Z0=；「(2罟：N—N；/2=(188—200.18)/、99.4296=—1.2215[1212“1/2N2(N—1)取a=0.05，查正態(tài)分布表Z=1.96>|ZJ=1.22152故應(yīng)該接受這一假設(shè)。6、模型評(píng)價(jià)通過對(duì)線性混合同余法產(chǎn)生的隨機(jī)數(shù)進(jìn)行參數(shù)檢驗(yàn)、均勻性檢驗(yàn)和獨(dú)立性檢驗(yàn)，驗(yàn)證了線性混合同余法的優(yōu)越性，這種方法也不失為一種好的隨機(jī)數(shù)生成法。7、參考文獻(xiàn)鄧集賢等.概率論及數(shù)理統(tǒng)計(jì)(上、下冊(cè)).北京：高等教育出版社.第4版.2009.石博強(qiáng)等.MATLAB數(shù)學(xué)計(jì)算范例教程.北京：中國(guó)鐵道出版社.2004.高祖新等.統(tǒng)計(jì)計(jì)算與軟件實(shí)用統(tǒng)計(jì)計(jì)算.南京：南京大學(xué)出版社.1996.附錄：附錄1產(chǎn)生400個(gè)隨機(jī)數(shù)：u=rand(1,400)m=mean(u)v=var(u)x1=sqrt(12.*400).*(m-1./2)x2=sqrt(180.*400).*(v-1./12)0.01490.28820.81670.98550.01740.81940.62110.56020.24400.82200.26320.75360.65960.21410.60210.60490.65950.18340.63650.17030.53960.62340.68590.67730.87680.01290.31040.77910.30730.92670.67870.07430.07070.01190.22720.51630.45820.70320.58250.50920.07430.19320.37960.2764

0.77090.31390.63820.98660.50290.94770.82800.91760.11310.81210.90830.15640.12210.76270.72180.65160.75400.66320.88350.27220.41940.21300.03560.08120.85060.34020.46620.91380.22860.86200.65660.89120.48810.99260.37330.53140.18130.50190.42220.66040.67370.95730.19190.11120.56510.96920.02370.87020.02690.51950.19230.71570.25070.93390.13720.52160.89520.94240.33510.43740.47120.14930.13590.53250.72580.39870.35840.28530.86860.62640.24120.97810.64050.22980.68130.66580.13470.02250.26220.11650.06930.85290.18030.03240.73390.53650.27600.36850.01290.88920.86600.25420.56950.15930.59440.33110.65860.86360.56760.98050.79180.15260.83300.19190.63900.66900.77210.37980.44160.48310.60810.17600.00200.79020.51360.21320.10340.15730.40750.40780.05270.94180.15000.38440.31110.16850.89660.32270.73400.41090.39980.50550.16930.52470.64120.01620.83690.80350.69780.46190.08260.82070.19300.44540.01300.30870.87540.83530.33310.88070.47970.56080.61590.66190.61660.68510.51020.71400.51520.60590.96670.82210.31780.58770.13020.25440.80300.66780.01360.56160.45460.90490.28220.06500.47660.98370.92230.56120.65230.77270.10620.00110.54180.00690.45130.19570.78710.61860.01550.89090.76170.90700.75860.38070.33110.50410.56460.76720.77990.48410.80220.47100.20280.57960.66650.67680.94250.77010.73740.86630.99090.50390.62910.79260.44860.52440.17150.13070.21880.10550.14140.45700.78810.28110.22480.90890.00730.58870.54210.65350.31340.23120.41610.29880.67240.93830.34310.56300.11890.16900.27890.55680.48560.95220.23190.47870.52650.79270.19300.90960.92220.01330.76750.94730.81330.92380.19900.67430.92710.34380.594