計(jì)算機(jī)中產(chǎn)生隨機(jī)數(shù)列的方法

計(jì)算機(jī)上的數(shù)學(xué)方法生成隨機(jī)數(shù)列是目前的一般方法。它的特點(diǎn)是內(nèi)存少，速度快。用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)列是根據(jù)確定的算法推算出來的,嚴(yán)格說來并不是隨機(jī)的,因此一般稱用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)列為偽隨機(jī)數(shù)列.不過只要用數(shù)學(xué)公式產(chǎn)生出來的偽隨機(jī)數(shù)列通過統(tǒng)計(jì)檢驗(yàn)符合一些統(tǒng)計(jì)要求,如均勻性、抽樣的隨機(jī)性等,也就是說只要具有真正隨機(jī)數(shù)列的一些統(tǒng)計(jì)特征,就可以把偽隨機(jī)數(shù)列當(dāng)作真正的隨機(jī)數(shù)列使用.由于是用算法產(chǎn)生的,因而本質(zhì)上是決定性的,再加上計(jì)算機(jī)字長有限,所以無論用什么算法產(chǎn)生的數(shù)列,在統(tǒng)計(jì)特征上都不可能完全與從均勻分布中抽樣所得的子樣完全相同,因而只能要求盡可能地近似.這種在計(jì)算機(jī)上用算法得到的統(tǒng)計(jì)性質(zhì)上近似于上均勻分布的數(shù),一般就稱為偽隨機(jī)數(shù),以區(qū)別于真正從上均勻分布中抽樣所得到的隨機(jī)數(shù).為了快速產(chǎn)生統(tǒng)計(jì)特征良好的偽隨機(jī)數(shù),數(shù)學(xué)上一般采用遞推公式xn+1=f(xn,xn-1,?,xn-k).xn+1=f(xn,xn?1,?,xn?k).在給定一組初值x0,x-1,x-2,…,x-k后,即可逐步求出x1,x2,x3,….在此基礎(chǔ)上產(chǎn)生隨機(jī)數(shù).1隨機(jī)隨機(jī)數(shù)生成算法1.1基于積分法的偽隨機(jī)數(shù)列生成產(chǎn)生偽隨機(jī)數(shù)列最早的方法是平方取中法,即將一個(gè)2s位十進(jìn)制隨機(jī)數(shù)平方后得到的一個(gè)4s位數(shù),去頭截尾取中間2s位數(shù)作為一個(gè)新的隨機(jī)數(shù),重復(fù)上述過程便得到一偽隨機(jī)數(shù)列.平方取中法的遞推公式xn+1=[x2n/10s](mod102s),(1)產(chǎn)生偽隨機(jī)數(shù)列rn=xn/102s=xn10-2s.(2)平方取中法的優(yōu)點(diǎn)為在計(jì)算機(jī)上易于實(shí)現(xiàn),內(nèi)存占用少,但仍存在對小數(shù)目偏倚的現(xiàn)象,均勻性不好,數(shù)列的長度和周期難以確定,對初始數(shù)據(jù)的依賴很大.對平方取中法進(jìn)行改進(jìn),產(chǎn)生了乘積取中法,如果要產(chǎn)生具有10進(jìn)制2s位的偽隨機(jī)數(shù)列,任取兩個(gè)初始隨機(jī)數(shù)x0,x1,用遞推公式xn+2=[xn+1xn/10s](mod102s),(3)產(chǎn)生偽隨機(jī)數(shù)列rn=xn10-2s.(4)乘積取中法與平方取中法比較,從產(chǎn)生的偽隨機(jī)數(shù)列長度及均勻性方面都有改善,但是數(shù)列長度還是不夠,而且對初始值的依賴性很大.1.2偽隨機(jī)數(shù)列的產(chǎn)生電子計(jì)算機(jī)善于進(jìn)行移位等邏輯運(yùn)算,應(yīng)用機(jī)器的這個(gè)特點(diǎn)有一類產(chǎn)生偽隨機(jī)數(shù)列的方法,該方法稱為移位法.如果字長為32位的計(jì)算機(jī),取一初始值x0,將x0左移7位得x01,右移7位得x02.將x01和x02進(jìn)行指令相加得到x1,再對x1進(jìn)行上述運(yùn)算過程得到x2,這樣多次重復(fù)可得正整數(shù)序列xn,取rn=xn2-32為上均勻分布的偽隨機(jī)數(shù)列通項(xiàng).這一算法遞推公式為xn+1=xn27+xn2-7(mod232),產(chǎn)生偽隨機(jī)數(shù)列rn+1=xn+12-32.(5)移位法運(yùn)算速度快,但是對初始值的依賴性也很大,一般地初始值不能取得太小,選得不好會(huì)使偽隨機(jī)數(shù)列長度較短.1.3u3000隨機(jī)數(shù)列的求解目前產(chǎn)生偽隨機(jī)數(shù)列比較好的方法為同余法,其中有混合同余法、加同余法和乘同余法.采用線性遞推公式xn+1=λxn+c(modΜ),0≤xn+1＜Μ,產(chǎn)生偽隨機(jī)數(shù)列rn=xn/Μ.(6)其中λ,c,M以及初值x0都是正整數(shù).容易看出rn滿足0≤rn＜1.若c≠0且λ≠1,則稱這種方法為混合同余法.當(dāng)c=0時(shí),xn+1=λxn(modΜ),0≤xn+1＜Μ?rn=xn/Μ?(7)稱為乘同余法.當(dāng)λ=1,c≠0時(shí),xn+1=xn+c(modΜ),0≤xn+1＜Μ,rn=xn/Μ?(8)稱為加同余法.雖然加同余法產(chǎn)生的序列周期長,電子計(jì)算機(jī)實(shí)現(xiàn)也很方便,只要進(jìn)行加法及移位運(yùn)算即可完成,但從理論上看,所得到隨機(jī)數(shù)列的性質(zhì)一般不如乘同余法和混合同余法.若λ=c=1,則有xn+1=xn+1(modM).這樣構(gòu)成的序列雖然周期可以達(dá)到M,但不是隨機(jī)的.因此乘子λ,增量c的選取是非常重要的,下面給出λ和c的所有使得周期為最大的選取方法,以便從中挑選符合統(tǒng)計(jì)性質(zhì)的偽隨機(jī)數(shù)列.本文不加證明地引用有關(guān)周期的一些結(jié)果:一個(gè)混合同余數(shù)列,達(dá)到周期M的充要條件是:1)c與M互素;2)對每一個(gè)M的素因子p,λ-1為p的倍數(shù);3)若M是4的倍數(shù),那么λ-1是4的倍數(shù).由上面的結(jié)論易知,混合同余法達(dá)到最大周期的條件為:在2進(jìn)制計(jì)算機(jī)上,若M=2s,應(yīng)取λ=4q1+1,c=2a1+1,x0為任意非負(fù)整數(shù),其中q1,a1為正整數(shù).乘同余法達(dá)到最大周期的條件為:當(dāng)M=2s,若取x0與M互素,λ為M的本原元素,應(yīng)取x0=2a2+1,λ=8qq±3,a2,q2為正整數(shù),此時(shí)乘同余法的最大周期為2s-2.2數(shù)值計(jì)算模型的檢驗(yàn)Monte-Carlo方法的要點(diǎn)是:對要解決的數(shù)值計(jì)算問題,構(gòu)造適當(dāng)?shù)母怕誓Ｐ?使要得到的解正好重合于概率模型中隨機(jī)變量的概率分布或數(shù)字特征,其后在計(jì)算機(jī)上用偽隨機(jī)數(shù)列對隨機(jī)變量進(jìn)行模擬得到一個(gè)大子樣的觀測數(shù)據(jù),進(jìn)行統(tǒng)計(jì)整理以后,給出問題的一個(gè)近似估計(jì).因此,Monte-Carlo方法是雙重近似,一是將數(shù)值計(jì)算問題用概率模型作近似,二是在計(jì)算機(jī)上用偽隨機(jī)數(shù)作近似抽樣值進(jìn)行統(tǒng)計(jì)整理作出一些估計(jì).為了應(yīng)用Monte-Carlo方法對上述的3類(6種)生成偽隨機(jī)數(shù)的方法進(jìn)行比較,設(shè)計(jì)如下的數(shù)學(xué)實(shí)驗(yàn):設(shè)g(x)是上的連續(xù)函數(shù),且0≤g(x)≤1.考慮定積分s=∫10g(x)dx(9)的值.如圖1所示,在單位正方形內(nèi),曲線y=g(x)下面的陰影A的面積就是積分值s.如果在圖1所示的單位正方形內(nèi)均勻地投點(diǎn)(x,y),則該隨機(jī)點(diǎn)落入曲線y=g(x)下陰影的概率為p=Ρ{y≤g(x)}=∫10∫g(x)0dydx=∫10g(x)ds=s.(10)由此想法構(gòu)造計(jì)算定積分的投點(diǎn)模型.向正方形0≤x≤1,0≤y≤1內(nèi)均勻投點(diǎn)(ξi,ηi),ξi,ηi是相互獨(dú)立的均勻隨機(jī)數(shù)列,第i次試驗(yàn)成功,即(ξi,ηi)落入A中,也就是滿足ηi≤g(ξi),若每次成功的概率為p,進(jìn)行n次試驗(yàn)成功了k次,則由大數(shù)定理知limn→∞kn=p(a.s).即n充分大時(shí),k/n依概率收斂到p.由于p=s,因此常取s≈k/n.3投點(diǎn)算法中的化合物2.下面應(yīng)用上述Monte-Carlo方法對第一部分所述的3類(6種)生成偽隨機(jī)數(shù)方法進(jìn)行比較,其步驟是:首先用算法產(chǎn)生隨機(jī)數(shù),然后用所生成的隨機(jī)數(shù)進(jìn)行簡單的定積分計(jì)算.以計(jì)算積分Ι2=∫10x2dx,Ι3=∫10x3dx,Ι4=∫10x4dx,Ι9=∫10x9dx為例.顯然由直接積分得:Ι2=∫10x2dx=13=0.3333,Ι3=∫10x3dx=14=0.2500,Ι4=∫10x4dx=15=0.2000,Ι9=∫10x9dx=110=0.1000.現(xiàn)在用投點(diǎn)方法來計(jì)算,以k/n作為Ii的近似值,這里k是n次投點(diǎn)試驗(yàn)中,使不等式ξim≥ηm(i=2,3,4,9)成立的個(gè)數(shù).下面將用混合同余法、乘同余法、加同余法、平方取中法、乘積取中法、移位法這6種方法分別產(chǎn)生隨機(jī)數(shù)列,并用投點(diǎn)法來計(jì)算上述4個(gè)積分.以下程序在VC++6.0上編譯通過.程序運(yùn)行結(jié)果:分別設(shè)定程序中的n為2,3,4,9即為分別求Ι2=∫10x2dx,Ι3=∫10x3dx,Ι4=∫10x4dx,Ι9=∫10x9dx.1n.2，n.2試驗(yàn)次數(shù)N=1600000,不同隨機(jī)數(shù)產(chǎn)生算法的實(shí)驗(yàn)結(jié)果:2n.3試驗(yàn)次數(shù)N=1600000,不同隨機(jī)數(shù)產(chǎn)生算法的實(shí)驗(yàn)結(jié)果:3n.4試驗(yàn)次數(shù)N=1600000,不同隨機(jī)數(shù)產(chǎn)生算法的實(shí)驗(yàn)結(jié)果:4初始x0的結(jié)果試驗(yàn)次數(shù)N=1600000,不同隨機(jī)數(shù)產(chǎn)生算法的實(shí)驗(yàn)結(jié)果:由上面的幾組試驗(yàn)可以看出,3種同余法都有比較好的性質(zhì),算出的結(jié)果精度高.而取中法(平方取中法、乘積取中法)和移位法的計(jì)算結(jié)果明顯是錯(cuò)誤的.在程序運(yùn)行時(shí)通過變量跟蹤會(huì)發(fā)現(xiàn),取中法(平方取中法、乘積取中法等)很快退化為0或產(chǎn)生周期較短的循環(huán).使用平方取中法,當(dāng)2s=4:初值x0=9835時(shí),x27=0;初值x0=6406時(shí),x19=8100,x20=6100,x21=2100,x22=4100,x23=8100,這時(shí)進(jìn)入了一個(gè)循環(huán);初值x0=5829時(shí),x36=8100,由上面的結(jié)果可以看出此時(shí)也進(jìn)入了循環(huán);初值x0=1234時(shí),x56=0.由此可以看出,平方取中法有很大的缺陷.而采用改進(jìn)了的乘積取中法,當(dāng)2s=4:初值x0=6513,x1=3245時(shí),x287=0;初值x0=4158,x1=3023時(shí),x93=0.由此可以看出,乘積取中法對平方取中法的改進(jìn)還是有一些效果.而采用移位法也很快退化為0或產(chǎn)生周期較短的循環(huán),并且對初值的依賴性很大,產(chǎn)生的隨機(jī)數(shù)列長度不能滿足一般的實(shí)驗(yàn)要求.在上面的程序中,移位為7位,M=221:初值x0=797152時(shí),x1=1378387,x1