Presentation is loading. Please wait.

Presentation is loading. Please wait.

2.3 线性乘同余方法 (Linear Congruential Method)

Similar presentations


Presentation on theme: "2.3 线性乘同余方法 (Linear Congruential Method)"— Presentation transcript:

1 2.3 线性乘同余方法 (Linear Congruential Method)
Monte Carlo 模拟 第二章 均匀分布随机数的产生 2.3 线性乘同余方法 (Linear Congruential Method)

2 2.3 线性乘同余方法(Linear Congruential Method)
1948年由Lehmer提出的一种产生伪随机数的方法,是最常用的方法。 1、递推公式: 其中: I0: 初始值(种子seed) a: 乘法器 (multiplier) c: 增值(additive constant) m: 模数(modulus) mod:取模运算:(aIn+c)除以m后的余数 a, c和m皆为整数 产生整型的随机数序列,随机性来源于取模运算 如果c=0  乘同余法:速度更快,也可产生长的随机数序列 mod:取模运算:(aIn+c)除以m后的余数 实型随机数序列:

3 2.3 线性乘同余方法(Linear Congruential Method)
2、实型随机数序列: 3、特点: 1)最大容量为m: 2)独立性和均匀性取决于参数a和c的选择 例:a=c=I0=7, m=10  7,6,9,0,7,6,9,0,…

4 2.3 线性乘同余方法(Linear Congruential Method)
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,..

5 2.3 线性乘同余方法(Linear Congruential Method)
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; }

6 2.3 线性乘同余方法(Linear Congruential Method)
存在严重的问题: 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);

7 2.3 线性乘同余方法(Linear Congruential Method)
for(int i=0; i < 10000; i++) { float x = randu(); float y = randu(); float z = randu(); h1->Fill(x); h2->Fill(x,y); h3->Fill(x,y,z); } pad1->cd(); h1->Draw(); pad2->cd(); h2->Draw(); pad3->cd(); h3->Draw();

8 2.3 线性乘同余方法(Linear Congruential Method)

9 2.3 线性乘同余方法(Linear Congruential Method)

10 2.3 线性乘同余方法(Linear Congruential Method)
1968年,Marsaglia对这一问题进行了研究,认为: 任何的乘同余产生器都存在这一问题:在三维和三维以上的空间中,所产生的随机数总是集聚在一些超平面上 随机数序列是关联的 对于32位的计算机,在d-维空间中超平面的最大数目为 d= d= d= d= 改进措施:将递推公式修改为 特点:1)需要两个初始值(种子); 2)周期可大于m;

11 2.3 线性乘同余方法(Linear Congruential Method)
#include <math.h> 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; }

12 2.3 线性乘同余方法(Linear Congruential Method)
a=b=65539, seed0=9, seed1=11

13 2.3 线性乘同余方法(Linear Congruential Method)
如何获取[0,1]区间均匀分布的随机数产生器: 每一个Monte Carlo模拟程序软件包都有自带的产生器: Jetset(LUND Monte Carlo模拟系列):利用Marsaglia等所提出的算法,周期可达1043 函数用法:r=rlu(idummy) Geant3(探测器模拟程序,FORTRAN): 周期=1018 Call grndm(vec*,len) …. 利用CERN程序库: Y=rndm(x): 周期:5x108 Y=rn32(dummy):乘同余法,a=69069,i0=65539 Call ranmar(vec,len): 周期:1043 Call ranecu(vec,len,isq)

14 2.3 线性乘同余方法(Linear Congruential Method)
利用CLHEP中的随机数产生器软件包: CLHEP(Class Library for High Energy Physics)中的随机数产生器

15 2.3 线性乘同余方法(Linear Congruential Method)
FORTRAN中使用随机数产生器应注意的问题: 在FORTRAN中,如果随机数产生器是带dummy变量的函数: 其中变量idum在函数中不使用,应注意以下问题: X=RAND(idum) FORTRAN编译器在对程序进行优化时: X=RAND(IDUM)+RAND(IDUM)  X=2.0*RAND(IDUM) DO I=1,10 X=RAND(IDUM) END DO …. 解决办法: DO I=1,10 IDUM = IDUM +1 X=RAND(IDUM) END DO


Download ppt "2.3 线性乘同余方法 (Linear Congruential Method)"

Similar presentations


Ads by Google