Monte carlo
This presentation is the property of its rightful owner.
Sponsored Links
1 / 15

Monte Carlo 模拟 PowerPoint PPT Presentation


  • 67 Views
  • Uploaded on
  • Presentation posted in: General

Monte Carlo 模拟. 第二章 均匀分布随机数的产生. 2.3 线性乘同余方法 ( Linear Congruential Method). 1948 年由 Lehmer 提出的一种产生伪随机数的方法,是最常用的方法。. 2.3 线性乘同余方法 ( Linear Congruential Method). 1 、递推公式:. mod: 取模运算: ( aI n +c ) 除以 m 后的余数 实型随机数序列:. 其中:. I 0 : 初始值(种子 seed) a : 乘法器 ( multiplier)

Download Presentation

Monte Carlo 模拟

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Monte carlo

Monte Carlo 模拟

第二章均匀分布随机数的产生

2.3 线性乘同余方法

(Linear Congruential Method)


2 3 linear congruential method

1948年由Lehmer提出的一种产生伪随机数的方法,是最常用的方法。

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

1、递推公式:

mod:取模运算:(aIn+c)除以m后的余数

实型随机数序列:

其中:

I0: 初始值(种子seed)

a: 乘法器 (multiplier)

c: 增值(additive constant)

m: 模数(modulus)

mod:取模运算:(aIn+c)除以m后的余数

a, c和m皆为整数

产生整型的随机数序列,随机性来源于取模运算

如果c=0  乘同余法:速度更快,也可产生长的随机数序列


2 3 linear congruential method1

1)最大容量为m:

2)独立性和均匀性取决于参数a和c的选择

例:a=c=I0=7, m=10  7,6,9,0,7,6,9,0,…

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

2、实型随机数序列:

3、特点:


2 3 linear congruential method2

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

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 3 linear congruential method3

1961年由IBM提出

RANDU随机数产生器:

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

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 3 linear congruential method4

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);


2 3 linear congruential method5

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();

}


2 3 linear congruential method6

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


2 3 linear congruential method7

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

如果取a=69069,将极大地改善结果


2 3 linear congruential method8

对于32位的计算机,在d-维空间中超平面的最大数目为

d=3 2953

d=4 566

d=6 120

d=10 41

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

1968年,Marsaglia对这一问题进行了研究,认为:

任何的乘同余产生器都存在这一问题:在三维和三维以上的空间中,所产生的随机数总是集聚在一些超平面上

随机数序列是关联的

改进措施:将递推公式修改为

特点:1)需要两个初始值(种子);

2)周期可大于m;


Monte carlo

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;

}


2 3 linear congruential method9

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

a=b=65539, seed0=9, seed1=11


2 3 linear congruential method10

  • Jetset(LUND Monte Carlo模拟系列):利用Marsaglia等所提出的算法,周期可达1043

函数用法:r=rlu(idummy)

  • Geant3(探测器模拟程序,FORTRAN): 周期=1018

Call grndm(vec*,len)

  • Y=rndm(x): 周期:5x108

  • Y=rn32(dummy):乘同余法,a=69069,i0=65539

  • Call ranmar(vec,len): 周期:1043

  • Call ranecu(vec,len,isq)

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

如何获取[0,1]区间均匀分布的随机数产生器:

  • 每一个Monte Carlo模拟程序软件包都有自带的产生器:

….

  • 利用CERN程序库:


2 3 linear congruential method11

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

  • 利用CLHEP中的随机数产生器软件包:

CLHEP(Class Library for High Energy Physics)中的随机数产生器

http://hepg.sdu.edu.cn/~zhangxy/clhep/html/index.html


2 3 linear congruential method12

在FORTRAN中,如果随机数产生器是带dummy变量的函数:

X=RAND(idum)

其中变量idum在函数中不使用,应注意以下问题:

  • DO I=1,10

  • X=RAND(IDUM)

  • END DO

X=RAND(IDUM)

DO I=1,10

….

END DO

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

FORTRAN中使用随机数产生器应注意的问题:

FORTRAN编译器在对程序进行优化时:

X=RAND(IDUM)+RAND(IDUM)  X=2.0*RAND(IDUM)

解决办法:

  • DO I=1,10

  • IDUM = IDUM +1

  • X=RAND(IDUM)

  • END DO


  • Login