随机数产生原理及应用
EmilMatthew(EmilMatthew@126.com)
摘要:
本文简述了随机数的产生原理，并用C语言实现了迭代取中法，乘同余法等随机数产生方法，同时，还给出了在符合某种概率分布的随机变量的产生方法。
关键词: 伪随机数产生，概率分布
1前言:
在用计算机编制程序时，经常需要用到随机数，尤其在仿真等领域，更对随机数的产生提出了较高的要求，仅仅使用C语言类库中的随机函数已难以胜任相应的工作。本文简单的介绍随机数产生的原理及符合某种分布下的随机变量的产生，并用C语言加以了实现。当然，在这里用计算机基于数学原理生成的随机数都是伪随机数。
注:这里生成的随机数所处的分布为0-1区间上的均匀分布。我们需要的随机数序列应具有非退化性，周期长，相关系数小等优点。
2.1迭代取中法：
这里在迭代取中法中介绍平方取中法,其迭代式如下:
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)
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;
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;
}
前一百个测试生成的随机数序列:
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
各参数意义及各步的作用可参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()
{
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上的均匀分布,不过在某些数据段还有些起伏。
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为任意非负整数
相应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;
}
图3: 乘同余法生成的300随机数的产生序列图
图4: 乘同余法生成的300随机数的分布情况
由图4可以看出，该种随机数生成方法已相当接近0-1上的均匀分布。但在图3中可以看出它的一个致命的弱点，那就是随机数的生成在某一周期内成线性增长的趋势，显然，在大多数场合，这种极富“规律”型的随机数是不应当使用的。
下面的概率分布型随机变量的生成，均采用乘同余法或C函数库中的随机数来生成0-1区间上的随机数。
下面将C语言中的随机数生成序列图和Matlab中的随机数生成序列图列于下面，以作对比之用:
C语言生成的300个随机数的序列图
Matlab中rand函数生成的300个随机数的序列图
可以看出:乘同余法生成的随机数序列的随机性与上述两个标准库函数相接近。
3连续型随机变量的生成:
3.1反函数法
采用概率积分变换原理,对于随机变量X的分布函数F(X)可以求其反函数，得:
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
得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
试验结果:
图6:用反函数法生成的300随机数的指数分布情况
可以看出，生成的随机量较好的符合了指数分布特征。
3.2正态分布随机变量的生成:
正态分布在概率统计的理论及应用中占有重要地位，因此，能产生符合正态分布的随机变量就在模拟一类的工作中占有相当重要的地位。下面介绍两种方法。
3.2.1舍选法：