MonteCarlo模擬

1、2.3 線性乘同余方法線性乘同余方法（Linear Congruential Method)mcaIInnmod)(1caImca,0,01948年由年由Lehmer提出的一種產生偽隨機數(shù)的方法，是最常用的方提出的一種產生偽隨機數(shù)的方法，是最常用的方法。法。1 1、遞推公式：、遞推公式：其中：其中：I0： 初始值（種子初始值（種子seed) a： 乘法器乘法器 （multiplier) c： 增值（增值（additive constant) m： 模數(shù)（模數(shù)（modulus) mod：取模運算：取模運算：(aIn+c)除以除以m后的余數(shù)后的余數(shù)a, c和和m皆為整數(shù)皆為整數(shù) 產生整型的隨機數(shù)序

2、列產生整型的隨機數(shù)序列,隨機性來源于取模運算隨機性來源于取模運算如果如果c=0 乘同余法：速度更快，也可產生長的隨機數(shù)序列乘同余法：速度更快，也可產生長的隨機數(shù)序列1 ,0)1()1 ,0)(mfloatIrmfloatIrnnnn1mImInn2 2、實型隨機數(shù)序列：、實型隨機數(shù)序列：3、特點：、特點：1）最大容量為）最大容量為m：mIn02）獨立性和均勻性取決于參數(shù)）獨立性和均勻性取決于參數(shù)a和和c的選擇的選擇例：例：a=c=I0=7, m=10 7,6,9,0,7,6,9,0,4、模數(shù)、模數(shù)m的選擇：的選擇： m 應盡可能地大，因為序列的周期不可能大于應盡可能地大，因為序列的周期不可能大

3、于m; 通常將通常將m取為計算機所能表示的最大的整型量，在取為計算機所能表示的最大的整型量，在32位計算位計算機上，機上，m=231=2x1095、乘數(shù)因子、乘數(shù)因子a的選擇：的選擇：1961年，年，M. Greenberger證明：用線性乘同余方法產生的隨機證明：用線性乘同余方法產生的隨機數(shù)序列具有周期數(shù)序列具有周期m的條件是：的條件是：1. c和和m為互質數(shù)；為互質數(shù)；2. a-1是質數(shù)是質數(shù)p的倍數(shù)，其中的倍數(shù)，其中p是是a-1和和m的共約數(shù)；的共約數(shù)；3. 如果如果m是是4的倍數(shù)，的倍數(shù)，a-1也是也是4的倍數(shù)。的倍數(shù)。例：例：a=5,c=1,m=16,I0=1 周期周期=m=16 1

4、,6,15,12,13,2,11,8,9,14,7,4,5,10,3,0,1,6,15, 12,13,2,.RANDU隨機數(shù)產生器：隨機數(shù)產生器：3112mod)65539(nnII1961年由年由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(u

5、nsigned long i) seed = i;存在嚴重的問題：存在嚴重的問題： MarsagliaMarsaglia效用，存在于所有乘同余方法的產生器效用，存在于所有乘同余方法的產生器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 TPa

6、d(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); for(int i=0; i Fill(x); h2-Fill(x,y); h3-Fill(x,y,z

7、); pad1-cd(); h1-Draw(); pad2-cd(); h2-Draw(); pad3-cd(); h3-Draw(); 如果取如果取a=69069,a=69069,將極大地改善結果將極大地改善結果mIbIaInnnmod)(211968年，年，Marsaglia對這一問題進行了研究，認為：對這一問題進行了研究，認為：任何的乘同余產生器都存在這一問題：在三維和三維以任何的乘同余產生器都存在這一問題：在三維和三維以上的空間中，所產生的隨機數(shù)總是集聚在一些超平面上上的空間中，所產生的隨機數(shù)總是集聚在一些超平面上隨機數(shù)序列是關聯(lián)的隨機數(shù)序列是關聯(lián)的對于對于32位的計算機，在位的計算機

8、，在d-維空間中超平面的最大數(shù)目為維空間中超平面的最大數(shù)目為d=3 2953d=4 566d=6 120d=10 41改進措施：將遞推公式修改為改進措施：將遞推公式修改為特點：特點：1）需要兩個初始值（種子）；）需要兩個初始值（種子）； 2）周期可大于）周期可大于m;#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

9、 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;a=b=65539, seed0=9, seed1=11如何獲取如何獲取0,1區(qū)間均勻分布的隨機數(shù)產生器：區(qū)間均勻分布的隨機數(shù)產生器：1. 每一個每一個Monte Carlo模擬程序軟件包都有自帶的產生器：模擬程序軟件包都有自帶的產生器： Jet

10、set(LUND Monte Carlo模擬系列）：利用模擬系列）：利用Marsaglia等等所提出的算法，周期可達所提出的算法，周期可達1043函數(shù)用法：函數(shù)用法：r=rlu(idummy) Geant3(探測器模擬程序，探測器模擬程序，F(xiàn)ORTRAN): 周期周期=1018Call grndm(vec*,len).2. 利用利用CERN程序庫：程序庫： Y=rndm(x): 周期：周期：5x108 Y=rn32(dummy):乘同余法，乘同余法，a=69069,i0=65539 Call ranmar(vec,len): 周期：周期：1043 Call ranecu(vec,len,isq)CLHEP(Class Library for High Energy Physics)中的隨機數(shù)產中的隨機數(shù)產生器生器3. 利用利用CLH