随机数的产生

1、随机数的产生摘要本文研究了连续型随机数列的产生，先给出了均匀分布的随机数的产生算法，在通过均匀分布的随机数变换得到其他连续型随机数的产生算法.在vc环境下，我们给出了产生均匀分布随机数的算法，然后探讨了同余法的理论原理.通过均匀随机数产生其他分布的随机数，我们列举了几种通用算法，并讨论各个算法的优缺点，最后以正态分布为例验证高效舍选法的优势.正文一、 随机数与伪随机数随机变量的抽样序列，称为随机数列如果随机变量是均匀分布的，则的抽样序列，称为均匀随机数列；如果随机变量是正态分布的随机变量则称其抽样序列为正态随机数列比如在掷一枚骰子的随机试验中出现的点数x是一个随机变量，该随机变量就服从离散型均

2、匀分布，x取值为1，2，3，4，5，6，取每个数的概率相等均为16如何得到x的随机数?通过重复进行掷骰子的试验得到的一组观测结果就是x的随机数要产生取值为0，1，2，9的离散型均匀分布的随机数，通常的操作方法是把10个完全相同的乒乓球分别标上0，1，2，9，然后放在一个不透明的袋中，搅拦均匀后从中摸出一球记号码后放回袋中，接着仍将袋中的球搅拌均匀后从袋中再摸出一球记下号码后再放回袋中，依次下去，就得到随机序列通常称类似这种摸球的方法产生的随机数为真正的随机数但是，当我们需要大量的随机数时，这种实际操作方法需要花费大量的时间，通常不能满足模拟试验的需要，比如教师不可能在课堂上做10000次掷硬币

3、的试验，来观察出现正面的频率计算机可以帮助人们在很短时间产生大量的随机数以满足模拟的需要，那么计算机产生的随机数是用类似摸球方法产生的吗?不是计算机是用某种数学方法产生的随机数，实际上是按照一定的计算方法得到的一串数，它们具有类似随机数的性质，但是它们是依照确定算法产生的，便不可能是真正的随机数，所以称计算机产生的随机数为伪随机数在模拟计算中通常使用伪随机数对这些伪随机数，只要通过统计检验符合一些统计要求，如均匀性、随机性等，就可以作为真正的随机数来使用，我们将称这样产生的伪随机数为随机数在计算机上用数学方法产生随机数的一般要求如下：1)产生的随机数列要有均匀性、抽样的随机性、试验的独立性和前

4、后的一致性2)产生的随机数列要有足够长的周期，以满足模拟实际问题的要求3)产生随机数的速度要快，占用的内存少计算机产生随机数的方法内容是丰富的，在这里我们介绍几种方法，计算机通常是先产生0，1区间上均匀分布的随机数，然后再产生其他分布的随机数二、均匀分布随机数的产生2.1 算法1 在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数.该随机数是平均在0RAND_MAX之间平均分布的，RAND_MAX是一个常量，在VC6.0环境下是这样定义的：#define RAND_MAX 0x7fff它是一个short 型数据的最大值，如果要产生一个浮点型的随机数，可以将rand()/1000

5、.0这样就得到一个032.767之间平均分布的随机浮点数.如果要使得范围大一点，那么可以通过产生几个随机数的线性组合来实现任意范围内的平均分布的随机数.例如要产生-10001000之间的精度为四位小数的平均分布的随机数可以这样来实现.先产生一个0到10000之间的随机整数.方法如下 ： int a = rand()%10000;然后保留四位小数产生01之间的随机小数：double b = (double)a/10000.0;然后通过线性组合就可以实现任意范围内的随机数的产生，要实现-10001000内的平均分布的随机数可以这样做：double dValue = (rand()%10000)/1

6、0000.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=2000000

7、0，显然远远不够.下面提供的方法可以实现正确的结果：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 =

8、(int)(max*10000); int randInteger = rand()*rand(); int diffInteger = maxInteger - minInteger; int resultInteger = randInteger % diffInteger + minInteger; return resultInteger/10000.0; 但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随机种子，如果不对随机种子进行初始化，那么计算机有一个缺省的随机种子，这样每次递推的结果就

9、完全相同了，因此需要在每次程序运行时对随机种子进行初始化，在vc中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时间来作为随机种子，因为系统时间可以保证它的随机性. 调用方法是srand(GetTickCount(),但是又不能在每次调用rand()的时候都用srand(GetTickCount()来初始化，因为现在计算机运行时间比较快，当连续调用rand()时，系统的时间还没有更新，所以得到的随机种子在一段时间内是完全相同的，因此一般只在进行一次大批随机数产生之前进行一次随机种子的初始化.下面的代码产生了40

10、0个在-11之间的平均分布的随机数.double dValue400; srand(GetTickCount(); for(int i= 0;i 0 . 5 ， 则令X = I n U，否则X = - I n U. 在式(2) 中， 若i， 有p 0 ，则可用函数列的前有限项和逼近f ( x ).这是一种近似的方法，与通常的函数逼近原理相同.只要近似的精度 ( 在某种“精度”的意义之下) 达到要求，我们就可以采用近似的方法 .使用组合法时，各f( x ) 的抽样应该容易产生，故选用合适的概率密度函数族 f( x )把任意连续分布表为式(2) ，乃是使用组合法的关键.概率密度变换法 这是一种比较

11、新的通用随机数产生方法.其主要的目的是对一般的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( l ( x ) )将是f ( x ) 的一个较好的覆盖函数，图3

12、概率密度变换法原理图 这个方法在选择合适的T ( x ) ( l o g ( x ) 或1 / x等) 后，能产生随机数包括了较多的分布类型.这个方法有较短的预处理时间，但需要较多的函数计算，不太适合硬件实现. 此外，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 我们的方法 当前的通用算法的问题是效率不能任意提高，不够灵活. 通常产生每个

13、所需随机数X需以较大的概率计算f ( x )等函数.我们认为在速度要求非常高的场合，计算f ( x )是不利的，尤其以硬件进行函数计算是十分不利的.针对己有通用算法的不足，我们提出了基于组合法的通用算法.主要目的是尽可能地减少三角、指数、对数等超越函数的计算，以便硬件实现.产生任意连续分布随机数的高效舍选算法 本文提出一种通用算法，可视需要使效率接近1 , 而且f ( x ) 的计算概率可任意小. 这些优点的取得是以长的预处理时间为代价的.在需要产生大量随机样本的场合 ( 例如通信系统的误码率测试， 可能需要数小时乃至数天的仿真时间) ， 本算法将有很大的优势， 尽管有看法认为只有能用简单代码

14、实现的算法才会被经常使用.算法原理 假定预定的连续概率密度函数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 , . . .

15、 , K -1 )是分点.并令R=L =m, L = a，R= b .图4 均分f ( x ) 曲 线与X 轴所围面积 我们的想法是利用舍选法的几何意义，分别在上述 2 K个曲边梯形内均匀投点，从而使随机点在f ( x ) 曲线与x 轴所围的整个区域中均匀分布，这样即可产生f ( x ) 的抽样X . 而在曲边梯形内均匀投点可使用简单舍选法:先在各个矩形内均匀投点，再选出落于相应曲 边梯形内的点. 这种投点法浪费的点只位于各个矩形的一角， 显然效率大大高于简单舍选法.最为重要的是:随着K的增大，效率会不断提高.另外，只有当投点位于曲边梯形的曲边之下时， 才需计算f ( x ) ，而且计算f (

16、 x ) 概率是随着K的增加而减小的. 我们每次“ 按概率”随机选中一个曲 边梯形进行投点. 这需要两步完成:先选择左边还是右边，再于此边的K个曲边梯形中选择一个.这里的概率显然就是面积，这可以从以下的推导中看出来.为清晰起见，我们先阐述随机数的产生法，而把面积的均分这个预处理过程置于随后.算法推导令为左边面积.则左边各曲边梯形面积皆为 P / K，右边各曲边梯形面积皆为( ( I -P ) / K . f ( x ) 可表为: (4)诸 ( x ) ( j = 1 , 2 ; i = 1 , 2 . . . k ) 皆为一腰为曲边的梯形形状的概率密度函数.依如下步骤可产生分布为f ( x )

17、 的一次抽样: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 .证明利用全概率公式，有:故x服从 f ( x ) 分布.下面完整地描述这个方法: S t e p l( 产生J ) : S t e p l . l产生 0 , 1 上的均匀随机数U 1 ; S t e

18、 p 1 . 2若U 1 L ，返回 X = W; e l s e 产生 O ， l 上的均匀随机数V ; i f f ( W) - f ( L) ( f ( L ) - f ( L ) ) 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 - R o ) U 3 + R o ; i f W0）为常数，它们分别为数学期望和均方差，如果读者对数学期望和均方差的概念还不大清楚，请查阅有关概率论的书.如果取 =0， =0.2,则其曲线为 图7 正态分布的概率密度函数曲线 从图中可以看出

19、，在附近的概率密度大，远离的地方概率密度小，我们要产生的随机数要服从这种分布，就是要使产生的随机数在附近的概率要大，远离处小，怎样保证这一点呢，可以采用如下的方法：在图7的大矩形中随机产生点，这些点是平均分布的，如果产生的点落在概率密度曲线的下方，则认为产生的点是符合要求的，将它们保留，如果在概率密度曲线的上方，则认为这些点不合格，将它们去处.如果随机产生了一大批在整个矩形中均匀分布的点，那么被保留下来的点的横坐标就服从了正态分布.可以设想，由于在处的f(x)的值比较大，理所当然的在附近的点个数要多，远离处的少，这从面积上就可以看出来.我们要产生的随机数就是这里的横坐标. 基于以上思想，我们可

20、以用程序实现在一定范围内服从正态分布的随机数.程序如下：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; 参数说明：d