随机数产生原理及应用

mPoint.x=miu1+sigma1*x; mPoint.y=miu2+sigma2*x;
Xn=fmod(Lamda*Xn,M);//here can‘s use % Rn=Xn/M; return Rn; } 初始化段,应有: M=pow(2,35)-31;
图 1：乘同余法生成的 300 随机数的产生序列图 图 2: 乘同余法生成的 300 随机数的分布情况
可以看到，该随机数生成方法所生成的随机序列比较符合 0-1 上的均匀分布, 不过在某些数据段还有些起伏。
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));
图 6:用反函数法生成的 300 随机数的指数分布情况 可以看出，生成的随机量较好的符合了指数分布特征。 3.2 正态分布随机变量的生成: 正态分布在概率统计的理论及应用中占有重要地位，因此，能产生符合正态 分布的随机变量就在模拟一类的工作中占有相当重要的地位。下面介绍两种方 法。 3.2.1 舍选法： 这种方法便捷而有效，且具有一定的代表性，其基本思路是: 在概率密度的函数图像的外围画一个大框，然后在这个框内部产生随机点 （rx,ry），根据是否落在概率密度函数的下方，来决定是否要留下这个点。
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();
经过一定的计算变行，符合二维的正态分布的随机变量的生成可按下面的方 法进行：
1)产生位于 0-1 区间上的两个随机数 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)正态分布,则按上述方法产生 x 后，x’=miu+sigma*x 由于采用基于乘同余法生成的 0-1 上的随机数的正态分布随机数始终无法 能过正态分布总体均值的假设检验。而采用 C 语言的库函数中的随机数生成函 数 rand()来产生 0-1 上的随机数，效果较为理想。 关键程序段(funNorm 返回一维的正态分布，而 funNorm2 则生成二维的随机 分布): float funNorm(float miu,float sigma) {
得 Fi=G(R)=2*PI*R ,其中，R 为 0-1 区间上的均匀分布的随机数. 程序略
试验结果:
图 5:用反函数法生成的 300 随机数的平均分布情况 由于这里相当对 0-1 上的分布进行线性变换，所以变换后仍呈均匀分布是显 然的。 3.1.2 指数分布: 指数分布的分布函数为: x<0 时,F(x)=0 ; x>=0,F(x)=1-exp(-lamda*x) 利用反函数法，可以求得: x=-lnR/lamda 试验结果:
Xn+1=(Xn^2/10^s)(mod 10^2s) Rn+1=Xn+1/10^2s 其中，Xn+1 是迭代算子，而 Rn+1 则是每次需要产生的随机数 。 第一个式子表示的是将 Xn 平方后右移 s 位，并截右端的 2s 位。 而第二个式子则是将截尾后的数字再压缩 2s 倍，显然:0=<Rn+1<=1. 这样的式子的构造需要深厚的数学(代数，统计学，信息学)功底，这里只是拿来 用一下而已，就让我们站在大师的肩膀上前行吧。 迭代取中法有一个不良的性就是它比较容易退化成 0. 平方取中法的实现: #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)
随机数产生原理及应用
摘要:
EmilMatthewΒιβλιοθήκη BaiduEmilMatthew@126.com)
本文简述了随机数的产生原理，并用 C 语言实现了迭代取中法，乘同余法等
随机数产生方法，同时，还给出了在符合某种概率分布的随机变量的产生方法。
关键词: 伪随机数产生，概率分布
1 前言: 在用计算机编制程序时，经常需要用到随机数，尤其在仿真等领域，更对随 机数的产生提出了较高的要求，仅仅使用 C 语言类库中的随机函数已难以胜任 相应的工作。本文简单的介绍随机数产生的原理及符合某种分布下的随机变量的 产生，并用 C 语言加以了实现。当然，在这里用计算机基于数学原理生成的随 机数都是伪随机数。 注:这里生成的随机数所处的分布为 0-1 区间上的均匀分布。我们需要的随机 数序列应具有非退化性，周期长，相关系数小等优点。 2.1 迭代取中法： 这里在迭代取中法中介绍平方取中法,其迭代式如下:
return 0; } 前一百个测试生成的随机数序列:
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
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));
Xi=G(Ri) 其中,Ri 为一个 0-1 区间内的均匀分布的随机变量. F(X)较简单时，求解较易，当 F(X)较复杂时，需要用到较为复杂的变换技巧。 3.1.1 平均分布: 例:已知炮弹对目标的方位角 Fi 在 0-2*P 内均匀分布，试用（0，1）均匀随 机数变换，模拟弹着点方位角的抽样值 Fi. 解: R=F(Fi)=Fi/2*PI
下面的概率分布型随机变量的生成，均采用乘同余法或 C 函数库中的随机数 来生成 0-1 区间上的随机数。
下面将 C 语言中的随机数生成序列图和 Matlab 中的随机数生成序列图列于 下面，以作对比之用:
C 语言生成的 300 个随机数的序列图
Matlab 中 rand 函数生成的 300 个随机数的序列图 可以看出:乘同余法生成的随机数序列的随机性与上述两个标准库函数相接 近。 3 连续型随机变量的生成: 3.1 反函数法 采用概率积分变换原理,对于随机变量 X 的分布函数 F(X)可以求其反函数， 得:
2.3 混合同余法： 混合同余法是加同余法和乘同余法的混合形式,其迭代式如下: Xn+1=(Lamda*Xn+Miu)%M Rn+1=Xn/M 经前人研究表明，在 M=2^q 的条件下，参数 lamda,miu,X0 按如下选取，周 期较大，概率统计特性好: Lamda=2^c+1,c 取 q/2 附近的数 Miu=(1/2+sqrt(3))/M X0 为任意非负整数
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; } /*测试主程序，注意，这里只列举一次测试主程序，以下不再重复*/ int main() { int i; FILE * debugFile;
相应 C 程序关键代码段: //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; }
容易看出其易退化成 0 的缺点.
2.2 乘同余法： 乘同余法的迭代式如下: Xn+1=Lamda*Xn(mod M) Rn+1=Xn/M 各参数意义及各步的作用可参 2.1 当然，这里的参数的选取是有一定的理论基础的，否则所产生的随机数的周 期将较小，相关性会较大。 经过前人检验的两组性能较好的素数取模乘同余法迭代式的系数为: 1） lamda=5^5,M=2^35-31 2） lamda=7^5,M=2^31-1 相应 C 程序关键代码段: double long M;//请注意，这里一定要用到 double long,否则计算 2^32 会溢 出 float MyRnd() {
图 3: 乘同余法生成的 300 随机数的产生序列图 图 4: 乘同余法生成的 300 随机数的分布情况
由图 4 可以看出，该种随机数生成方法已相当接近 0-1 上的均匀分布。但在 图 3 中可以看出它的一个致命的弱点，那就是随机数的生成在某一周期内成线性 增长的趋势，显然，在大多数场合，这种极富“规律”型的随机数是不应当使用 的。
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) {