一维均匀分布随机数序列的产生方法.doc

一维均匀分布随机数序列的产生方法

一维均匀分布随机数序列的产生方法
【摘要】 利用混沌的随机数产生算法和线性同余发生器以及MATLAB产生一维均匀分布随机数序列.经过检验,随机数列的统计性质有了很大提高，
【关键词】 混沌；线性同余发生器；MATLAB；随机数
1 引言
随机数在信息加密、数值运算及医学中基因序列分析等研究中有着广泛的应用。比如数值运算中,Monte Carlo方法占有重要的地位,随机数是该方法的基础.随机数的质量影响了信息的安全和计算结果的精度。特别是一些安全级别比较高的应用,对随机数提出了很高的要求。随机数可由硬件和软件两种方式产生。在计算机中广泛使用的是软件方式,通过计算机利用数学模拟随机过程产生随机数。此方法有着自身的不足,数据之间有着关联性,存在周期,并非真正的随机数,因此被成为伪随机数。
生成随机数的方法繁多，从产生机理来说，可分为数学方法和物理方法两种，其所产生的随机数分别被称之为伪随机数和真随机数，前者易被破解，后者取自物理世界的真实随机源，难以破解，但这并不代表基于真随机源产生的随机数质量就很高，要取决于产生算法如何利用这个真随机源，相反的，许多用数学方法产生的随机数质量比较好。因此，若能将数学方法和物理方法结合起来，则可能产生高质量的真随机数。常见的产生随机数的方法有【1】线性同余法（LCG,Linear Congruent Generators）、Tarsworthe位移计数器法、Fibonacci延迟产生器法等。为了克服以上方法的缺陷,人们还发展了许多新的方法。组合发生器就是著名的一种。它是将两个随机数发生器进行组合,以一种发生器产生一个随机数列,再用另一个随机数发生器对随机数列进行重修排列,得到一个更为独立,周期更长的随机数列。已有一些利用混沌序列转换伪随机数列的报道【2】,文献【3】虽然提出了一种由logistic映射构造具有均匀性数列的好方法,但数据之间的独立性较差。本研究中提出了一种新的方法,利用混沌算法【4】和线性同余发生器相组合得到随机数列,并就数据的均匀性和独立性进行了检验。
从实现方法来说，有以软件为主、以硬件为主以及软硬结合等方法【5】。
相比于伪随机数发生器的研究而言，真随机数发生器的研究还相当初步。设计一个真随机数发生器包括两步：首先是获取真随机源；然后是利用真随机源依照特定的数学方法获得真随机数。
2 理论基础
一维均匀分布随机数的产生
2.1 算法1
在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数.该随机数是平均在0~RAND_MAX之间平均分布的，RAND_MAX是一个常量，在VC6.0环境下是这样定义的：
#define RAND_MAX 0x7fff
它是一个short 型数据的最大值，如果要产生一个浮点型的随机数，可以将rand()/1000.0这样就得到一个0~32.767之间平均分布的随机浮点数.如果要使得范围大一点，那么可以通过产生几个随机数的线性组合来实现任意范围内的平均分布的随机数.例如要产生-1000~1000之间的精度为四位小数的平均分布的随机数可以这样来实现.先产生一个0到10000之间的随机整数.方法如下 ：
int a = rand()%10000;

然后保留四位小数产生0~1之间的随机小数：
double b = (double)a/10000.0;
然后通过线性组合就可以实现任意范围内的随机数的产生，要实现-1000~1000内的平均分布的随机数可以这样做：
double dValue = (rand()%10000)/10000.0*1000-(rand()%10000)/10000.0*1000;
则dValue就是所要的值.
但是，上面的式子化简后就变为：
double dValue = (rand()%10000)/10.0-(rand()%10000)/10.0;
这样一来，产生的随机数范围是正确的，但是精度不正确了，变成了只有一位正确的小数的随机数了，后面三位的小数都是零，显然不是我们要求的，什么原因呢，又怎么办呢.
先找原因，rand()产生的随机数分辨率为32767，两个就是65534，而经过求余后分辨度还要减小为10000，两个就是20000而要求的分辨率为1000*10000*2=20000000，显然远远不够.下面提供的方法可以实现正确的结果：
double a = (rand()%10000) * (rand()%1000)/10000.0;
double b = (rand()%10000) * (rand()%1000)/10000.0;
double dValue = a-b;
则dValue就是所要求的结果.在下面的函数中可以实现产生一个在一个区间之内的平均分布的随机数，精度是4位小数.
double AverageRandom(double min,double max)
{
int minInteger = (int)(min*10000);
int maxInteger = (int)(max*10000);
int randInteger = rand()*rand();
int diffInteger = maxInteger - minInteger;
int resultInteger = randInteger % diffInteger + minInteger;
return resultInteger/10000.0;
}
但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随机种子，如果不对随机种子进行初始化，那么计算机有一个缺省的随机种子，这样每次递推的结果就完全相同了，因此需要在每次程序运行时对随机种子进行初始化，在vc中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时间来作为随机种子，因为系统时间可以保证它的随机性.
调用方法是srand(GetTickCount()),但是又不能在每次调用rand()的时候都用srand(GetTickCount())来初始化，因为现在计算机运行时间比较快，当连续调用rand()时，系统的时间还没有更新，所以得到的随机种子在一段时间内是完全相同的，因此一般只在进行一次大批随机数产生之前进行一次随机种子的初始化.下面的代码产生了400个在-1~1之间的平均分布的随机数.
double dValue[400];
srand(GetTickCount());

for(int i= 0;i < 400; i++)
{
double dValue[i] = AverageRandom(-1,1);
}
用该方法产生的随机数运行结果如图1所示：

图1 400个-1~1之间平均分布的随机数
2.2 算法2：用同余法产生随机数
同余法简称为LCG(Linear Congruence Gener-ator)，它是Lehmer于1961年提出来的．同余法利用数论中的同余运算原理产生随机数．同余法是目前发展迅速且使用普遍的方法之一．
同余法(LCG)递推公式为
（n=1，2，…)， (1)
其中 ，a，c均为正整数．只需给定初值x.，就可以由式(1)得到整数序列{ }，对每一 ，作变换 = ／m，则{ }(n=1，2，…)就是[0，1)上的一个序列．如果{ }通过了统计检验，那么就可以将 作为[0，1)上的均匀分布随机数．
在式(1)中，若c=0，则称相应的算法为乘同余法，并称口为乘子；若c≠0，则称相应的算法为混合同余法．同余法也称为同余发生器，其中称为种子．
由式(1)可以看出，对于十进制数，当取模m=10 (k为正整数)时，求其同余式运算较简便．例如36=31236(mod10 )，只要对21236从右截取k=2位数，即得余数36．同理，对于二进制数，取模m=2 时，求其同余式运算更简便了．
电子计算机通常是以二进制形式表示数的．在整数尾部字长为L位的二进制计算机上，按式(1)求以m为模的同余式时，可以利用计算机具有的整数溢出功能．设L为计算机的整数尾部字长，取模m=2 ，若按式(1)求同余式时，显然有

这里[x]是取x的整数部分．在电子计算机上由 求 时，可利用整数溢出原理．不进行上面的除法运算．实际上，由于计算机的整数尾部字长为L，机器中可存放的最大整数为2 -1，而此时a +c≥m≥2 -1，因此a +c在机器里存放时占的位数多于L位，于是发生溢出，只能存放 的右后L位．这个数值恰是模m=2 的剩余，即xn,这就减少了除法运算，而实现了求同余式．经常取模m=2 (L为计算机尾部字长)，正是利用了溢出原理来减少除法运算.
由式(1)产生的 (n=1,2,……)，到一定长度后，会出现周而复始的周期现象，即{ }可以由其某一子列的重复出现而构成，这种重复出现的子列的最短长度称为序列 的周期．由式(1)不难看出，{ }中两个重复数之间的最短距离长度就是它的周期，用T代表周期．周期性表示一种规律性，它与随机性是矛盾的．
因此，通常只能取{ }的一个周期作为可用的随机序列．这样一来，为了产生足够多的随机数，就必须{ }的周期尽可能地大．由前所述，一般取m=2 ，这就是说模m已取到计算机能表示的数的最大数值，意即使产生的随机数列{ }的周期达到可能的最大数值，如适当地选取参数 ，a，c等，还可能使随机数列{ }达到满周期.
2.3 混沌算法
混沌是一种确定性系统中出现的类似随机的过程。它首先是由洛伦兹在流体热对流的简化模型计算中观察到的。混沌现象可以用它某些非线性函数的迭代（映射）来表示类似随机的输出。一维迭代Logistic方程是一个很简单的非线性抛物线函数,可以表示为：
xi+1=λxi(1-xi), (0≤xi≤1, 0≤λ≤4) (1)
xi 是i 次迭代的值。研究表明当λ>3.57 时,进入混沌状态。特征是表现为初值敏感性,既使初值相差非常小,经过几十次迭代，得到的值也会相差非常大，并且毫无关联。当λ=4 时,序列xi 的概率分布服从
P(xi)=1πxi(1-xi)(2)
所以xi 序列分布并不均匀，通过代换
ri=2arcsinxiπ(3)
就可以得到均匀的ri 数列。
2.4 线性同余发生器
线性同余发生器是现在使用最多的随机数发生器之一，该方法利用同余运算产生随机数。
Ii+1=(aIi+b)mod m
xi=Iim (4)
其中m 为模数, a为乘子, b为增量,并且a,b,m,I1 均为非负整数。如何选择他们就决定了产生器的质量。在计算中取a=16807,b=0,m=231-1,I1=12546863 。显然,由公式(4)得到的Ii 满足：0≤Ii≤ m,所以0≤xi≤1 。线性同余法的优点就是速度快,每次调用只需几个操作步骤即可,因此最为常用.缺点就是数据序列的有着关联性,表现为高维稀疏网格结构【6】。特别是当a 和m 的值较小时,这种现象更为明显。

2.5 组合发生器
20世纪70年代就有人提出用用组合发生器来产生高质量的伪随机数。有文献【7】对此进行了理论分析并指出：几个独立且近似的随机变量的线性组合也是一个近似均匀的随机变量,其分布比组成它的任何一个变量更接近U（0，1）。
该算法是分别用两个发生器各自产生一组随机数列,然后再把它们进行相互组合,以期得到一个统计结果更好的随机数列。具体步骤如下：
①N个不同的随机数发生器各自产生一组随机数列{ Xi（j)(i=1,2,…,L, j=1,2,…,N)。
②令c1,c2,…,cN 为N个非负整数。Ui=∑Nj=1c(j)Xi(j)ri=Ui mod 1 i=1,2,… (5)
按上述算法把混沌算法和线性同余产生的随机数组合在一起。在计算中,取c1=2,c2=5。每个发生器都有其各自的周期Period(Xi) ,如果它们互素,则组合后周期是各自周期的最小公倍数（LCM）。Period(Xi)=LCM(Period(Xi(1)）,Period(Xi(2),…)。因此组合发生器产生数据的周期会大大增加,实际应用中可视为无穷大。
3 随机数的统计检验
随机数的好坏是通过各种统计检验来确定的,这些检验包括均匀性检验、独立性检验、参数检验等。其中最基本的是均匀性检验和独立性检验。在随机数的检验中独立性检验是首要的,如果不存在独立性就无所谓随机。也希望得到的随机数在区间［0，1］上出现的概率是相同的,所以均匀性检验也起着决定作用。
3.1 均匀性检验
假设区间［0，1］上的随机数列{Xi}(i=1,2,…,n)。如果随机数分布均匀,则把区间分成k 个相等的子区间后,落在每个子区间的随机数个数ni 近似等于n / k 。统计量χ2 按定义为:
χ2=∑ki=1(ni-n / k)2n / k (6)
χ2 应服从χ2(k-1) 分布,在给定一个显著水平α 下,可查表得到临界值χα2(k-1)。如果χ2>χ2α ,则拒绝均匀假设。
3.2 独立性检验
对独立性检验中的相关检验基于相关系数r，把随机数列{ri}(i=1,2,…,2n) 分成两部分{Xi}(i=1,3,…,2n-1),{Yi}(i=2,4,…,2n} ,相关系数r 为:
r=1n∑ni=1(Xi-)(Yi-)1n∑ni=1(Xi-)2 1n∑ni=1(Yi-)2 (7)
当ρ=(1.96)2+n-2 |r|≥1.96 时, 认为线性回归效果显著,拒绝独立检验。
3.3 参数检验
随机数列{Xi}(i=1,2,…,n) ,其一阶矩、二阶矩和二阶中心矩分别为:
=1n∑ni=1Xi, 2=1n∑ni=1Xi2, s2=1n∑ni=1(Xi-12)2 (8)
对于均匀分布的数列,其值分别为12, 13,112如果与理论值没有显著差别，则通过参数检验。
检验中,3种发生器分别生成100000个数据,对其统计性质进行了分析。取检验的置信度α=0.05 ,在均匀检验中，检验区间k=500 ,χ2(500)>552.78 时拒绝均匀假设【8】；独立性检验中ρ=(1.96)2+n-2 |r|≥1.96,拒绝独立假设。表1 随机数列的统计结果 从表1的统计结果可以看出,组合后的发生器表现出了更为优秀的统计结果.在数据的独立性和均匀性上都大大增加。
混沌算法虽然可以通过均匀性检验,但独立性较差,没有通过独立性检验,不宜单独作为随机数发生器使用.如果把数列中相邻两个数据作为（x,y) ,分别作出各个数列的二维散布图。不难看出,图1混沌算法得到的随机数列,数据之间有着强烈的关联,数据之间并不独立。图2和图3数据分布则比较均匀。

4 结论
计算机模拟、信息加密或Monte Carlo方法成功的基础是实现真正的随机抽样,随机抽样的基础是随机数。从以上统计结果可以看出,两种发生器组合之后,得到的随机数列具有更好的统计结果,在保持均匀性的基础上,独立性有了较大的提高。并且组合后,随机数列的周期可视为无穷大。因此该方法具有良好的统计性质,符合随机数的要求。
【参考文献】
【1】杨自强,魏公毅. 产生伪随机数的若干新方法.数值计算机应用,
2001,3:210～216.
【2】 孙霞,吴自勤,黄畇,分形原理及其应用.中国科学技术大学出版社 2003,1～22.
【3】 盛利元,肖燕予,盛喆,将混沌序列变换成均匀伪随机序列的普适算法.
物理学报,2008,57：4407～4412.
【4】周燕,覃朝玲用混沌法产生随机数发生器的研究.西南师范大学学报,
2000,25：150～154.
【5】王相生, 甘骏人. 一种基于混沌的序列密码生成方法[J]. 计算机学报,
2002, 25(4): 351-356.
【6】罗平. 线性同余发生器的缺陷及其改进.计算机工程,1995,21：295～297.
【7】Deng L Y, George E O. Generation of Uniform Variate from Several Nearly
Uniformly Distributed Variables. Comm Statist Simu, 1990,19:145～154.
【8】王萍,许海洋.一种新的随机数组合发生器的研究 计算机技术与发展，
2006,16：79～81.