蒙特卡罗随机数.ppt

随机数的定义及产生方法 伪随机数产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性,2.2 MC随机数,随机数的定义及性质随机数表物理方法,1.MC随机数的定义及产生方法,什么是随机数？,单个的数字不是随机数,是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；,在一个均匀分布的随机数中，每一个体出现的概率是均等的；,例如：在0,1区间上均匀分布的随机数序列中，0.00001与0.5出现的机会均等,1.1 MC随机数定义及特性,由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。,在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为0，1上的均匀分布，其分布密度函数为：分布函数为：,1.1 MC随机数定义及特性,由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号表示。由随机数序列的定义可知，1，2，是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（n+1，n+2，n+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai，如下等式成立：,1.1 MC随机数定义及特性,其中P()表示事件发生的概率。反之，如果随机变量序列1,2对于任意自然数s，由s个元素所组成的s维空间上的点（n+1，n+s）在Gs上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。,1.1 MC随机数定义及特性,为了产生随机数，可以使用随机数表。随机数表是由0，1，9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。,1.2 MC随机数表,用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：其中i（i=1，2，m）或者为0，或者为1。,1.3 MC随机数物理方法,因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。,1.3 MC随机数物理方法,0,1区间上均匀分布的随机数是Monte Carlo模拟的基础:,0,1均匀分布的随机数的产生方法：,利用一些具有内在的随机性的过程,缺点：模拟的结果不可再现，使得模拟程序的找错困难,利用事先制订好的随机数表,缺点：表的容量有限，不适合需要大量随机数的应用,服从任意分布的随机数序列可以用0,1区间均匀分布的随机数序列作适当的变换或舍选后求得,1.3 MC随机数小结,伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量,2 伪随机数,在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：产生随机数序列。对于给定的初始值1,2，k，确定n+k，=1，2，。经常使用的是k=1的情况，其递推公式为：对于给定的初始值1，确定n+1，=,2 伪随机数,用数学方法产生的随机数，存在两个问题：递推公式和初始值1,2，k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的0，1上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n，n(n n)，使得下面等式成立：随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。,2 伪随机数两个问题,由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。,2 伪随机数两个问题,发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为nn。从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n。,2 伪随机数周期和最大容量,乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定：其中a为常数。,2 伪随机数乘同余方法,对于任意正整数M，根据数论中的标准分解定理，总可以分解成如下形式：其中P0=2，P1，Pr表示不同的奇素数，0表示非负整数，1，r表示正整数。a无论取什么值，乘同余方法的最大容量的上界为：的最小公倍数。其中：,2 伪随机数乘同余方法(上限）,关于a与x1的取值,如果a与x1满足如下条件：对于，x1与M互素，则乘同余方法产生的伪随机数序列的最大容量达到最大可能值(M)。,2 伪随机数乘同余方法,乘同余方法在计算机上的使用,为了便于在计算机上使用，通常取：=2s 其中s为计算机中二进制数的最大可能有效位数x1=奇数 a=52k+1 其中k为使52k+1在计算机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，a=513；s=48，a=515等。伪随机数序列的最大容量(M)=2s-2。乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。,2 伪随机数乘同余方法,1948年由Lehmer提出的一种产生伪随机数的方法，是最常用的方法。,1、递推公式：,其中：,I0：初始值（种子seed)a：乘法器（multiplier)c：增值（additive constant)m：模数（modulus)mod：取模运算：(aIn+c)除以m后的余数,a,c和m皆为整数 产生整型的随机数序列,随机性来源于取模运算,如果c=0 乘同余法：速度更快，也可产生长的随机数序列,2 伪随机数乘加同余方法,2、实型随机数序列：,3、特点：,1）最大容量为m：,2 伪随机数乘加同余方法,4、模数m的选择：,m 应尽可能地大，因为序列的周期不可能大于m;,通常将m取为计算机所能表示的最大的整型量，在32位计算机上，m=231=2x109,5、乘数因子a的选择：,1961年，M.Greenberger证明：用线性乘同余方法产生的随机数序列具有周期m的条件是：,c和m为互质数；a-1是质数p的倍数，其中p是a-1和m的共约数；如果m是4的倍数，a-1也是4的倍数。,例：a=5,c=1,m=16,I0=1 周期=m=16 1,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,1,6,15,12,13,2,.,2 伪随机数乘加同余方法,RANDU随机数产生器：,1961年由IBM提出,unsigned long seed=9;float randu()const unsigned long a=65539;const unsigned long m=pow(2,31);unsigned long i1;i1=(a*seed)%m;seed=i1;return(float)i1/float(m);,void SetSeed(unsigned long i)seed=i;,2 伪随机数,存在严重的问题：Marsaglia效用，存在于所有乘同余方法的产生器,void test()c1=new TCanvas(c1,“Test of random number generator,200,10,700,900);pad1=new TPad(pad1,“one,0.03,0.62,0.50,0.92,21);pad2=new TPad(pad2,“one vs one,0.51,0.62,0.98,0.92,21);pad3=new TPad(pad3,“one vs one vs one,0.03,0.02,0.97,0.57,21);pad1-Draw();pad2-Draw();pad3-Draw();TH1F*h1=new TH1F(h1,h1,100,0.0,1.0);TH2F*h2=new TH2F(h2,h2,100,0.0,1.0,100,0.0,1.0);TH3F*h3=new TH3F(h3,h3,100,0.0,1.0,100,0.0,1.0,100,0.0,1.0);,2 伪随机数,for(int i=0;i Fill(x);h2-Fill(x,y);h3-Fill(x,y,z);pad1-cd();h1-Draw();pad2-cd();h2-Draw();pad3-cd();h3-Draw();,2 伪随机数,2 伪随机数,如果取a=69069,将极大地改善结果,2 伪随机数,1968年，Marsaglia对这一问题进行了研究，认为：,任何的乘同余产生器都存在这一问题：在三维和三维以上的空间中，所产生的随机数总是集聚在一些超平面上,随机数序列是关联的,改进措施：将递推公式修改为,特点：1）需要两个初始值（种子）；2）周期可大于m;,2 伪随机数,#include unsigned long seed0=9;unsigned long seed1=11;float randac()const unsigned long a=65539;const unsigned long b=65539;unsigned long i2;unsigned long m=pow(2,31);i2=(a*seed1+b*seed0)%m;seed0=seed1;seed1=i2;return(float)i1/float(m);,void SetSeed(unsigned long i0,unsigned long i1)seed0=i0;seed1=i1;,2 伪随机数,a=b=65539,seed0=9,seed1=11,2 伪随机数,如何获取0,1区间均匀分布的随机数产生器：,每一个Monte Carlo模拟程序软件包都有自带的产生器：,.,利用CERN程序库：,2 伪随机数,CLHEP(Class Library for High Energy Physics)中的随机数产生器,利用CLHEP中的随机数产生器软件包：,http:/,2 伪随机数,FORTRAN中使用随机数产生器应注意的问题：,FORTRAN编译器在对程序进行优化时：,X=RAND(IDUM)+RAND(IDUM)X=2.0*RAND(IDUM),解决办法：,DO I=1,10 IDUM=IDUM+1 X=RAND(IDUM)END DO,2 伪随机数,判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。,2 伪随机数均匀性和独立性,伪随机数的均匀性,这里只考虑伪随机数序列1,2，n全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对于任意的0 x1，令Nn(x)表示伪随机数序列1,2，n中适合不等式i x i=1，2，n 的个数，则 标志伪随机数序列1,2，n的均匀程度，称为均匀偏度。,2 伪随机数均匀性和独立性,将伪随机数序列1,2，n从小至大重新排列 并令，则由(n)的定义，容易证明 很明显，对于固定的，(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立：当 时，所对应的伪随机数序列为最佳分布。,2 伪随机数均匀性和独立性,可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列 的所有值。对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数序列的均匀偏度 对于乘加同余方法 对于部分伪随机数的均匀性问题通常用统计检验方法检验。,2 伪随机数均匀性和独立性,伪随机数的独立性,对于任意，令 表示(1,2)，(2,3)，(n,n+1)中适合不等式 的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令 则(n)标志伪随机数序列1,2，n的独立程度，简称为独立偏度。对于固定的n,(n)的值越接近于零，伪随机数序列的独立性越好。,2 伪随机数均匀性和独立性,对于乘同余方法，对于乘加同余方法，因此，这两种方法的独立性都是很好的。同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出(i)是相当因难的。因此，伪随机数序列的独立性问题的统计检验方法同样是非常重要的。,2 伪随机数均匀性和独立性,40,例子,41,例子,