蒙特卡罗方法 (Monte Carlo simulation)

2014-11-12
均匀分布随机数的产生
15
2.3 线性乘同余方法
2014-11-12
均匀分布随机数的产生
16
2.3 线性乘同余方法
如果取a=69069,将极大地改善结果
2014-11-12
均匀分布随机数的产生
17
2.3 线性乘同余方法
1968年，Marsaglia对这一问题进行了研究，认为： 任何的乘同余产生器都存在这一问题：在三维和三维以 上的空间中，所产生的随机数总是集聚在一些超平面上 随机数序列是关联的 对于32位的计算机，在d-维空间中超平面的最大数目为
2014-11-12 均匀分布随机数的产生 12
2.3 线性乘同余方法
RANDU随机数产生器：
1961年由IBM提出
I n1 (65539 I n ) mod231
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);
算法：产生[0，M]区间上的整数In，然后利用公式 rn=In/M返回[0,1]区间上的实数
优点： – 占用计算机的内存少； – 产生速度快； – 可以重复前次的模拟结果，便于程序的找错；
2014-11-12 均匀分布随机数的产生 7
2.2 随机数的产生
缺点：
• 不满足随机数之间相互独立的要求：公式和初值确定 后，序列就唯一地确定了； • 不满足均匀性：计算机能表示的[0,1]区间内的数是有限 的（由字长确定）
2014-11-12
均匀分布随机数的产生
2
2.1 随机数的定义和特性
随机数应具有的基本特性
考虑一个对高能粒子反应过程的模拟：需用随机数确定：
出射粒子的属性：能量、方向、… 粒子与介质的相互作用 对这一过程的模拟应满足以下要求（相空间产生过程）：
出射粒子的属性应是互不相关的，即每一粒 子的属性的确定独立于其它的粒子的属性的 确定； 粒子的属性的分布应满足物理所要求的理论 分布；
void SetSeed(unsigned long i0, unsigned long i1) { seed0 = i0; seed1 = i1; }
unsigned long i2;
unsigned long m = pow(2,31); i2 = (a * seed1 + b * seed0 ) % m; seed0 = seed1; seed1 = i2; return (float) i1/float(m); }
9
2.3 线性乘同余方法
1948年由Lehmer提出的一种产生伪随机数的方法，是最常用的方 法。 1、递推公式：
I n1 (aIn c) modm
其中： I0： 初始值（种子seed) a： 乘子 （multiplier) c： 常数（additive constant) m： 模数（modulus) mod： 取模运算：(aIn+c)除以m后的余数
d=3 2953 d=4 566 d=6 120 d=10 41
改进措施：将递推公式修改为
I n (a I n1 b I n2 ) modm
特点：1）需要两个初始值（种子）； 2）周期可大于m;
2014-11-12 均匀分布随机数的产生 18
2.3 线性乘同余方法
#include <math.h> unsigned long seed0 = 9; unsigned long seed1 = 11; float randac() { const unsigned long a = 65539; const unsigned long b = 65539;
Monte Carlo模拟
第二章 均匀分布随机数的产生
2.1 随机数的定义和特性 2.2 随机数的产生 2.3 线形乘同余方法
2014-11-12
均匀分布随机数的产生
1
2.1 随机数的定义和特性
什么是随机数？ 单个的数字不是随机数 是指一个数列，其中的每一个体称为随机数，其值与数列中 的其它数无关； 在一个均匀分布的随机数中，每一个体出现的概率是均等的； 例如：在[0,1]区间上均匀分布的随机数序列中， 0.00001与0.5出现的机会均等
}
2014-11-12 均匀分布随机数的产生 13
void SetSeed(unsigned long i) { seed = i; }
2.3 线性乘同余方法
存在严重的问题： 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. 利用CERN程序库：
• Y=rndm(x): 周期：5x108 • Y=rn32(dummy):乘同余法，a=69069,i0=65539 • Call ranmar(vec,len): 周期：1043 • Call ranecu(vec,len,isq)
2014-11-12 均匀分布随机数的产生 21
利用事先制订好的随机数表:
缺点：表的容量有限，不适合需要大量随机数的应用
2014-11-12 均匀分布随机数的产生 6
2.2 随机数的产生
利用数学递推公式在计算机中产生随机数
rn k T (rn , rn1 ,..., rn k 1 )
其中：T为某个函数，给定初值r1,r2,…,rk,可按上式确 定rn+1, n=1,2,… 随机数序列.
其中变量idum在函数中不使用，应注意以下问题： FORTRAN编译器在对程序进行优化时：
X=RAND(IDUM)+RAND(IDUM) X=2.0*RAND(IDUM)
2014-11-12
均匀分布随机数的产生
19
2.3 线性乘同余方法
a=b=65539, seed0=9, seed1=11
2014-11-12
均匀分布随机数的产生
20
2.3 线性乘同余方法
如何获取[0,1]区间均匀分布的随机数产生器：
1. 每一个Monte Carlo模拟程序软件包都有自带的产生器： • Jetset(LUND Monte Carlo模拟系列）：利用Marsaglia等 所提出的算法，周期可达1043 函数用法：r=rlu(idummy) • Geant3(探测器模拟程序，FORTRAN): 周期=1018 Call grndm(vec*,len)
3. 均匀分布的随机数应满足均匀性(Uniformity)：
随机数序列应是均匀的、无偏的，即：如果两个子区间的“面积”相等， 则落于这两个子区间内的随机数的个数应相等。
4. 有效性（Efficiency):
模拟结果可靠模拟产生的样本容量大所需的随机数的数量大 随机数的产生必须快速、有效，最好能够进行并 行计算。
2014-11-12
均匀分布随机数的产生
3
2.1 随机数的定义和特性
所模拟的物理过程要求随机数应具有下列特性：
1. 随机数序列应是独立的、互不相关的(uncorrelated)：
即序列中的任一子序列应与其它的子序列无关；
2. 长的周期(long period)：
实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期 性，即当序列达到一定长度后会重复；
递推到一定次数后，出现周期性的重复现象
伪随机数（Pseudo-Random Number)
2014-11-12
均匀分布随机数的产生
8
Monte Carlo模拟
第二章 均匀分布随机数的产生
2.1 随机数的定义和特性 2.2 随机数的产生 2.3 线形乘同余方法
2014-11-12
均匀分布随机数的产生
2014-11-12 均匀分布随机数的产生 14
2.3 线性乘同余方法
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(); }
In m I n m 1
3、特点： 1）最大容量为m：0 I n m
2）独立性和均匀性取决于参数a和c的选择 例：a=c=I0=7, m=10 7,6,9,0,7,6,9,0,…
2014-11-12
均匀分布随机数的产生
11
2.3 线性乘同余方法
4、模数m的选择： • m 应尽可能地大，因为序列的周期不可能大于m; • 通常将m取为计算机所能表示的最大的整型量，在32位计算 机上，m=231=2x109 5、乘数因子a的选择： 1961年，M. Greenberger证明：用线性乘同余方法产生的随机 数序列具有周期m的条件是： 1. c和m为互质数； 2. a-1是质数p的倍数，其中p是a-1和m的共约数； 3. 如果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 线性乘同余方法
3. 利用CLHEP中的随机数产生器软件包： CLHEP(Class Library for High Energy Physics)中的随机数产 生器
2014-11-12
均匀分布随机数的产生
22
2.3 线性乘同余方法
FORTRAN中使用随机数产生器应注意的问题： 在FORTRAN中，如果随机数产生器是带dummy变量的函数： X=RAND(idum)
a, c 0 m I 0 , a, c
a, c和m皆为整数 产生整型的随机数序列,随机性来源于取模运算
如果c=0 乘同余法：速度更快，也可产生长的随机数序列
2014-11-12 均匀分布随机数的产生 10
2.3 线性乘同余方法
2、实型随机数序列：
In rn [0,1) float(m) In rn [0,1] float(m 1)
2014-11-12 均匀分布随机数的产生 4
Monte Carlo模拟
第二章 均匀分布随机数的产生
2.1 随机数的定义和特性 2.2 随机数的产生 2.3 线形乘同余方法
2014-11-12
均匀分布随机数的产生
5
Leabharlann Baidu
2.2 随机数的产生
• [0,1]区间上均匀分布的随机数是Monte Carlo模拟的基础: 服从任意分布的随机数序列可以用[0,1]区间均匀分布的 随机数序列作适当的变换或舍选后求得 • [0,1]均匀分布的随机数的产生方法： 利用一些具有内在的随机性的过程： 放射性衰变过程（radioactive decay)； 热噪声(thermal noise); 宇宙线的到达时间（cosmic ray arrival); … 缺点：模拟的结果不可再现，使得模拟程序的找错困难