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

一维均匀分布随机数序列的产生方法【摘要】 利用混沌的随机数产生算法和线性同余发生器以及 MATLAB 产生一维 均匀分布随机数序列.经过检验,随机数列的统计性质有了很大提高，【关键词】混沌；线性同余发生器；MATLAB；随机数1 引言随机数在信息加密、数值运算及医学中基因序列分析等研究中有着广泛的应 用。比如数值运算中,Monte Carlo方法占有重要的地位，随机数是该方法的基础. 随机数的质量影响了信息的安全和计算结果的精度。特别是一些安全级别比较高 的应用,对随机数提出了很高的要求。随机数可由硬件和软件两种方式产生。在 计算机中广泛使用的是软件方式 ,通过计算机利用数学模拟随机过程产生随机 数。此方法有着自身的不足,数据之间有着关联性,存在周期,并非真正的随机数, 因此被成为伪随机数。生成随机数的方法繁多，从产生机理来说，可分为数学方法和物理方法两种， 其所产生的随机数分别被称之为伪随机数和真随机数，前者易被破解，后者取自 物理世界的真实随机源，难以破解，但这并不代表基于真随机源产生的随机数质 量就很高，要取决于产生算法如何利用这个真随机源，相反的，许多用数学方法 产生的随机数质量比较好。因此，若能将数学方法和物理方法结合起来，则可能 产生高质量的真随机数。常见的产生随机数的方法有【1】线性同余法(LCG,Linear Con grue nt Ge nerators)、Tarsworthe 位移计数器法、Fib on acci 延迟产生器法等。 为了克服以上方法的缺陷,人们还发展了许多新的方法。组合发生器就是著名的 一种。它是将两个随机数发生器进行组合,以一种发生器产生一个随机数列,再用 另一个随机数发生器对随机数列进行重修排列,得到一个更为独立,周期更长的随 机数列。已有一些利用混沌序列转换伪随机数列的报道【 2】 ,文献【3】虽然提 出了一种由 logistic 映射构造具有均匀性数列的好方法,但数据之间的独立性较 差。本研究中提出了一种新的方法,利用混沌算法【4】和线性同余发生器相组合 得到随机数列,并就数据的均匀性和独立性进行了检验。从实现方法来说，有以软件为主、以硬件为主以及软硬结合等方法【5】。相比于伪随机数发生器的研究而言，真随机数发生器的研究还相当初步。设 计一个真随机数发生器包括两步：首先是获取真随机源；然后是利用真随机源依 照特定的数学方法获得真随机数。2 理论基础 一维均匀分布随机数的产生2.1 算法 1在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数.该随机 数是平均在 0RAND_MAX 之间平均分布的， RAND_MAX 是一个常量，在 VC6.0 环境下是这样定义的：#define RAND_MAX 0x7fff它是一个 short 型数据的最大值，如果要产生一个浮点型的随机数，可以将 rand()/1000.0这样就得到一个032.767之间平均分布的随机浮点数.如果要使得 范围大一点，那么可以通过产生几个随机数的线性组合来实现任意范围内的平均 分布的随机数.例如要产生-10001000 之间的精度为四位小数的平均分布的随机 数可以这样来实现.先产生一个 0 到 10000 之间的随机整数.方法如下 ： int a = rand()%10000;然后保留四位小数产生 01 之间的随机小数：double b = (double)a/10000.0; 然后通过线性组合就可以实现任意范围内的随机数的产生，要实现 -10001000 内的平均分布的随机数可以这样做：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()的时候都用 sra nd(GetTickCou nt()来初始化，因为现在计算机运行时间比较快，当连续调用 ran d()时，系统的时间还没有更新，所以得到的随机种子在一段时间内是完全相 同的，因此一般只在进行一次大批随机数产生之前进行一次随机种子的初始化 . 下面的代码产生了400个在-11 之间的平均分布的随机数.double dValue400;srand(GetTickCount();for(int i= 0;i 400; i+)double dValuei = AverageRandom(-1,1); 用该方法产生的随机数运行结果如图 1 所示：图 1 400 个 -11 之间平均分布的随机数2.2 算法 2：用同余法产生随机数同余法简称为 LCG(L in ear Con grue nee Gen er-ator),它是 Lehmer 于 1961 年提 出来的同余法利用数论中的同余运算原理产生随机数同余法是目前发展迅速 且使用普遍的方法之一同余法 (LCG) 递推公式为x 二(ax + c ) ( m md(n=1，2,)，(1)nn-1其中，a, c均为正整数.只需给定初值x.,就可以由式得到整数序列 , 对每一，作变换=/m,贝肝(n=1, 2,)就是0, 1)上的一个序列.如果 通过了统计检验,那么就可以将 作为0, 1)上的均匀分布随机数.在式中，若c=0,则称相应的算法为乘同余法，并称口为乘子；若cH0, 贝称相应的算法为混合同余法.同余法也称为同余发生器,其中称为种子.由式(1)可以看出,对于十进制数,当取模 m=10 (k 为正整数)时,求其同余 式运算较简便.例如36=31236(mod10 ),只要对21236从右截取k=2位数，即得 余数 36.同理,对于二进制数,取模 m=2 时,求其同余式运算更简便了.电子计算机通常是以二进制形式表示数的.在整数尾部字长为L位的二进制 计算机上,按式(1)求以 m 为模的同余式时,可以利用计算机具有的整数溢出功 能.设L为计算机的整数尾部字长，取模m=2，若按式(1)求同余式时，显然有当ax + c m时，贝Ux 二 ax + c-m(ax + c)/m.n -1nn -1n-1这里X是取x的整数部分.在电子计算机上由求 时，可利用整数溢出 原理.不进行上面的除法运算.实际上，由于计算机的整数尾部字长为L,机 器中可存放的最大整数为2 -1,而此时a +c$m$2 -1,因此a +c在机器里存放 时占的位数多于L位，于是发生溢出，只能存放的右后L位.这个数值恰是模 m=2的剩余，即xn,这就减少了除法运算，而实现了求同余式.经常取模m=2 (L 为计算机尾部字长)，正是利用了溢出原理来减少除法运算.由式产生的(n=1,2,)，到一定长度后，会出现周而复始的周期现象, 即 可以由其某一子列的重复出现而构成，这种重复出现的子列的最短长度 称为序列 的周期.由式(1)不难看出， 中两个重复数之间的最短距离长度就是 它的周期，用 T 代表周期.周期性表示一种规律性，它与随机性是矛盾的.因此，通常只能取 的一个周期作为可用的随机序列.这样一来，为了产 生足够多的随机数，就必须的周期尽可能地大.由前所述，一般取m=2，这 就是说模m已取到计算机能表示的数的最大数值，意即使产生的随机数列的 周期达到可能的最大数值，如适当地选取参数，a，c等，还可能使随机数列 达到满周期.2.3 混沌算法 混沌是一种确定性系统中出现的类似随机的过程。它首先是由洛伦兹在流体热对流的简化模型计算中观察到的。混沌现象可以用它某些非线性函数的迭代(映射)来表示类似随机的输出。一维迭代Logistic方程是一个很简单的非线性 抛物线函数,可以表示为：xi+1=Axi(1-xi), (0xi1, 0X3.57时，进入混沌状态。特征是表现为初值敏 感性,既使初值相差非常小,经过几十次迭代，得到的值也会相差非常大，并且毫 无关联。当入=4时，序列xi的概率分布服从P(xi)=1nxi(1-xi)(2)所以 xi 序列分布并不均匀，通过代换ri=2arcs in xin(3)就可以得到均匀的 ri 数列。2.4 线性同余发生器 线性同余发生器是现在使用最多的随机数发生器之一，该方法利用同余运算产生随机数。Ii+1=(aIi+b)mod mxi=Iim (4)其中 m 为模数, a 为乘子, b 为增量,并且 a,b,m,I1 均为非负整数。如何选择他们 就决定了产生器的质量。在计算中取a=16807,b=0,m=231-1,l1=12546863。显然, 由公式得到的Ii满足：0ii m,所以0xix2a，则拒绝均匀假设。3.2 独立性检验对独立性检验中的相关检验基于相关系数r,把随机数列ri(i=1,2,.,2n)分成两部分Xi(i=1,3,.,2n-1),Yi(i=2,4,.,2n，相关系数 r 为：r=1 nn i=1(Xi-)(Yi-)1 nn i=1(Xi-)2 1nn i=1(Yi-)2 (7)当p=(1.96)2+n-2 |r|1.96时，认为线性回归效果显著，拒绝独立检验。3.3 参数检验随机数列Xi(i=1,2,.,n)，其一阶矩、二阶矩和二阶中心矩分别为：=1nni=1Xi, 2=1 nni=1Xi2, s2=1 nni=1(Xi-12)2 (8)对于均匀分布的数列,其值分别为12, 13,112如果与理论值没有显著差别,则通过 参数检验。检验中,3种发生器分别生成100000个数据,对其统计性质进行了分析。取检验的置信度a=0.05，在均匀检验中，检验区间k=500 ,x2(500)552.78时拒绝均匀假设 【8】；独立性检验中p=(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,122. 【 3】 盛利元,肖燕予,盛喆,将混沌序列变换成均匀伪随机序列的普适算法. 物理学报,2008,57：44074412.【4】周燕,覃朝玲用混沌法产生随机数发生器的研究.西南师范大学学报, 2000,25：150154.【5】王相生，甘骏人.一种基于混沌的序列密码生成方法J.计算机学报，2002, 25(4): 351-356.【6】罗平. 线性同余发生器的缺陷及其改进.计算机工程,1995,21：295297.【 7 】 Deng L Y, George E O. Generation of Uniform Variate from Several Nearly Uniformly Distributed Variables. Comm Statist Simu, 1990,19:145154. 【8】王萍,许海洋.一种新的随机数组合发生器的研究 计算机技术与发展， 2006,16：7981.