哈工大概率论课程设计

1、概率论课程设计随机数的产生 摘要：随机数是概率论与数理统计中一个重要的概念。本文研究了随机数的产生，先给出了均匀分布的随机数的产生算法，再通过均匀分布的随机数变换得到其他连续型随机数的产生算法.利用编程给出了产生均匀分布随机数的算法，探讨了同余法的理论原理.通过均匀随机数产生其他分布的随机数，我们列举了几种通用算法，并讨论各个算法的优缺点，最后以正态分布为例验证高效舍选法的优势.关键词：随机数；概率论；均匀分布；算法；目录：一 随机数与伪随机数二 均匀分布随机数的产生三 非均匀分布随机数的产生正文一、 随机数与伪随机数随机变量的抽样序列，称为随机数列如果随机变量是均匀分布的，则的抽样序列，称为

2、均匀随机数列；如果随机变量是正态分布的随机变量则称其抽样序列为正态随机数列比如在掷一枚骰子的随机试验中出现的点数x是一个随机变量，该随机变量就服从离散型均匀分布，x取值为1，2，3，4，5，6，取每个数的概率相等均为16如何得到x的随机数?通过重复进行掷骰子的试验得到的一组观测结果就是x的随机数要产生取值为0，1，2，9的离散型均匀分布的随机数，通常的操作方法是把10个完全相同的乒乓球分别标上0，1，2，9，然后放在一个不透明的袋中，搅拦均匀后从中摸出一球记号码后放回袋中，接着仍将袋中的球搅拌均匀后从袋中再摸出一球记下号码后再放回袋中，依次下去，就得到随机序列通常称类似这种摸球的方法产生的随机

3、数为真正的随机数但是，当我们需要大量的随机数时，这种实际操作方法需要花费大量的时间，通常不能满足模拟试验的需要，比如教师不可能在课堂上做10000次掷硬币的试验，来观察出现正面的频率计算机可以帮助人们在很短时间产生大量的随机数以满足模拟的需要，那么计算机产生的随机数是用类似摸球方法产生的吗?不是计算机是用某种数学方法产生的随机数，实际上是按照一定的计算方法得到的一串数，它们具有类似随机数的性质，但是它们是依照确定算法产生的，便不可能是真正的随机数，所以称计算机产生的随机数为伪随机数在模拟计算中通常使用伪随机数对这些伪随机数，只要通过统计检验符合一些统计要求，如均匀性、随机性等，就可以作为真正的

4、随机数来使用，我们将称这样产生的伪随机数为随机数在计算机上用数学方法产生随机数的一般要求如下：1)产生的随机数列要有均匀性、抽样的随机性、试验的独立性和前后的一致性2)产生的随机数列要有足够长的周期，以满足模拟实际问题的要求3)产生随机数的速度要快，占用的内存少计算机产生随机数的方法内容是丰富的，在这里我们介绍几种方法，计算机通常是先产生0，1区间上均匀分布的随机数，然后再产生其他分布的随机数二、均匀分布随机数的产生2.1 算法1    在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数.该随机数是平均在0RAND_MAX之间平均分布的，RAND_

5、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

6、/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;   这样一来，产生的随机数范围是正确的，但是精度不正确了，变成了只有一位正确的小数的随机数了，后面三位的小数都是零，显然不是我们要求的，

7、什么原因呢，又怎么办呢.   先找原因，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就是所要求的结果.在下面的函数中可以实现产生一个在一

8、个区间之内的平均分布的随机数，精度是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.

9、0;    但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随机种子，如果不对随机种子进行初始化，那么计算机有一个缺省的随机种子，这样每次递推的结果就完全相同了，因此需要在每次程序运行时对随机种子进行初始化，在vc中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时间来作为随机种子，因为系统时间可以保证它的随机性.   2.2 算法2：用同余法产生随机数同余法简称为LCG(Line

10、ar Congruence Gener-ator)，它是Lehmer于1951年提出来的同余法利用数论中的同余运算原理产生随机数同余法是目前发展迅速且使用普遍的方法之一同余法(LCG)递推公式为 （n=1，2，)， (1)其中，a，c均为正整数只需给定初值x.，就可以由式(1)得到整数序列，对每一，作变换=m，则(n=1，2，)就是0，1)上的一个序列如果通过了统计检验，那么就可以将作为0，1)上的均匀分布随机数在式(1)中，若c=0，则称相应的算法为乘同余法，并称口为乘子；若c0，则称相应的算法为混合同余法同余法也称为同余发生器，其中称为种子由式(1)可以看出，对于十进制数，当取模m=10(

11、k为正整数)时，求其同余式运算较简便例如36=31236(mod10)，只要对21236从右截取k=2位数，即得余数36同理，对于二进制数，取模m=2时，求其同余式运算更简便了电子计算机通常是以二进制形式表示数的在整数尾部字长为L位的二进制计算机上，按式(1)求以m为模的同余式时，可以利用计算机具有的整数溢出功能设L为计算机的整数尾部字长，取模m=2，若按式(1)求同余式时，显然有这里x是取x的整数部分在电子计算机上由求时，可利用整数溢出原理不进行上面的除法运算实际上，由于计算机的整数尾部字长为L，机器中可存放的最大整数为2-1，而此时a+cm2-1，因此a+c在机器里存放时占的位数多于L位，

12、于是发生溢出，只能存放的右后L位这个数值恰是模m=2的剩余，即这就减少了除法运算，而实现了求同余式经常取模m=2(L为计算机尾部字长)，正是利用了溢出原理来减少除法运算.由式(1)产生的(n=1,2,)，到一定长度后，会出现周而复始的周期现象，即 可以由其某一子列的重复出现而构成，这种重复出现的子列的最短长度称为序列的周期由式(1)不难看出，中两个重复数之间的最短距离长度就是它的周期，用T代表周期周期性表示一种规律性，它与随机性是矛盾的因此，通常只能取的一个周期作为可用的随机序列这样一来，为了产生足够多的随机数，就必须的周期尽可能地大由前所述，一般取m=2，这就是说模m已取到计算机能表示的数的

13、最大数值，意即使产生的随机数列的周期达到可能的最大数值，如适当地选取参数，a，c等，还可能使随机数列达到满周期.三、非均匀分布随机数的产生3.1 一般通用方法3.1.1组合法组合法的基本思想是把预定概率密度函数f ( x ) 表为其它一些概率密度的线性组合.而这些概率密度的随机抽样容易产生.通过这种避难就易的手段我们也许可以达到较高的输出速度和较好的性能.若分布密度函数f ( x ) 能表为如下式(2)所示的函数项级数的和， (2)其 中，诸f( x )皆为概率密度函数.则依如下步骤可产生分布为f ( x )一次抽样.( 1 ) 产生一个随机自然数I ， 使I服从如下分布律:P ( I = i

14、 ) = p i = 1 , 2 , 3( 2 ) 产生服从f( x )的随机数证明利用全概率公式，有:故X服从f ( x ) 分布.我们以产生双指数(或拉普拉斯)分布的随机数为例来简单说明这种方法.双指数分布具有 概率密度函数f ( x ) = 0 . 5 f ( x ) 可表为: (3)其中是指数分布，是指数分布的对称分布.故产生双指数分布的抽样可按如下方法: 产生U , UU ( 0 , 1 ) ;若U > 0 . 5 ， 则令X = I n U，否则X = - I n U. 在式(2) 中， 若i， 有p 0 ，则可用函数列的前有限项和逼近f ( x ).这是一种近似的方法，与通

15、常的函数逼近原理相同.只要近似的精度 ( 在某种“精度”的意义之下) 达到要求，我们就可以采用近似的方法 .使用组合法时，各f( x ) 的抽样应该容易产生，故选用合适的概率密度函数族 f( x )把任意连续分布表为式(2) ，乃是使用组合法的关键.3.3.2 概率密度变换法 这是一种比较新的通用随机数产生方法.其主要的目的是对一般的f(x)找出较好的覆盖函数以达到较高的效率.我们知道，对某一特定的概率密度f(x)，我们可以使用最优化技术找到好的覆盖函数.但对于一般情况，我们只能期望产生效率尚可的覆盖函数. H O R M A N N用概率密度变换的方法生成一曲边梯形作为覆盖函数.其原理如下:

16、使用一个变换函数T (x)把预定密度函数f ( x ) 变换为h ( x ) = T ( f ( x ) ) ，用一个分段线性函数l ( x )覆盖h ( x )，如图2 - 4 左图; h ( x ) 若是上凸的，则T( l ( x ) )将是f ( x ) 的一个较好的覆盖函数 这个方法在选择合适的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

17、f - u n i f o r m s 的方法也是一个通用算法.还有一种近似的方法，其产生的随机数与指定分布的随机数具有相同的前四阶矩，但概率分布不一定相同.这里就不详细介绍了.3.2 我们的方法 当前的通用算法的问题是效率不能任意提高，不够灵活. 通常产生每个所需随机数X需以较大的概率计算f ( x )等函数.我们认为在速度要求非常高的场合，计算f ( x )是不利的，尤其以硬件进行函数计算是十分不利的.针对己有通用算法的不足，我们提出了基于组合法的通用算法.主要目的是尽可能地减少三角、指数、对数等超越函数的计算，以便硬件实现.产生任意连续分布随机数的高效舍选算法 本文提出一种通用算法，可视

18、需要使效率接近1 , 而且f ( x ) 的计算概率可任意小. 这些优点的取得是以长的预处理时间为代价的.在需要产生大量随机样本的场合 ( 例如通信系统的误码率测试， 可能需要数小时乃至数天的仿真时间) ， 本算法将有很大的优势， 尽管有看法认为只有能用简单代码实现的算法才会被经常使用.3.2.1 算法原理 假定预定的连续概率密度函数f ( x ) 为单峰的( 这是实际的大多数情况) ， 已知其峰值点为m .一般f ( x ) 不关于x =m 对称，如图2 -5 .我们假定f ( x ) 定义在有界的区间 a , b 上( 上文说过， 对正态分布这类定义区间无限的情况， 我们把这个区间取得足够

19、大就可以了) . 直线X=m把f ( x ) 曲线与X轴所围面积分为左右两部分， 我们把左右两部分各等分为K份，一共得到2 K个曲边梯形.并用2 K个矩形各自覆盖相应的曲边梯形. 我们的想法是利用舍选法的几何意义，分别在上述 2 K个曲边梯形内均匀投点，从而使随机点在f ( x ) 曲线与x 轴所围的整个区域中均匀分布，这样即可产生f ( x ) 的抽样X . 而在曲边梯形内均匀投点可使用简单舍选法:先在各个矩形内均匀投点，再选出落于相应曲 边梯形内的点. 这种投点法浪费的点只位于各个矩形的一角， 显然效率大大高于简单舍选法.最为重要的是:随着K的增大，效率会不断提高.另外，只有当投点位于曲边

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

21、一腰为曲边的梯形形状的概率密度函数.依如下步骤可产生分布为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 .证明利用全概率公式，有:故x服从 f ( x ) 分布.下面完整地描述这个方法: S t e p l( 产生J ) : S t

22、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 表示不大于x的最大整数. 产生 ( x ) 的样本需区别j = 1 与j = 2 两种情况. 图2 - 6 示出j = 2 时一 典型的 ( x ) ， 用简单舍选法产生其抽样，覆盖函数为矩形. 首先产生一个 0 , R 的均匀数， 如它属于 0 , R 小无需再产生y 轴方向的均匀随机数，接受此均匀数

23、即可;否则还需产生一个Y 轴方向的均匀随机数进行投点，那些落在曲边下方的点被接受，投在矩形右上角的点被舍弃.同理易得j = 1 时的产生法. 整个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 - L ) U 3 + L : i f W> 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

24、p ; e l s e l o o p : 产生 0 , 1 上的 均匀随机数U 3 , W = ( R - R o ) U 3 + R o ; i f W< R，返回 X = W; e l s e 产生 0 , 1 上的均匀随机数V ; i f f ( W) - f ( R) C ( f ( R ) - f ( R ) ) 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的硬件来说， 减少均匀随机数意味着减少