随机数的产生

随机数的产生摘要本文研究了连续型随机数列的产生，先给出了均匀分布的随机数的产生算法，在通过均匀分布的随机数变换得到其他连续型随机数的产生算法.在v c 环境下，我们给出了产生均匀分布随机数的算法，然后探讨了同余法的理论原理.通过均匀随机数产生其他分布的随机数，我们列举了几种通用算法，并讨论各个算法的优缺点，最后以正态分布为例验证高效舍选法的优势. 正文一、 随机数与伪随机数随机变量η的抽样序列12,,n ηηη ，…称为随机数列．如果随机变量η是均匀分布的，则η的抽样序列12,,n ηηη ，…称为均匀随机数列；如果随机变量η是正态分布的随机变量则称其抽样序列为正态随机数列．比如在掷一枚骰子的随机试验中出现的点数x 是一个随机变量，该随机变量就服从离散型均匀分布，x 取值为1，2，3，4，5，6，取每个数的概率相等均为1／6．如何得到x 的随机数?通过重复进行掷骰子的试验得到的一组观测结果12,,,n x x x 就是x 的随机数．要产生取值为0，1，2，…，9的离散型均匀分布的随机数，通常的操作方法是把10个完全相同的乒乓球分别标上0，1，2，…，9，然后放在一个不透明的袋中，搅拦均匀后从中摸出一球记号码1x 后放回袋中，接着仍将袋中的球搅拌均匀后从袋中再摸出一球记下号码2x 后再放回袋中，依次下去，就得到随机序列12,,,n x x x ．通常称类似这种摸球的方法产生的随机数为真正的随机数．但是，当我们需要大量的随机数时，这种实际操作方法需要花费大量的时间，通常不能满足模拟试验的需要，比如教师不可能在课堂上做10000次掷硬币的试验，来观察出现正面的频率．计算机可以帮助人们在很短时间产生大量的随机数以满足模拟的需要，那么计算机产生的随机数是用类似摸球方法产生的吗?不是．计算机是用某种数学方法产生的随机数，实际上是按照一定的计算方法得到的一串数，它们具有类似随机数的性质，但是它们是依照确定算法产生的，便不可能是真正的随机数，所以称计算机产生的随机数为伪随机数．在模拟计算中通常使用伪随机数．对这些伪随机数，只要通过统计检验符合一些统计要求，如均匀性、随机性等，就可以作为真正的随机数来使用，我们将称这样产生的伪随机数为随机数．在计算机上用数学方法产生随机数的一般要求如下：1)产生的随机数列要有均匀性、抽样的随机性、试验的独立性和前后的一致性．2)产生的随机数列要有足够长的周期，以满足模拟实际问题的要求．3)产生随机数的速度要快，占用的内存少．计算机产生随机数的方法内容是丰富的，在这里我们介绍几种方法，计算机通常是先产生[0，1]区间上均匀分布的随机数，然后再产生其他分布的随机数．二、均匀分布随机数的产生2.1 算法1在v c的环境下，为我们提供了库函数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;}但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随机种子，如果不对随机种子进行初始化，那么计算机有一个缺省的随机种子，这样每次递推的结果就完全相同了，因此需要在每次程序运行时对随机种子进行初始化，在v c中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时间来作为随机种子，因为系统时间可以保证它的随机性.调用方法是srand(GetT ickCount()),但是又不能在每次调用rand()的时候都用srand(GetT ickCount())来初始化，因为现在计算机运行时间比较快，当连续调用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 于1951年提出来的．同余法利用数论中的同余运算原理产生随机数．同余法是目前发展迅速且使用普遍的方法之一．同余法(LCG)递推公式为1()(mod )n n x ax c m -=+ （n=1，2，…)， (1)其中n x ，a ，c 均为正整数．只需给定初值x.，就可以由式(1)得到整数序列{n x }，对每一n x ，作变换n u =n x ／m ，则{n u }(n=1，2，…)就是[0，1)上的一个序列．如果{n u }通过了统计检验，那么就可以将n u 作为[0，1)上的均匀分布随机数．在式(1)中，若c=0，则称相应的算法为乘同余法，并称口为乘子；若c ≠0，则称相应的算法为混合同余法．同余法也称为同余发生器，其中0x 称为种子．由式(1)可以看出，对于十进制数，当取模m=10k(k 为正整数)时，求其同余式运算较简便．例如36=31236(mod102)，只要对21236从右截取k=2位数，即得余数36．同理，对于二进制数，取模m=2k 时，求其同余式运算更简便了．电子计算机通常是以二进制形式表示数的．在整数尾部字长为L 位的二进制计算机上，按式(1)求以m 为模的同余式时，可以利用计算机具有的整数溢出功能．设L 为计算机的整数尾部字长，取模m=2L ，若按式(1)求同余式时，显然有11111;[()/].n n n n n n n ax c m x ax c ax c m x ax c m ax c m -----+<=++≥=+-+当时，则当时，则这里[x]是取x 的整数部分．在电子计算机上由1n x -求n x 时，可利用整数溢出原理．不进行上面的除法运算．实际上，由于计算机的整数尾部字长为L ，机器中可存放的最大整数为2L -1，而此时a 1n x -+c ≥m ≥2L -1，因此a 1n x -+c 在机器里存放时占的位数多于L 位，于是发生溢出，只能存放n x 的右后L 位．这个数值恰是模m=2L 的剩余，即n x ．这就减少了除法运算，而实现了求同余式．经常取模m=2L (L 为计算机尾部字长)，正是利用了溢出原理来减少除法运算.由式(1)产生的n x (n=1,2,……)，到一定长度后，会出现周而复始的周期现象，即{n x }可以由其某一子列的重复出现而构成，这种重复出现的子列的最短长度称为序列n x 的周期．由式(1)不难看出，{n x }中两个重复数之间的最短距离长度就是它的周期，用T 代表周期．周期性表示一种规律性，它与随机性是矛盾的．因此，通常只能取{n x }的一个周期作为可用的随机序列．这样一来，为了产生足够多的随机数，就必须{n x }的周期尽可能地大．由前所述，一般取m=2L ，这就是说模m 已取到计算机能表示的数的最大数值，意即使产生的随机数列{n x }的周期达到可能的最大数值，如适当地选取参数0x ，a ，c 等，还可能使随机数列{n x }达到满周期.三、非均匀分布随机数的产生3.1 一般通用方法 3.1.1组合法组合法的基本思想是把预定概率密度函数f ( x ) 表为其它一些概率密度的线性组合.而这些概率密度的随机抽样容易产生.通过这种避难就易的手段我们也许可以达到较高的输出速度和较好的性能.若分布密度函数f ( x ) 能表为如下式(2)所示的函数项级数的和，1()()i i i f x p f x ∞==∑(2)其 中1i i p ∞=∑，诸f( x )皆为概率密度函数.则依如下步骤可产生分布为f ( x )一次抽样.( 1 ) 产生一个随机自然数I ， 使I 服从如下分布律:P ( I = i ) = p i i = 1 , 2 , 3……( 2 ) 产生服从f I ( x )的随机数0X 证明利用全概率公式，有:11()()()()()i i i i P x X x dx P Ii P x X x dx I i p f x dxf x dx∞=∞=<≤+==<≤+|===∑∑故X 服从f ( x ) 分布.我们以产生双指数(或拉普拉斯)分布的随机数为例来简单说明这种方法.双指数分布具有 概率密度函数f ( x ) = 0 . 5xe-， 如图2 :图2 双指数密度函数 f ( x ) 可表为:()0.5()0.5()l r f x f x f x =+ (3)其中()r f x 是指数分布，()l f x 是指数分布的对称分布.故产生双指数分布的抽样可按如下方法: 产生U 1 , U 2~U ( 0 , 1 ) ;若U 1 > 0 . 5 ， 则令X = I n U 2，否则X = - I n U 2. 在式(2) 中， 若i →∞， 有p i → 0 ，则可用函数列{()}i i p f x 的前有限项和逼近f ( x ).这是一种近似的方法，与通常的函数逼近原理相同.只要近似的精度 ( 在某种“精度”的意义之下) 达到要求，我们就可以采用近似的方法 .使用组合法时，各f i ( x ) 的抽样应该容易产生，故选用合适的概率密度函数族{ f i ( x )}把任意连续分布表为式(2) ，乃是使用组合法的关键.3.3.2 概率密度变换法这是一种比较新的通用随机数产生方法.其主要的目的是对一般的f(x)找出较好的覆盖函数以达到较高的效率.我们知道，对某一特定的概率密度f(x)，我们可以使用最优化技术找到好的覆盖函数.但对于一般情况，我们只能期望产生效率尚可的覆盖函数. H O R M A N N 用概率密度变换的方法生成一曲边梯形作为覆盖函数.其原理如下:使用一个变换函数T (x)把预定密度函数f ( x ) 变换为h ( x ) = T ( f ( x ) ) ，用一个分段线性函数l ( x )覆盖h ( x )，如图2 - 4 左图; h ( x ) 若是上凸的，则T 1-( l ( x ) )将是f ( x ) 的一个较好的覆盖函数，图3 概率密度变换法原理图这个方法在选择合适的T ( x ) ( l o g ( x ) 或1 / x a等) 后，能产生随机数包括了较多的分布类型.这个方法有较短的预处理时间，但需要较多的函数计算，不太适合硬件实现.此外，A h r e n s l用每段为常数的分段函数作为覆盖函数.L e y d o l d基于r a t i o - o f - u n i f o r m s 的方法也是一个通用算法.还有一种近似的方法，其产生的随机数与指定分布的随机数具有相同的前四阶矩，但概率分布不一定相同.这里就不详细介绍了.3.2 我们的方法当前的通用算法的问题是效率不能任意提高，不够灵活. 通常产生每个所需随机数X需以较大的概率计算f ( x )等函数.我们认为在速度要求非常高的场合，计算f ( x )是不利的，尤其以硬件进行函数计算是十分不利的.针对己有通用算法的不足，我们提出了基于组合法的通用算法.主要目的是尽可能地减少三角、指数、对数等超越函数的计算，以便硬件实现.产生任意连续分布随机数的高效舍选算法本文提出一种通用算法，可视需要使效率接近1 , 而且f ( x ) 的计算概率可任意小. 这些优点的取得是以长的预处理时间为代价的.在需要产生大量随机样本的场合( 例如通信系统的误码率测试，可能需要数小时乃至数天的仿真时间) ，本算法将有很大的优势，尽管有看法认为只有能用简单代码实现的算法才会被经常使用.3.2.1 算法原理假定预定的连续概率密度函数f ( x ) 为单峰的( 这是实际的大多数情况) ，已知其峰值点为m .一般f ( x ) 不关于x =m 对称，如图2 -5 .我们假定f ( x ) 定义在有界的区间[ a , b ] 上( 上文说过，对正态分布这类定义区间无限的情况，我们把这个区间取得足够大就可以了) . 直线X=m把f ( x ) 曲线与X轴所围面积分为左右两部分，我们把左右两部分各等分为K份，一共得到2 K个曲边梯形.并用2 K个矩形各自覆盖相应的曲边梯形.如图4所示( 图中各均分为四份) .R i , L i ( i = 1 , 2 , . . . , K -1 )是分点.并令R0=L=m,Lk = a，Rk= b .图4 均分f ( x ) 曲 线与X 轴所围面积我们的想法是利用舍选法的几何意义，分别在上述 2 K 个曲边梯形内均匀投点，从而使随机点在f ( x ) 曲线与x 轴所围的整个区域中均匀分布，这样即可产生f ( x ) 的抽样X . 而在曲边梯形内均匀投点可使用简单舍选法:先在各个矩形内均匀投点，再选出落于相应曲 边梯形内的点. 这种投点法浪费的点只位于各个矩形的一角， 显然效率大大高于简单舍选法.最为重要的是:随着K 的增大，效率会不断提高.另外，只有当投点位于曲边梯形的曲边之下时， 才需计算f ( x ) ，而且计算f ( x ) 概率是随着K 的增加而减小的. 我们每次“ 按概率”随机选中一个曲 边梯形进行投点. 这需要两步完成:先选择左边还是右边，再于此边的K 个曲边梯形中选择一个.这里的概率显然就是面积，这可以从以下的推导中看出来.为清晰起见，我们先阐述随机数的产生法，而把面积的均分这个预处理过程置于随后. 3.2.2 算法推导令()mP f x dx -∞=⎰为左边面积.则左边各曲边梯形面积皆为 P / K ，右边各曲边梯形面积皆为( ( I -P ) / K . f ( x ) 可表为:12111()()()KKi i i i P P f x f x f x KK==-=+∑∑(4)诸ji f ( x ) ( j = 1 , 2 ; i = 1 , 2 . . . k ) 皆为一腰为曲边的梯形形状的概率密度函数.依如下步骤可产生分布为f ( x ) 的一次抽样:S t e p l :产生一个随机自然数J ，使J 服从如下两点分布: P ( J = 1 ) = P , P ( J = 2 ) = 1 - P : S t e p 2 :产生一个随机自然数I ， 使I 服从如下均匀分布律:P ( I = i ) = 1 / K , i = 1 , 2 . . . . K ; S t e p 3 : 用基本舍选法产生概率密度为f ( x ) 的随机数X . 证明利用全概率公式，有:2111211()()()(,)1(()())()Kj i KKi i i i P x X x dx P Jj P I i P x X x dx J j I i P P f x f x dxKKf x dx====<≤+===<≤+∣==-=+=∑∑∑∑故x 服从 f ( x ) 分布.下面完整地描述这个方法: S t e p l( 产生J ) :S t e p l . l 产生[ 0 , 1 ] 上的均匀随机数U 1 ;S t e p 1 . 2若U 1 < P ，则返回J = 1 ， 否则返回J = 2 ; S t e p 2( 产生I ) :S t e p 2 . l 产生 [ 0 , I ] 上的均匀随机数U 2 ;S t e p 2 . 2 21;I kU x =+⎢⎥⎢⎥⎣⎦⎣⎦表示不大于x 的最大整数.产生 ji f ( x ) 的样本需区别j = 1 与j = 2 两种情况. 图2 - 6 示出j = 2 时一 典型的ji f ( x ) ， 用简单舍选法产生其抽样，覆盖函数为矩形. 首先产生一个[ 0 , R i ] 的均匀数， 如它属于[ 0 , R 1i -] 小无需再产生y 轴方向的均匀随机数，接受此均匀数即可;否则还需产生一个Y 轴方向的均匀随机数进行投点，那些落在曲边下方的点被接受，投在矩形右上角的点被舍弃.同理易得j = 1 时的产生法.图5 典型的ji f ( x ) j=2整个S t e p 3 如下: S t e p 3( 产生X ) : i f J = =1{ l o o p :产生[ 0 , 1 ] 上的均匀随机数U 3 , W = ( L 0 - L 1 ) U 3 + L 1 : i f W> L1i -，返回 X = W;e l s e { 产生[ O ， l ] 上的均匀随机数V ; if f ( W) - f ( L 1) < ( f ( L1j - ) - f ( L 1 ) ) V 返回X = W;e l s e 舍弃W ，重复l o o p ;} } e l s e{ l o o p : 产生[ 0 , 1 ] 上的 均匀随机数U 3 , W = ( R 1 - R o ) U 3 + R o ;i f W< R 1i -，返回 X = W; e l s e {产生[ 0 , 1 ] 上的均匀随机数V ;i f f ( W) - f ( R 1) C ( f ( R 1I - ) - f ( R 1) ) V , 返回X = W; e l s e 舍弃W ，重复l o o p ;} } 均匀随机数U 2 实际上可由U 1 变换得到， U 3 可由均匀数U2变换得到. 例如从U1 产生U 2 的方法是:当J = l 时， U 1 在[ 0 , P ] 上均匀分布， 故可令U 2 = U l / P ;当J = 2 时， U 1在[ P , 1] 上均匀分布， 故可令U 2 = ( U 1 - P ) / ( 1 - P ) . 从U 2 产生U 3 的方法是:当I = i 时， U 2 在 [ i / K , ( i + l ) / K ]上均匀分布， 故可令U 3 = K ( U 2 - i / K ) . 这样的做法节省了均匀随机数，增加了一些乘法和除法运算.对F P G A 等并行处理的硬件来说，产生均匀随机数是便宜的，除法运算是耗费的，所以我们不提倡减少均匀数的做法. 而对有C P U 的硬件来说， 减少均匀随机数意味着减少了过程调用，也许是值得的. 再介绍预处理过程.各分点需解下列递推方程求得:从i=1开始求解，直至i = K - 1 .这些方程可事先利用软件求解. 3.2.3 算法性能分析影响随机数产生速度的主要因素之一是f ( x ) 的计算，故把产生每个抽样平均计算f ( x )的次数 ( 计算概率)做为一个性能指标.另外舍选法的平均效率也作为一个性能指标，这个指标反映了每产生一个随机数所需的均匀数个数.产生每个样点X 需计算f ( x ) 的平均概率P f 可利用全概率公式计算:其中10i i iL L L L ---的分母是左边第i 个曲边梯形的下底长，分子是下底与上底的差，这个比值就是在此曲边梯形内投点时计算f ( x ) 的概率.10i i i R R R R ---的意义相仿.舍选法的平均效率” 可利用全概率公式计算:11()()11(1)()()L R KKi i L R A i A i P P K B i K B i η===+-∑∑诸(),(),(),()L L R R A i B i A i B i 分别表示左边各曲边梯形面积、左边各矩形面积、右边各曲边梯形面积和右边各矩形面积.在不同的K 值下，计算了算法用于产生正态分布、 指数分布、 瑞利分布三种标准分布时的上述两个性能参数.各个概率密度函数如下: 正态分布:2())2xf x =-指数分布：()xf x e-=瑞利分布:2()exp()48x xf x =-结果如下图6 :左图反映出概率密度函 数的计算概率P f 随K 的增大而减小， 最终趋于零，例如当K = 1 0 2 4 时， P f 已 非常小;右图反映出 舍选法的平均效率随K 的 增加而提高， 最终趋于 1 ， 也就是三个均匀随机数产生一个预期的随机数.我们可根据实际情况选择合适的K 值.图6 计算概率密度函数的概率3.3正态分布的随机数的产生下面提出了一种已知概率密度函数的分布的随机数的产生方法，以典型的正态分布为例来说名任意分布的随机数的产生方法.如果一个随机数序列服从一维正态分布，那么它有有如下的概率密度函数：22()2()xf xμσ--=其中μ，σ（>0）为常数，它们分别为数学期望和均方差，如果读者对数学期望和均方差的概念还不大清楚，请查阅有关概率论的书.如果取μ =0，σ =0.2,则其曲线为图7 正态分布的概率密度函数曲线从图中可以看出，在μ附近的概率密度大，远离μ的地方概率密度小，我们要产生的随机数要服从这种分布，就是要使产生的随机数在μ附近的概率要大，远离μ处小，怎样保证这一点呢，可以采用如下的方法：在图7的大矩形中随机产生点，这些点是平均分布的，如果产生的点落在概率密度曲线的下方，则认为产生的点是符合要求的，将它们保留，如果在概率密度曲线的上方，则认为这些点不合格，将它们去处.如果随机产生了一大批在整个矩形中均匀分布的点，那么被保留下来的点的横坐标就服从了正态分布.可以设想，由于在μ处的f(x)的值比较大，理所当然的在μ附近的点个数要多，远离μ处的少，这从面积上就可以看出来.我们要产生的随机数就是这里的横坐标.基于以上思想，我们可以用程序实现在一定范围内服从正态分布的随机数.程序如下：double Normal(double x,double miu,double sigma) //概率密度函数{return 1.0/sqrt(2*PI*sigma) *exp(-1*(x-miu)*(x-miu)/(2*sigma*sigma));}double NormalRandom(double miu,double sigma,double min,double max)//产生正态分布随机数{double x;double dScope;double y;do{x = AverageRandom(min,max);y = Normal(dResult, miu, sigma);dScope = AverageRandom(0, Normal(miu,miu,sigma));}while( dScope > y);return x;}参数说明：double m iu：μ，正态函数的数学期望double s igma：σ，正态函数的均方差double m in,double max，表明产生的随机数的范围用如上方法，取μ=0，σ=0.2，范围是-1~1产生400个正态随机数如图8所示：图8 μ=0，σ=0.2，范围在-1~1时的400个正态分布的随机数分布图取μ=0，σ=0.05，范围是-1~1产生400个正态随机数如图9所示：图9 μ=0，σ=0.05，范围在-1~1时的400个正态分布的随机数分布图从图8和图9的比较可以看出，σ越小，产生的随机数靠近μ的数量越多，也说明了产生的随机数靠近μ的概率越大.我们，先产生4000个在0到4之间的正态分布的随机数，取μ=0，σ=0.2，再把产生的数据的数量做个统计，画成曲线，如下图10所示：图10 μ=0，σ=0.2，范围在0~4时的4000个正态分布的随机数统计图从图10中也可以看出，在靠近2处的产生的个数多，远离2处的产生的数量少，该图的轮廓线和概率密度曲线的形状刚好吻合.也就验证了该方法的正确性.有了以上基础，也就用同样的方法，只要知道概率密度函数，也就不难产生任意分布的随机数，方法都是先产生一个点，然后进行取舍，落在概率密度曲线下方的点就满足要求，取其横坐标就是所要获取的随机数参考文献：[1] 肖云茹．概率统计计算方法[M]．天津：南开大学出版社，1994．[2]程兴新.曹敏．统计计算方法EM3．北京：北京大学出版社，1989．[3]王永德等.随机信号分析基础.北京:电子工业出版社，2 0 0 3.[4]皇甫堪等. 现代数字信号处理. 国防科技大学电子科学与工程学院内部印刷，2 0 0 2.。