随机数的算法分析 
摘 要:
我们知道,不管是在计算机编程,还是网站设计、分析实际问题,随机数都有广泛的应用。本文首先讨论了生成0-1之间均匀分布随机数的一些算法,进而给出了由0-1之间均匀分布生成指数分布、正态分布、χ2分布、二项分布、泊松分布的一般算法,并通过方差分析、均值检验、χ2检验对所得数据进行分析,最后得出满足一般要求的一系列随机数。
 
关键词:
随机数,均匀分布,线性同余法,人字映射,平方取中法,正态分布,指数分布,χ2分布,二项分布,泊松分布,皮尔逊χ2检验法,C++
 
1-0:Microsoft VC++产生随机数的原理:
Srand ( )和Rand( )函数。它本质上是利用线性同余法,y=ax+b(mod m)。其中a,b,m都是常数。因此rand的产生决定于x,x被称为Seed。Seed需要程序中设定,一般情况下取系统时间作为种子。它产生的随机数之间的相关性很小,取值范围是0—32767(int),即双字节(16位数),若用unsigned int 双字节是65535,四字节是4294967295,一般可以满足要求。
1-1: 线性同余法:
?/P>
其中M是模数,A是乘数,C是增量,为初始值,当C=0时,称此算法为乘同余法;若C≠0,则称算法为混合同余法,当C取不为零的适当数值时,有一些优点,但优点并不突出,故常取C=0。模M大小是发生器周期长短的主要标志,常见有M为素数,取A为M的原根,则周期T=M-1。例如:

a=1220703125        
a=32719            (程序中用此组数)   
a=16807          
代码:
void main( )
{
const int n=100;
double a=32719,m=1,f[n+1],g[n],seed;
m=pow(2,31);
cout<<"设置m值为  "<<m-1<<endl;
cout<<"输入种子"<<endl;  //输入种子
cin>>seed;
f[0]=seed;    
    for(int i=1;i<=n;i++)    //线性同余法生成随机数
      {
         f[i]=fmod((a*f[i-1]),(m-1));
         g[i-1]=f[i]/(m-1);
         cout.setf(ios::fixed);cout.precision(6); //设置输出精度
         cout<<i<<"   "<<'\t'<<g[i-1]<<endl;
      }
}


结果分析:统计数据的平均值为:0.485653
统计数据的方差为:0.320576
 
1-2:人字映射
递推公式
?/P>
就是有名的混沌映射中的“人字映射”或称“帐篷映射”,它的非周期轨道点的分布密度函数:人字映射与线性同余法结合,可产生统计性质优良的均匀随机数。

for(int i=1;i<=n;i++)    //线性同余法生成随机数
      {
         f[i]=fmod((a*f[i-1]),m);
             if(f[i]<=m/2)     //与人字映射结合生成随机数
             {
                    f[i]=2*f[i];
             }
             else
             {
                    f[i]=2*(m-f[i])+1;
             }


1-3:平方取中法——冯·诺伊曼
1946年前后,由冯·诺伊曼提出,他的办法是去前面的随机数的平方,并抽取中部的数字。例如要生成10位数字,而且先前的值是5772156649,平方后得到33317792380594909201,所以下一个数是7923805949。

for(j=1;j<=n;j++)
      {
             i[j]=i[j-1]*i[j-1];  
        i[j]=i[j]/pow(10,5); 
        i[j]=fmod(i[j],pow(10,10));
        g[j]=i[j]/pow(10,10);
        cout.setf(ios::fixed);cout.precision(6); //设置输出精度
        cout<<j<<'\t'<<g[j]<<endl;
      }


二:任意分布随机数的生成
     利用(0,1)均匀分布的随机数可以产生任意分布的随机数。主要的方法有反函数法,舍选法,离散逼近法,极限近似法和随机变量函数法等。这里主要讨论了反函数法,当然对于具体分布函数可以采用不同的方法。
设随机变量X具有分布函数F(X),则对一个给定的分布函数值,X的值为
                                             
其中inv表示反函数。现假设r是(0,1)均匀分布的随机变量R的一个值,已知R的分布函数为
                             
因此,如果r是R的一个值,则X具有概率
 
也就是说如果 (r1,r2,...,rn)是R的一组值,则相应可得到的一组值
                   
具有分布。从而,如果我们已知分布函数的反函数,我们就可以从(0,1)分布的均匀分布随机数得到所需分布的随机数了。
1-4:指数分布:
指数分布的分布函数为:
x<0时,F(x)=0    ; ,F(x)=1-exp
利用上面所述反函数法,可以求得:  x= ln(1-y),这里不妨取常数 为1.

for(int j=0;j<n;j++)
 { 
              i=rand()%100;//产生从0-32767的任意一个值
             a[j]=double(i)/double(100); 
             a[j]=-log(a[j]);//  常数大于0,这里取1
          、、、、、、、


1-5:正态分布:
正态分布的概率密度是:
正态分布的分布函数是:
对于正态分布,利用反函数的方法来获取正态分布序列显然是很麻烦的,牵涉到很复杂的积分微分运算,同时为了方便,我们取,即标准正态分布。因此这里介绍了两种算法:
第一种:
Box和Muller在1958年给出了由均匀分布的随机变量生成正态分布的随机变量的算法。设U1, U2是区间 (0, 1)上均匀分布的随机变量,且相互独立。令 

X1=sqrt(-2*log(U1)) * cos(2*PI*U2); 
  X2=sqrt(-2*log(U1)) * sin(2*PI*U2);  
那么X1, X2服从N(0,1)分布,且相互独立。
             p=rand()%100;//产生从0-32767的任意一个值
             b[j]=double(p)/double(100);
             a[j]=sqrt(-2*log(a[j]))*cos(2*3.1415926*b[j]);


第二种:
近似生成标准正态分布,独立同分布的多个随机变量和的分布趋近于正态分布,取k个均匀分布的(0,1)随机变量,,…… ,则它们的和近似服从正态分布。
  实践中,取k=12,(因为D( )=1/12),则新的随机变量y=x1+x2+...+x12-6,可以求出数学期望E(y)=0,方差D(y)=12*1/12=1,因此可以近似描述标准正态分布。