隨機數(shù)產(chǎn)生原理及應(yīng)用
EmilMatthew(EmilMatthew@126.com)
摘要:
本文簡述了隨機數(shù)的產(chǎn)生原理，并用C語言實現(xiàn)了迭代取中法，乘同余法等隨機數(shù)產(chǎn)生方法，同時，還給出了在符合某種概率分布的隨機變量的產(chǎn)生方法。
關(guān)鍵詞: 偽隨機數(shù)產(chǎn)生，概率分布
1前言:
在用計算機編制程序時，經(jīng)常需要用到隨機數(shù)，尤其在仿真等領(lǐng)域，更對隨機數(shù)的產(chǎn)生提出了較高的要求，僅僅使用C語言類庫中的隨機函數(shù)已難以勝任相應(yīng)的工作。本文簡單的介紹隨機數(shù)產(chǎn)生的原理及符合某種分布下的隨機變量的產(chǎn)生，并用C語言加以了實現(xiàn)。當(dāng)然，在這里用計算機基于數(shù)學(xué)原理生成的隨機數(shù)都是偽隨機數(shù)。
注:這里生成的隨機數(shù)所處的分布為0-1區(qū)間上的均勻分布。我們需要的隨機數(shù)序列應(yīng)具有非退化性，周期長，相關(guān)系數(shù)小等優(yōu)點。
2.1迭代取中法：
這里在迭代取中法中介紹平方取中法,其迭代式如下:
Xn+1=(Xn^2/10^s)(mod 10^2s)
Rn+1=Xn+1/10^2s
其中，Xn+1是迭代算子，而Rn+1則是每次需要產(chǎn)生的隨機數(shù) 。
第一個式子表示的是將Xn平方后右移s位，并截右端的2s位。
而第二個式子則是將截尾后的數(shù)字再壓縮2s倍，顯然:0=<Rn+1<=1.
這樣的式子的構(gòu)造需要深厚的數(shù)學(xué)(代數(shù)，統(tǒng)計學(xué)，信息學(xué))功底，這里只是拿來用一下而已，就讓我們站在大師的肩膀上前行吧。
迭代取中法有一個不良的性就是它比較容易退化成0.
平方取中法的實現(xiàn):
#include <stdio.h>
#include <math.h>
#define S 2
float Xn=12345;//Seed & Iter
float Rn;//Return Val
void InitSeed(float inX0)
{
Xn=inX0;
}
/*
Xn+1=(Xn^2/10^s)(mod 10^2s)
Rn+1=Xn+1/10^2s
*/
float MyRnd()
{
Xn=(int)fmod((Xn*Xn/pow(10,S)),pow(10,2*S));//here can‘s use %
Rn=Xn/pow(10,2*S);
return Rn;
}
/*測試主程序，注意，這里只列舉一次測試主程序，以下不再重復(fù)*/
int main()
{
int i;
FILE * debugFile;
if((debugFile=fopen("outputData.txt","w"))==NULL)
{
fprintf(stderr,"open file error!");
return -1;
}
printf("\n");
for(i=0;i<100;i++)
{
tempRnd=MyRnd();
fprintf(stdout,"%f ",tempRnd);
fprintf(debugFile,"%f ",tempRnd);
}
getchar();
return 0;
}
前一百個測試生成的隨機數(shù)序列:
0.399000 0.920100 0.658400 0.349000 0.180100 0.243600 0.934000 0.235600 0.550700 0.327000 0.692900 0.011000 0.012100 0.014600 0.021300 0.045300 0.205200 0.210700 0.439400 0.307200 0.437100 0.105600 0.115100 0.324800 0.549500 0.195000 0.802500 0.400600 0.048000 0.230400 0.308400 0.511000 0.112100 0.256600 0.584300 0.140600 0.976800 0.413800 0.123000 0.512900 0.306600 0.400300 0.024000 0.057600 0.331700 0.002400 0.000500 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
容易看出其易退化成0的缺點.
2.2乘同余法：
乘同余法的迭代式如下:
Xn+1=Lamda*Xn(mod M)
Rn+1=Xn/M
各參數(shù)意義及各步的作用可參2.1
當(dāng)然，這里的參數(shù)的選取是有一定的理論基礎(chǔ)的，否則所產(chǎn)生的隨機數(shù)的周期將較小，相關(guān)性會較大。
經(jīng)過前人檢驗的兩組性能較好的素數(shù)取模乘同余法迭代式的系數(shù)為:
1） lamda=5^5,M=2^35-31
2） lamda=7^5,M=2^31-1
相應(yīng)C程序關(guān)鍵代碼段:
double long M;//請注意，這里一定要用到double long,否則計算2^32會溢出
float MyRnd()
{
Xn=fmod(Lamda*Xn,M);//here can‘s use %
Rn=Xn/M;
return Rn;
}
初始化段,應(yīng)有:
M=pow(2,35)-31;
圖1： 乘同余法生成的300隨機數(shù)的產(chǎn)生序列圖
圖2: 乘同余法生成的300隨機數(shù)的分布情況
可以看到，該隨機數(shù)生成方法所生成的隨機序列比較符合0-1上的均勻分布,不過在某些數(shù)據(jù)段還有些起伏。
2.3混合同余法：
混合同余法是加同余法和乘同余法的混合形式,其迭代式如下:
Xn+1=(Lamda*Xn+Miu)%M
Rn+1=Xn/M
經(jīng)前人研究表明，在M=2^q的條件下，參數(shù)lamda,miu,X0按如下選取，周期較大，概率統(tǒng)計特性好:
Lamda=2^c+1,c取q/2附近的數(shù)
Miu=(1/2+sqrt(3))/M
X0為任意非負整數(shù)
相應(yīng)C程序關(guān)鍵代碼段:
//init proper argu number
M=pow(2,32);
Lamda=pow(2,16)+1;
Miu=(0.5+sqrt(3)/6)/M;
float MyRnd()
{
Xn=fmod(Lamda*Xn+Miu,M);
Rn=Xn/M;
return Rn;
}
圖3: 乘同余法生成的300隨機數(shù)的產(chǎn)生序列圖
圖4: 乘同余法生成的300隨機數(shù)的分布情況
由圖4可以看出，該種隨機數(shù)生成方法已相當(dāng)接近0-1上的均勻分布。但在圖3中可以看出它的一個致命的弱點，那就是隨機數(shù)的生成在某一周期內(nèi)成線性增長的趨勢，顯然，在大多數(shù)場合，這種極富“規(guī)律”型的隨機數(shù)是不應(yīng)當(dāng)使用的。
下面的概率分布型隨機變量的生成，均采用乘同余法或C函數(shù)庫中的隨機數(shù)來生成0-1區(qū)間上的隨機數(shù)。
下面將C語言中的隨機數(shù)生成序列圖和Matlab中的隨機數(shù)生成序列圖列于下面，以作對比之用:
C語言生成的300個隨機數(shù)的序列圖
Matlab中rand函數(shù)生成的300個隨機數(shù)的序列圖
可以看出:乘同余法生成的隨機數(shù)序列的隨機性與上述兩個標(biāo)準(zhǔn)庫函數(shù)相接近。
3連續(xù)型隨機變量的生成:
3.1反函數(shù)法
采用概率積分變換原理,對于隨機變量X的分布函數(shù)F(X)可以求其反函數(shù)，得:
Xi=G(Ri)
其中,Ri為一個0-1區(qū)間內(nèi)的均勻分布的隨機變量.
F(X)較簡單時，求解較易，當(dāng)F(X)較復(fù)雜時，需要用到較為復(fù)雜的變換技巧。
例:已知炮彈對目標(biāo)的方位角Fi在0-2*P內(nèi)均勻分布，試用（0，1）均勻隨機數(shù)變換，模擬彈著點方位角的抽樣值Fi.
解: R=F(Fi)=Fi/2*PI
得Fi=G(R)=2*PI*R ,其中，R為0-1區(qū)間上的均勻分布的隨機數(shù).
程序略
試驗結(jié)果:
圖5:用反函數(shù)法生成的300隨機數(shù)的平均分布情況
由于這里相當(dāng)對0-1上的分布進行線性變換，所以變換后仍呈均勻分布是顯然的。
指數(shù)分布的分布函數(shù)為:
x<0時,F(x)=0 ; x>=0,F(x)=1-exp(-lamda*x)
利用反函數(shù)法，可以求得: x=-lnR/lamda
試驗結(jié)果:
圖6:用反函數(shù)法生成的300隨機數(shù)的指數(shù)分布情況
可以看出，生成的隨機量較好的符合了指數(shù)分布特征。
3.2正態(tài)分布隨機變量的生成:
正態(tài)分布在概率統(tǒng)計的理論及應(yīng)用中占有重要地位，因此，能產(chǎn)生符合正態(tài)分布的隨機變量就在模擬一類的工作中占有相當(dāng)重要的地位。下面介紹兩種方法。
這種方法便捷而有效，且具有一定的代表性，其基本思路是:
在概率密度的函數(shù)圖像的外圍畫一個大框，然后在這個框內(nèi)部產(chǎn)生隨機點（rx,ry），根據(jù)是否落在概率密度函數(shù)的下方，來決定是否要留下這個點。
經(jīng)過一定的計算變行，符合二維的正態(tài)分布的隨機變量的生成可按下面的方法進行：
1)產(chǎn)生位于0-1區(qū)間上的兩個隨機數(shù)r1和r2.
2)計算u=2*r1-1,v=2*r2-1及w=u^2+v^2
3)若w>1，則返回1)
4) x=u[(-lnw)/w]^(1/2)
y=v[(-lnw)/w]^(1/2)
如果為(miu,sigma^2)正態(tài)分布,則按上述方法產(chǎn)生x后，x’=miu+sigma*x
由于采用基于乘同余法生成的0-1 上的隨機數(shù)的正態(tài)分布隨機數(shù)始終無法能過正態(tài)分布總體均值的假設(shè)檢驗。而采用C語言的庫函數(shù)中的隨機數(shù)生成函數(shù)rand()來產(chǎn)生0-1 上的隨機數(shù)，效果較為理想。
關(guān)鍵程序段(funNorm返回一維的正態(tài)分布，而funNorm2則生成二維的隨機分布):
float funNorm(float miu,float sigma)
{
float r1,r2;
float u,v,w;
float x,y;
do
{
r1=MyRnd();
r2=MyRnd();
u=2*r1-1;
v=2*r2-1;
w=u*u+v*v;
}while(w>1);
x=u*sqrt(((-log(w))/w));
y=v*sqrt(((-log(w))/w));
return miu+sigma*x; //also could return miu+sigma*y;
}
typedef struct
{
float x;
float y;
}sPoint;
sPoint funNorm2(float miu1,float sigma1,float miu2,float sigma2)
{
float r1,r2;
float u,v,w;
float x,y;
sPoint mPoint;
do
{
r1=rand()/(float)32767;
r2=rand()/(float)32767;
u=2*r1-1;
v=2*r2-1;
w=u*u+v*v;
}while(w>1);
x=u*sqrt(((-log(w))/w));
y=v*sqrt(((-log(w))/w));
mPoint.x=miu1+sigma1*x;
mPoint.y=miu2+sigma2*x;
return mPoint;
}
列出在Matlab下對某次試驗（生成分布為N(0,1)的隨機數(shù)）的檢測結(jié)果:
[muhat,sigmahat,muci,sigmaci]=normfit(a)
[h,sig,ci]=ztest(a,0,1)
Output:
muhat =-0.0246 %顯著性為0.95的條件下，均值的點估計。
h = 0 %接受方差為1，均值估計為0的假設(shè)檢驗,不能拒絕假設(shè)。
sig =0.6700 %假設(shè)的成功概率.
由此可見，在這種條件下生成的正態(tài)隨機數(shù)序列基本能符合使用的要求,在精度上仍有該進的余地。
根據(jù)獨立同分布的中心極限定理，有:
這里，其實只要取n=12（這里，亦即生成12個0-1上的隨機數(shù)序列）就會有比較好的效果。
經(jīng)驗證，用該種方法生成生的隨機數(shù)序列同樣能比較好的符合正態(tài)分布特性。
由于生成的都是標(biāo)準(zhǔn)正態(tài)分布，所以，當(dāng)需要生成N(a,b)的正態(tài)分布隨機量時，根據(jù)正態(tài)分布的線性變換特性，只要用x=a*x0+b即可。（其中，x0表示生成的符合N(0,1)分布的正態(tài)隨機變量。方法3.1亦是如此）
4離散型隨機變量的生成:
離散型隨機變量的生成主要是希望得到在已知X符合某類分布的條件下，生成相應(yīng)的X。
4.1符合泊松分布的隨機變量的生成:
這里，采用”上限攔截”法進行測試，其本的思想是這樣的：
1）在泊松分布中，求出X取何值時，p(X=k)取最大值時，設(shè)為Pxmax.
其時，這樣當(dāng)于求解f(x)=lamda^k/k! 在k取何值時有最大值，雖然，這道題有一定的難度，但在程序中可以能過一個循環(huán)來得到取得f(x)取最大值時的整數(shù)自變量Xmax。
2) 通過迭代，不斷生成0-1區(qū)間上的隨機數(shù)，當(dāng)隨機數(shù)<Pxmax時，則終止迭代，否則重復(fù)(2)
3) 記錄迭代過程的次數(shù)，即為所需要得到的符何泊松分布的隨機量。
顯然，這種方法較為粗糙，在試驗的過程中發(fā)現(xiàn)：生成的的隨機量只能算是近似的服從
泊松分布，所以，更為有效的算法還有待嘗試。
4.2符合二項分布的隨機變量的生成:
符合二項分布的隨機變量產(chǎn)生類似上限攔截法，不過效果要好許多，這是由二項分布的特點決定的。
具體方法如下:
設(shè)二項分布B(p,n)，其中，p為每個單獨事件發(fā)生的概率：
關(guān)鍵算法:
i=0;reTimes=0
while(i<n)
{
temp=MyRnd();//生成0-1區(qū)間上的隨機變量
if(temp>1-p)
reTimes++;
i++;
}
顯然，直觀的看來，這種算法將每個獨立的事件當(dāng)作一個0-1分布來做，生成的0-1區(qū)間上的隨機數(shù)，若小于1-p則不發(fā)生，否則認為發(fā)生，這樣的生成方式較為合理。實驗結(jié)果也驗證了其合理性。
5實驗結(jié)論:
通過本實驗，使我感受到，要生成符合要求的隨機數(shù)序列，絕不是一件很輕松的事，除了要有相當(dāng)?shù)闹R儲備外，更應(yīng)當(dāng)有嚴謹求實的態(tài)度在其中；否則，光憑主觀感覺說某些隨機數(shù)的隨機性好，是會在實際應(yīng)用中是要栽跟頭的。
有些隨機數(shù)的生成仍未達到要求，有機會仍應(yīng)繼續(xù)被充，加深該方面理論和應(yīng)用的理解。
參考文獻:
[1]王可定,計算機摸擬及其應(yīng)用,東南大學(xué)出版社,1997
完成日:
附錄:
1全文測試程序下載:
http://emilmatthew.51.net/EmilPapers/06_14RandomNum/code.rar
若直接點擊無法下載(或瀏覽)，請將下載(或瀏覽)的超鏈接粘接至瀏覽器地（ 推薦MYIE或GREENBORWSER）址欄后按回車.若不出意外，此時應(yīng)能下載.
若下載中出現(xiàn)了問題，請參考:
http://blog.csdn.net/emilmatthew/archive/2006/04/08/655612.aspx
Trackback: http://tb.blog.csdn.net/TrackBack.aspx?PostId=672276
聯(lián)系客服
