Course_4(随机数生成).ppt

1、随机数的生成,郑红星 大连海事大学交通运输管理学院 物流管理与工程系,Management System Simulation,2,2020/8/10,内容,问题的提出 均匀分布随机数发生器及其性质 均匀分布随机数的生成 均匀分布随机数的检验 小结,Management System Simulation,3,2020/8/10,问题的提出,Management System Simulation,4,2020/8/10,管理系统模拟多为随机系统模拟,对随机系统模拟必须得到服从各种分布的随机变量。 问题1: 如何得到这些随机变量生成问题 问题2: 如何证明所得随机变量是所要求的验证问题 研究结

2、果表明: 可以通过对独立同分布的随机数(IID U(0,1)(identical and independent distribution )加以某种转换得到所需的各种随机变量,也就是说可以利用U(0,1)求得所需要的随机数。,问题的提出,Management System Simulation,5,2020/8/10,问题转化为: 问题1: 如何得到这种随机数生成问题 问题2: 如何证明所得随机数是真正的随机数验证问题 因此可以说: 随机变量为管理系统模拟的基础 随机数为得到所需随机变量基础,问题的提出(Continued),Management System Simulation,6,20

3、20/8/10,均匀随机数发生器及其性质,Management System Simulation,7,2020/8/10,随机数的定义： 设:u1,u2, , un是(a, b)上的连续的均匀分布的随机数序列, 取a = 0, b = 1则得(0, 1)区间上均匀分布的随机数,通常表示为U(0,1)。 均匀随机数发生器的主要问题 这种随机数发生器应该具有什么特点； 原有的生成方法及方法存在的问题； 可能的解决方法； 依然存在的问题； 到底需要怎样的随机数发生器。,均匀分布随机数发生器,Management System Simulation,8,2020/8/10,U(0, 1)的统计特性,

4、密度函数 分布函数 期望值 方差,Management System Simulation,9,2020/8/10,U(0, 1)的均匀性和独立性,均匀性： N个样本,均匀落入n等份的(0,1)区间,其期望值为N/n。 独立性： 每一个采样值所落区间与前一个的无关。,Management System Simulation,10,2020/8/10,物理方法: 电子噪声发生器，放射源激励计数器，等。 问题: 物理设备；真随机；但不可重复。 计算机方法： 伪随机数不是真正的随机数 可以重现周期性的重复 统计上满足某些要求统计上满足随机数的数字特性。,如何来得到随机数,Management Sys

5、tem Simulation,11,2020/8/10,不是真正的随机数,而是统计上满足使用需求的随机数,产生方便,使用方便，易于产生。 分布均匀，尽可能接近U(0,1) 统计上独立 产生速度足够快 可以重现 有足够长的周期 用内存尽可能的少 严格意义上来说，现在所有的随机数发生器都不是在概率论意义下的真正随机数，而只能是称为伪随机数，因为无论哪一种随机数发生器都是采用递推算法即有一定规律可言，这已经违反了概率论意义下随机的含义。 但是，如果算法选择合适，则通过该算法得到的数据，可以在一定程度下满足统计检验的要求，具有较好的统计特性，也就可以用于模拟。,我们所真正需要的随机数应该是,Manag

6、ement System Simulation,12,2020/8/10,均匀随机数的生成,Management System Simulation,13,2020/8/10,历史上常用方法介绍,平方取中法 乘法取中法 常数取中法 菲波那契法,Management System Simulation,14,2020/8/10,冯诺依曼于40年代中期提出的。 思路： 取2s位整数做种子； 种子平方后得4s位，如不够4s位，前面补0； 取中间2s位做下一个种子； 归一化后得(0, 1)上的随机数。 通式: Xi+1 = Xi2/10Smod 102S，Ui+1 = Xi+1/102s,平方取中法,

7、Management System Simulation,15,2020/8/10,x0 = 1349,特点与退化问题,平方取中法方法简单，易于实现。 但是会出现退化问题。,Management System Simulation,16,2020/8/10,思路: 取2s位整数做种子； 种子*前次种子后得4s位，如不够4s位，前面补0； 取中间2s位做下一个种子； 归一化后得(0,1)上的随机数。 通式: Xi+2= Xi+1 Xi /10Smod 102S Ui+2=Xi+2/102s,乘法取中法,Management System Simulation,17,2020/8/10,思路 取2

8、s位整数做种子； 种子*某系数后得4s位，如不够4s位，前面补0； 取中间2s位做下一个种子； 规格化后得(0, 1)上的随机数。 通式: Xi+1= k* Xi /10Smod 102S Ui+1= Xi+1/102s,常数取中法,Management System Simulation,18,2020/8/10,思路: 以斐波那契序列运算(0,1,1,2,3,5,0,5,5,2,7,1,0,1,1,.) 通式: Xi+2= Xi+1 +Ximod m Ui+2=Xi+2/m 例：取x0=0; x1=1; m=23=8; 斐波那契序列运算结果 (0,1,1,2,3,5,0,5,5,2,7,1

9、,0,1,1,.) 周期: 3m/2 = 12,菲波那契法,Management System Simulation,19,2020/8/10,历史上常用方法的特点,简单、直观、易于实现； 数据内部特性难以控制，容易出现退化等问题。,Management System Simulation,20,2020/8/10,目前常用的方法线性同余法,混和线性同余法 乘法同余法 二次同余法 加法同余法,Management System Simulation,21,2020/8/10,通式: Xi+1 = aXi + C mod m Ui+1 = Xi+1 /m 其中: a乘子；c增量；x0初值；m模

10、通常4个参数都取整数 参数的选择不同对产生的随机序列和统计特性有极大的影响. 关键问题：如何选择参数以得到好的结果。 Xi的范围 参数X0，a，c，m。,混和线性同余法,Management System Simulation,22,2020/8/10,混和线性同余法的性质简析,显然，根据通式Xi+1 = aXi + C mod m，可得： 0 Xi m 1； 由于Xi为整数，因此由此所得到的随机数肯定是有限个数； 且Xi肯定是不随机的，因为可以根据递归调用发现，只要m，a，C，X0一旦确定，则Xi是函数可得的，本质上不存在任何随机性。,Management System Simulation

11、,23,2020/8/10,通式: Xi+1 = aXi + C mod m = (aXi + C)(aXi + C)/m m 令 N = aXi + c 则: Xi+1 = N - N/mm Xi的范围 0 Xi+1 m； Ui+1 = Xi+1/m 1；Xi为整数，所以max (Xi+1) = m-1，0 = Xi+1 = m-1 而且, 其数据周期的重复出现； 因为Xi取值有限，因此肯定会在m次迭代内出现重复值； 又因为Xi+1取值，只与Xi相关，则一旦Xi重复出现，后续Xi+1都会完全重复。,Xi的范围,Management System Simulation,24,2020/8/10

12、,Xi+1为整数，所以只能为0,1,2,.,m-1中的一个；则Ui+1也为有限长, 为0,1/m,2/m,.,m-1/m之一。 带来的问题：可能会不均匀。即最大的周期长度可能是m，但是并不一定。 m越大，Xi+1的选择范围越大,数字出现的频率越均匀, Ui+1的可选择范围也越大。反之如果m取的小，则Ui+1很难为均匀分布。 通常,取m尽可能的大。如:取 m = 2b其中b为计算机最大位长。如果取b=31, m 2*109, 即该随机数发生器能够产生21亿以上个不重复随机数. Ui+1为周期式的重复，在每个m段中, 计算公式是相同的.,几点说明,Management System Simulat

13、ion,25,2020/8/10,适当选择X0, a, c, m参数, 以便获得满周期的伪随机数 Xn完全由4个参数确定，一般4个参数取整数, 还应有 m 0; a m; c m ; X0 m 对X0, a, c, m的选择与Xi+1的结果紧密相关; 如下: X1 = (aX0+c) mod m X2 = (aX1+c) mod m = a(aX0+c)+c mod m = a2X0+c(a+1)mod m X3 = (aX2+c) mod m = a(a2X0+c(a+1)+cmod m = a3X0+c (a2+a+1) mod m . Xn = (aXn-1 + c) mod m = a

14、nX0 + c*(an-1)/(a-1) mod m 即： Xn完全由a, c, X0, m四个参数确定,满周期的伪随机数,Management System Simulation,26,2020/8/10,定理：当且仅当满足下列条件时，线性同余法是满周期的，即p=m： c与m互为质数，即能同时整除c和m的唯一正整数为1； 如q为整除m的质数，则a-1也能被q整除，即m/q 为整数，(a-1)/q也为整数； 若m可被4整除，则 a-1也可被4整除，即m/4为整数，(a-1)/4也为整数。 当a=zr+1，其中z是计算机中用于表示数字的基数，m=zb时（b为计算机字长），且有z=r=b。,获取满

15、周期随机数的判定定理,Management System Simulation,27,2020/8/10,例：取m=16，a=5，c=3，X0=7。通式为: Xi+1 =(5Xi+3) mod 16 则c=3和m=16是互质的。 令q=2，m/2为整数，(a-1)/2=4/2也为整数 m/4=16/4为整数，(a-1)/4=4/4也为整数 即这组参数可以生成满周期的随机序列。 计算结果如下，显然周期为16。,获取满周期随机数示例,Management System Simulation,28,2020/8/10,m的选取 越大越好，一般取m=2b，b为计算机字长(32位)，产生随机数总数为2b

16、-1 21亿。m可取一个足够大的质数。 a和c的选取 c与m互质，一般为奇数；a-1可为4整除；则a一定为4的倍数+1。 x0的选取 因为xi可以是0m之间的任一数，所以x0的选取不是很重要，但是如果取x0=0有时确实会使通式的结果退化。,一些惯例,Management System Simulation,29,2020/8/10,由Ui+1=Xi+1/m可知，当m较大时，每得到一个 Ui+1都要进行一次大数的除法。有人利用计算机位数的截断的特点，简化除大数的问题，提高计算机计算的效率。 称之为“整数溢出法”。 思路: 计算机的寄存器本身就是一个MOD函数, 如果将一个超出字长的整数装入, 多

17、出的部分会自动溢出, 实现实际上的取模运算, 从而简化了计算量. 例：假想有一台计算机，字长Xi+1= (5Xi+3) mod 16 当5Xi+3=16时，二进制为10000, Xi+1=0 当5Xi+3=1016时，二进制为11100, Xi+1=12,计算Ui+1(除大数的效率),Management System Simulation,30,2020/8/10,由Coveyou和Macpherson 给出的b = 35时的随机数发生器： Xi+1=(515Xi + 1) mod 235 该通式满足三个条件, 为满周期的随机数发生器. 由Kobayashi提出, Knuth介绍的b = 3

18、1时的随机数发生器 Xi+1=(314159269Xi + 453806245) mod 231 该通式也满足三个条件, 为满周期的随机数发生器. 是否是真正的随机数？ 公式计算，有通式： Xn = an X0 + c*( an -1)/(a-1) mod m 伪随机数。,例子及说明,Management System Simulation,31,2020/8/10,通式：(早期计算机中经常采用) Xi+1 =aXi mod m，Ui+1 =Xi+1 /m 特点: c=0 不满足前定理，可能得到满周期的随机序列。 如: D. H. Lehmer法可以求得周期为p=m-1随机序列。 如: Ban

19、ks等人的方法可以求最大可能周期的随机序列。 条件：取x0为奇数,且满足a = 3 + 8k(k为整数)； 当m=2b和c=0时, 可以获得最大可能周期p=m/4=2b-2, 也可以使用“整数溢出法”。,乘法同余法,Management System Simulation,32,2020/8/10,例：a=13, b=6, m=2b=26=64, x0=1,2,3,4。 结果：当x0=1, 3时, 可以获得最大可能周期为p=26-2 =16随机序列。而x0=2,4时,得不到最大可能周期。,乘法同余法示例,Management System Simulation,33,2020/8/10,二次同

20、余法 通式：Xi+1 = aXi2 +bXi +cmod m，Ui+1 =Xi+1 /m 特点 a和b 为互质的常数, c为一奇整数,且不为3和5除尽, 并与m互质。 m的取法和线性同余法相同. 加法同余法 通式：Xi = Xi-1+Xi-n mod m，Ui=Xi /m 特点 免除了乘法运算, 所以计算速度较快。 但事先需要若干个数, 从已知的x1,x2,.,xi计算xi+1,二次同余法和加法同余法,Management System Simulation,34,2020/8/10,均匀分布随机数的检验,Management System Simulation,35,2020/8/10,随机

21、数生成 还可以以许多其他方法生成随机数 特点：有周期性，伪随机 应用中，我们并不真正需要“真”的随机数； 满足“某些条件”的“伪”随机数已经足够了 面临的问题 如何验证所得的随机数是否满足使用要求。 由于“伪”随机数应该是随机数，又具有均匀、独立的特点，因此要进行以下三个方面的检验： 参数检验; 均匀性检验; 独立性检验。,检验均匀生成随机数的背景,Management System Simulation,36,2020/8/10,目的：检验样本在基本参数检验方面是否符合要求。也可视为是最为基本的统计上的检验。 对于随机数u1, u2, ., uN (N为数据长度)，一些常见的进行检验的参数包

22、括： 均值的估计值: u = ui / N； 方差的估计值: S2 =(1/(N-1) (u-ui)2； 均值估计值的期望: E(u)=1/2； 均值估计值的方差:V(u)=1/12N； 方差估计值的期望: E(S2)=1/12； 方差估计值的方差: V(S2)=1/180N。,参数检验,Management System Simulation,37,2020/8/10,可构造如下统计量： 注意到当N足够大时，V1，V2近似服从N(0, 1)分布，对给定的显著水平(假定为0.05)，查表可得临界值为1.96。 当|V1| 1.96时，拒绝E(u)=1/2和V(u)=1/12N的假设； 当|V2

23、| 1.96时，拒绝E(s2)=1/12和V(s2)=1/180N假设。 否则接受这两个假设. 特点：转化为标准正态分布的检验,统计量构造,Management System Simulation,38,2020/8/10,例：用excel产生100个随机数进行参数估计 均值估计：=average(A1:A100)； 方差估计：=stdev(A1:A100)；或者=var(A1:A100) 随机数发生器：=RAND( ) 检验：取显著水平=0.05；N=100；临界值为1.96 计算结果；该批数据通过参数检验。,参数检验示例,Management System Simulation,39,20

24、20/8/10,对于给定的随机数，其均匀性检验主要是采用频率检验法，手段是检验所产生的随机数落在各子区间的频率和理论频率之间的差异是否显著。常用的频率检验法有两种： X2检验 K-S检验(Kolmogorov-Smirnov),均匀性检验,Management System Simulation,40,2020/8/10,思路：将(0,1 )区间分成k个相等的子区间，检查总数为N的随机数(样本)落在各子区间的实际频数ni与理论频数的差异是否显著。 设落在第i个子区间的随机数的实际频数为ni，落在第i个子区间的理论频数为mi，mi = N / k。建立X2统计量 当样本量足够大时，该统计量近似服

25、从自由度为k-1的 X2分布，显著水平一般取0.010.05。 当统计量 临界值 时，则在该显著水平上拒绝该批样本为均匀的假设。,X2检验,Management System Simulation,41,2020/8/10,分区间k； 统计ni和计算mi； 计算统计量X2； 根据显著水平及自由度k-1从X2表查临界值； 比较X2和临界值, 确定是否接受。 特点为必须在大样本时才有效 本方法的问题：区间选择对统计量的影响。,X2检验步骤,Management System Simulation,42,2020/8/10,区间及频率：,X2检验示例,Management System Simula

26、tion,43,2020/8/10,K-S检验是将U(0,1)连续的分布函数F(x)同由观察而得到达到经验分布函数SN(x)进行比较。 连续的U(0,1)的分布函数的定义为: F(x) = x (0=x=1) 经验分布函数的定义为: SN(x) = nx / N，nx为=x的u的个数。 经验分布函数为（xi , SN(x)）阶梯状。 特点：可以小样本，但是通用性差。,K-S检验,Management System Simulation,44,2020/8/10,将样本从小到大排列； 计算D+和D-，其中： D+ = maxi/N - ui； D- = maxui - (i-1)/N (1=i=

27、N) (绝对值)； 计算样本统计值 D=maxD+, D-； 根据显著水平及样本数N查K-S表, 确定临界值 D； 将临界值同统计值进行比较，确定是否接受。,K-S检验步骤,Management System Simulation,45,2020/8/10,给定5个生成的随机数：0.44,0.81,0.14,0.05,0.93,K-S检验示例,Management System Simulation,46,2020/8/10,K-S检验示例,Management System Simulation,47,2020/8/10,检查随机数样本间是否存在相关性。常用的有下列三种方法： 相关性检验；

28、组合规律性检验； 连贯性检验。,独立性检验,Management System Simulation,48,2020/8/10,设给定随机数列u1,u2,.,uN，前后相距为k的样本的相关系数的估计值为k。则 若k阶自相关系数为0，则k=0。建立统计量： 当N充分大时，该统计量近似服从N(0,1)分布； 给定显著水平 = 0.05，计算出k=1,2,.,k,时各阶自相关系数的估计值； 若 |Vk|1.96，则拒绝相关系数k =0的假设。见Excel示例。,相关性检验,Management System Simulation,49,2020/8/10,用于检查随机数样本中重复数字出现的频率，从而检查该组样本的组合规律性。如，三位数： 三个数字相同出现的概率为: P(三个数字相同)=0.1*0.1=0.01 三个数字不同出现的概率为: P(三个数字不同)=0.9*0.8=0.72 有两个数字同时出现的概率为: P(两个数字相同)=C(3, 2)(0.1*0.9)=1-0.01-0.72=0.27 以此作为理论频率与实际频数作比较，然后对其进行X2检验。在给定的显著水平下，找出临界值，根据临界值确定差距是否显著。,组合规律性检验-Poker检验,Management System Simulation,50,2020/8/10,Poker检