1、实验数据处理方法 第二部分：Monte Carlo模拟,第七章 均匀分布随机数的产生,随机数的定义和特性,随机数的产生,线性乘同余方法,第七章 均匀分布随机数的产生,7.1 随机数的定义和特性,7.1 随机数的定义和特性,什么是随机数？,单个的数字不是随机数,是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；,在一个均匀分布的随机数中，每一个体出现的概率是均等的；,例如：在0,1区间上均匀分布的随机数序列中，0.00001与0.5出现的机会均等,7.1 随机数的定义和特性,随机数应具有的基本特性,考虑一个对高能粒子反应过程的模拟：需用随机数确定：,出射粒子的属性：能量、方向、,
2、粒子与介质的相互作用,对这一过程的模拟应满足以下要求（相空间产生过程）：,出射粒子的属性应是互不相关的，即每一粒子的属性的确定独立于其它的粒子的属性的确定；,粒子的属性的分布应满足物理所要求的理论分布；,所模拟的物理过程要求随机数应具有下列特性：,随机数序列应是独立的、互不相关的(uncorrelated)：,即序列中的任一子序列应与其它的子序列无关；,7.1 随机数的定义和特性,长的周期(long period)：,实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；,均匀分布的随机数应满足均匀性(Uniformity)：,随机数序列应是均匀的、无偏
3、的，即：如果两个子区间的“面积”相等，则落于这两个子区间内的随机数的个数应相等。,例如：对0,1)区间均匀分布的随机数，如果产生了足够多的随机数，而有一半的随机数落于区间0,0.1不满足均匀性,如果均匀性不满足，则会出现序列中的多组随机数相关的情况均匀性与互不相关的特性是有联系的,7.1 随机数的定义和特性,有效性（Efficiency):,模拟结果可靠,模拟产生的样本容量大,所需的随机数的数量大,随机数的产生必须快速、有效，最好能够进行并行计算。,第七章 均匀分布随机数的产生,7.2 随机数的产生,7.2 随机数的产生,0,1区间上均匀分布的随机数是Monte Carlo模拟的基础:,0,1
4、均匀分布的随机数的产生方法：,利用一些具有内在的随机性的过程：,缺点：模拟的结果不可再现，使得模拟程序的找错困难,利用事先制订好的随机数表:,缺点：表的容量有限，不适合需要大量随机数的应用,服从任意分布的随机数序列可以用0,1区间均匀分布的随机数序列作适当的变换或舍选后求得,7.2 随机数的产生,利用数学递推公式在计算机中产生随机数,其中：T为某个函数，给定初值r1,r2,rk,可按上式确定rn+1, n=1,2, 随机数序列.,算法：产生0，M区间上的整数In，然后利用公式rn=In/M返回0,1区间上的实数,优点：,占用计算机的内存少； 产生速度快； 可以重复前次的模拟结果，便于程序的找错
5、；,7.2 随机数的产生,缺点：,不满足随机数之间相互独立的要求：公式和初值确定后，序列就唯一地确定了；,伪随机数（Pseudo-Random Number),不满足均匀性：计算机能表示的0,1区间内的数是有限的（由字长确定）递推到一定次数后，出现周期性的重复现象,第七章 均匀分布随机数的产生,7.3 线性乘同余方法 （Linear Congruential Method),7.3 线性乘同余方法 （Linear Congruential Method),mod:取模运算：(aIn+c)除以m后的余数实型随机数序列：,1948年由Lehmer提出的一种产生伪随机数的方法，是最常用的方法。,1、
6、递推公式：,其中：,I0： 初始值（种子seed) a： 乘法器 （multiplier) c： 增值（additive constant) m： 模数（modulus) mod：取模运算：(aIn+c)除以m后的余数,a, c和m皆为整数产生整型的随机数序列,随机性来源于取模运算,如果c=0 乘同余法：速度更快，也可产生长的随机数序列,7.3 线性乘同余方法 （Linear Congruential Method),2、实型随机数序列：,3、特点：,1）最大容量为m：,7.3 线性乘同余方法 （Linear Congruential Method),4、模数m的选择：,m 应尽可能地大，因为
7、序列的周期不可能大于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=161,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,1,6,15, 12,13,2,7.3 线性乘同余方法 （Linear Congruential Method),RANDU随
8、机数产生器：,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; ,7.3 线性乘同余方法 （Linear Congruential Method),存在严重的问题：Marsaglia效用，存
9、在于所有乘同余方法的产生器,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-Dra
10、w(); 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);,7.3 线性乘同余方法 （Linear Congruential Method),for(int i=0; i Fill(x);h2-Fill(x,y);h3-Fill(x,y,z); pad1-cd(); h1-Draw();
11、pad2-cd(); h2-Draw();pad3-cd(); h3-Draw(); ,7.3 线性乘同余方法 （Linear Congruential Method),7.3 线性乘同余方法 （Linear Congruential Method),如果取a=69069,将极大地改善结果,7.3 线性乘同余方法 （Linear Congruential Method),1968年，Marsaglia对这一问题进行了研究，认为：,任何的乘同余产生器都存在这一问题：在三维和三维以上的空间中，所产生的随机数总是集聚在一些超平面上,随机数序列是关联的,改进措施：将递推公式修改为,特点：1）需要两个初
12、始值（种子）；2）周期可大于m;,7.3 线性乘同余方法 （Linear Congruential Method),#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;r
13、eturn (float) i1/float(m); ,void SetSeed(unsigned long i0,unsigned long i1) seed0 = i0;seed1 = i1; ,7.3 线性乘同余方法 （Linear Congruential Method),a=b=65539, seed0=9, seed1=11,7.3 线性乘同余方法 （Linear Congruential Method),如何获取0,1区间均匀分布的随机数产生器：,每一个Monte Carlo模拟程序软件包都有自带的产生器：,.,利用CERN程序库：,7.3 线性乘同余方法 （Linear Congruential Method),CLHEP(Class Library for High Energy Physics)中的随机数产生器,利用CLHEP中的随机数产生器软件包：,http:/ 线性乘同余方法 （Linear Congruential Method),FORTRAN中使用随机数产生器应注意的问题：,FORTRAN编译器在对程序进行优化时：,X=RAND(IDUM)+RAND(IDUM) X=2.0*RAND(IDUM),解决办法：,DO I=1,10IDUM = IDUM +1X=RAND(IDUM) END DO,
