哈工大概率论课程设计

概率论课程设计随机数的产生摘要：随机数是概率论与数理统计中一个重要的概念。本文研究了随机数的产生，先给出了均匀分布的随机数的产生算法，再通过均匀分布的随机数变换得到其他连续型随机数的产生算法.利用编程给出了产生均匀分布随机数的算法，探讨了同余法的理论原理.通过均匀随机数产生其他分布的随机数，我们列举了几种通用算法，并讨论各个算法的优缺点，最后以正态分布为例验证高效舍选法的优势.关键词：随机数；概率论；均匀分布；算法；目录：一随机数与伪随机数二均匀分布随机数的产生三非均匀分布随机数的产生正文一、随机数与伪随机数随机变量η的抽样序列12,,n，…称为随机数列．如果随机变量η是均匀分布的，则η的抽样序列12,,n，…称为均匀随机数列；如果随机变量η是正态分布的随机变量则称其抽样序列为正态随机数列．比如在掷一枚骰子的随机试验中出现的点数x是一个随机变量，该随机变量就服从离散型均匀分布，x取值为1，2，3，4，5，6，取每个数的概率相等均为1／6．如何得到x的随机数?通过重复进行掷骰子的试验得到的一组观测结果12,,,nxxx就是x的随机数．要产生取值为0，1，2，…，9的离散型均匀分布的随机数，通常的操作方法是把10个完全相同的乒乓球分别标上0，1，2，…，9，然后放在一个不透明的袋中，搅拦均匀后从中摸出一球记号码1x后放回袋中，接着仍将袋中的球搅拌均匀后从袋中再摸出一球记下号码2x后再放回袋中，依次下去，就得到随机序列12,,,nxxx．通常称类似这种摸球的方法产生的随机数为真正的随机数．但是，当我们需要大量的随机数时，这种实际操作方法需要花费大量的时间，通常不能满足模拟试验的需要，比如教师不可能在课堂上做10000次掷硬币的试验，来观察出现正面的频率．计算机可以帮助人们在很短时间产生大量的随机数以满足模拟的需要，那么计算机产生的随机数是用类似摸球方法产生的吗?不是．计算机是用某种数学方法产生的随机数，实际上是按照一定的计算方法得到的一串数，它们具有类似随机数的性质，但是它们是依照确定算法产生的，便不可能是真正的随机数，所以称计算机产生的随机数为伪随机数．在模拟计算中通常使用伪随机数．对这些伪随机数，只要通过统计检验符合一些统计要求，如均匀性、随机性等，就可以作为真正的随机数来使用，我们将称这样产生的伪随机数为随机数．在计算机上用数学方法产生随机数的一般要求如下：1)产生的随机数列要有均匀性、抽样的随机性、试验的独立性和前后的一致性．2)产生的随机数列要有足够长的周期，以满足模拟实际问题的要求．3)产生随机数的速度要快，占用的内存少．计算机产生随机数的方法内容是丰富的，在这里我们介绍几种方法，计算机通常是先产生[0，1]区间上均匀分布的随机数，然后再产生其他分布的随机数．二、均匀分布随机数的产生2.1算法1在vc的环境下，为我们提供了库函数rand()来产生一个随机的整数.该随机数是平均在0~RAND_MAX之间平均分布的，RAND_MAX是一个常量，在VC6.0环境下是这样定义的：#defineRAND_MAX0x7fff它是一个short型数据的最大值，如果要产生一个浮点型的随机数，可以将rand()/1000.0这样就得到一个0~32.767之间平均分布的随机浮点数.如果要使得范围大一点，那么可以通过产生几个随机数的线性组合来实现任意范围内的平均分布的随机数.例如要产生-1000~1000之间的精度为四位小数的平均分布的随机数可以这样来实现.先产生一个0到10000之间的随机整数.方法如下：inta=rand()%10000;然后保留四位小数产生0~1之间的随机小数：doubleb=(double)a/10000.0;然后通过线性组合就可以实现任意范围内的随机数的产生，要实现-1000~1000内的平均分布的随机数可以这样做：doubledValue=(rand()%10000)/10000.0*1000-(rand()%10000)/10000.0*1000;则dValue就是所要的值.但是，上面的式子化简后就变为：doubledValue=(rand()%10000)/10.0-(rand()%10000)/10.0;这样一来，产生的随机数范围是正确的，但是精度不正确了，变成了只有一位正确的小数的随机数了，后面三位的小数都是零，显然不是我们要求的，什么原因呢，又怎么办呢.先找原因，rand()产生的随机数分辨率为32767，两个就是65534，而经过求余后分辨度还要减小为10000，两个就是20000而要求的分辨率为1000*10000*2=20000000，显然远远不够.下面提供的方法可以实现正确的结果：doublea=(rand()%10000)*(rand()%1000)/10000.0;doubleb=(rand()%10000)*(rand()%1000)/10000.0;doubledValue=a-b;则dValue就是所要求的结果.在下面的函数中可以实现产生一个在一个区间之内的平均分布的随机数，精度是4位小数.doubleAverageRandom(doublemin,doublemax){intminInteger=(int)(min*10000);intmaxInteger=(int)(max*10000);intrandInteger=rand()*rand();intdiffInteger=maxInteger-minInteger;intresultInteger=randInteger%diffInteger+minInteger;returnresultInteger/10000.0;}但是有一个值得注意的问题，随机数的产生需要有一个随机的种子，因为用计算机产生的随机数是通过递推的方法得来的，必须有一个初始值，也就是通常所说的随机种子，如果不对随机种子进行初始化，那么计算机有一个缺省的随机种子，这样每次递推的结果就完全相同了，因此需要在每次程序运行时对随机种子进行初始化，在vc中的方法是调用srand（int）这个函数，其参数就是随机种子，但是如果给一个常量，则得到的随机序列就完全相同了，因此可以使用系统的时间来作为随机种子，因为系统时间可以保证它的随机性.2.2算法2：用同余法产生随机数同余法简称为LCG(LinearCongruenceGener-ator)，它是Lehmer于1951年提出来的．同余法利用数论中的同余运算原理产生随机数．同余法是目前发展迅速且使用普遍的方法之一．同余法(LCG)递推公式为1()(mod)nnxaxcm（n=1，2，…)，(1)其中nx，a，c均为正整数．只需给定初值x.，就可以由式(1)得到整数序列{nx}，对每一nx，作变换nu=nx／m，则{nu}(n=1，2，…)就是[0，1)上的一个序列．如果{nu}通过了统计检验，那么就可以将nu作为[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;[()/].nnnnnnnaxcmxaxcaxcmxaxcmaxcm当时，则当时，则这里[x]是取x的整数部分．在电子计算机上由1nx求nx时，可利用整数溢出原理．不进行上面的除法运算．实际上，由于计算机的整数尾部字长为L，机器中可存放的最大整数为2L-1，而此时a1nx+c≥m≥2L-1，因此a1nx+c在机器里存放时占的位数多于L位，于是发生溢出，只能存放nx的右后L位．这个数值恰是模m=2L的剩余，即nx．这就减少了除法运算，而实现了求同余式．经常取模m=2L(L为计算机尾部字长)，正是利用了溢出原理来减少除法运算.由式(1)产生的nx(n=1,2,……)，到一定长度后，会出现周而复始的周期现象，即{nx}可以由其某一子列的重复出现而构成，这种重复出现的子列的最短长度称为序列nx的周期．由式(1)不难看出，{nx}中两个重复数之间的最短距离长度就是它的周期，用T代表周期．周期性表示一种规律性，它与随机性是矛盾的．因此，通常只能取{nx}的一个周期作为可用的随机序列．这样一来，为了产生足够多的随机数，就必须{nx}的周期尽可能地大．由前所述，一般取m=2L，这就是说模m已取到计算机能表示的数的最大数值，意即使产生的随机数列{nx}的周期达到可能的最大数值，如适当地选取参数0x，a，c等，还可能使随机数列{nx}达到满周期.三、非均匀分布随机数的产生3.1一般通用方法3.1.1组合法组合法的基本思想是把预定概率密度函数f(x)表为其它一些概率密度的线性组合.而这些概率密度的随机抽样容易产生.通过这种避难就易的手段我们也许可以达到较高的输出速度和较好的性能.若分布密度函数f(x)能表为如下式(2)所示的函数项级数的和，1()()iiifxpfx(2)其中1iip，诸f(x)皆为概率密度函数.则依如下步骤可产生分布为f(x)一次抽样.(1)产生一个随机自然数I，使I服从如下分布律:P(I=i)=pii=1,2,3……(2)产生服从fI(x)的随机数0X证明利用全概率公式，有:11()()()()()iiiiPxXxdxPIiPxXxdxIipfxdxfxdx故X服从f(x)分布.我们以产生双指数(或拉普拉斯)分布的随机数为例来简单说明这种方法.双指数分布具有概率密度函数f(x)=0.5xef(x)可表为:()0.5()0.5()lrfxfxfx(3)其中()rfx是指数分布，()lfx是指数分布的对称分布.故产生双指数分布的抽样可按如下方法:产生U1,U2~U(0,1);若U10.5，则令X=InU2，否则X=-InU2.在式(2)中，若i→∞，有pi→0，则可用函数列{()}iipfx的前有限项和逼近f(x).这是一种近似的方法，与通常的函数逼近原理相同.只要近似的精度(在某种“精度”的意义之下)达到要求，我们就可以采用近似的方法.使用组合法时，各fi(x)的抽样应该容易产生，故选用合适的概率密度函数族{fi(x)}把任意连续分布表为式(2)，乃是使用组合法的关键.3.3.2概率密度变换法这是一种比较新的通用随机数产生方法.其主要的目的是对一般的f(x)找出较好的覆盖函数以达到较高的效率.我们知道，对某一特定的概率密度f(x)，我们可以使用最优化技术找到好的覆盖函数.但对于一般情况，我们只能期望产生效率尚可的覆盖函数.HORMANN用概率密度变换的方法生成一曲边梯形作为覆盖函数.其原理如下:使用一个变换函数T(x)把预定密度函数f(x)变换为h(x)=T(f(x))，用一个分段线性函数l(x)覆盖h(x)，如图2-4左图;h(x)若是上凸的，则T1(l(x))将是f(x)的一个较好的覆盖函数这个方法在选择合适的T(x)(log(x)或1/xa等)后，能产生随机数包括了较多的分布类型.这个方法有较短的预处理时间，但需要较多的函数计算，不太适合硬件实现.此外，Ahrensl用每段为常数的分段函数作为覆盖函数.Leydold基于ratio-of-uniforms的方法也是一个通用算法.还有一种近似的方法，其产生的随机数与指定分布的随机数具有相同的前四阶矩，但概率分布不一定相同.这里就不详细介绍了.3.2我们的方法当前的通用算法的问题是效率不能任意提高，不够灵活.通常产