蒙特卡羅方法課件(清華大學(xué)-林謙)

1、蒙特卡羅方法在核技術(shù)中的應(yīng)用,林謙,目 錄,第一章 蒙特卡羅方法概述 第二章 隨機(jī)數(shù) 第三章 由已知分布的隨機(jī)抽樣 第四章 蒙特卡羅方法解粒子輸運(yùn)問題,教材,蒙特卡羅方法在實(shí)驗核物理中的應(yīng)用 許淑艷 編著原子能出版社 蒙特卡羅方法 清華大學(xué),參考書,蒙特卡羅方法及其在粒子輸運(yùn)問題中的應(yīng)用 裴鹿成 張孝澤 編著科學(xué)出版社 蒙特卡羅方法 徐鐘濟(jì) 編著上海科學(xué)技術(shù)出版社,聯(lián)系方式,電話83918 電子郵件,第一章 蒙特卡羅方法概述,蒙特卡羅方法的基本思想 蒙特卡羅方法的收斂性，誤差 蒙特卡羅方法的特點(diǎn) 蒙特卡羅方法的主要應(yīng)用范圍 作 業(yè),第一章 蒙特卡羅方法概述,蒙特卡羅方法又稱隨機(jī)抽樣技巧或統(tǒng)計試

2、驗方法。半個多世紀(jì)以來，由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明 ，這種方法作為一種獨(dú)立的方法被提出來，并首先在核武器的試驗與研制中得到了應(yīng)用。蒙特卡羅方法是一種計算方法，但與一般數(shù)值計算方法有很大區(qū)別。它是以概率統(tǒng)計理論為基礎(chǔ)的一種方法。由于蒙特卡羅方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗過程，解決一些數(shù)值方法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。,蒙特卡羅方法的基本思想,二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨(dú)立的方法被提出來，并首先在核武器的試驗與研制中得到了應(yīng)用。但其基本思想并非新穎，人們在生產(chǎn)實(shí)踐和科學(xué)試驗中就已發(fā)現(xiàn)，并加以利用。 兩個

3、例子 例1. 蒲豐氏問題 例2. 射擊問題（打靶游戲） 基本思想 計算機(jī)模擬試驗過程,例1. 蒲豐氏問題,為了求得圓周率值，在十九世紀(jì)后期，有很多人作了這樣的試驗：將長為2l的一根針任意投到地面上，用針與一組相間距離為2a（ la）的平行線相交的頻率代替概率P，再利用準(zhǔn)確的關(guān)系式： 求出值 其中為投計次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。,一些人進(jìn)行了實(shí)驗，其結(jié)果列于下表 ：,例2. 射擊問題（打靶游戲）,設(shè)r表示射擊運(yùn)動員的彈著點(diǎn)到靶心的距離，(r)表示擊中r處相應(yīng)的得分?jǐn)?shù)（環(huán)數(shù)），f(r)為該運(yùn)動員的彈著點(diǎn)的分布密度函數(shù)，它反映運(yùn)動員的射擊水平。該運(yùn)動員的射擊成

4、績?yōu)?用概率語言來說，是隨機(jī)變量(r)的數(shù)學(xué)期望，即,現(xiàn)假設(shè)該運(yùn)動員進(jìn)行了次射擊，每次射擊的彈著點(diǎn)依次為r1，r2，rN，則次得分g(r1)，g(r2)，g(rN)的算術(shù)平均值 代表了該運(yùn)動員的成績。換言之，為積分的估計值，或近似值。 在該例中，用次試驗所得成績的算術(shù)平均值作為數(shù)學(xué)期望的估計值（積分近似值）。,基本思想,由以上兩個例子可以看出，當(dāng)所求問題的解是某個事件的概率，或者是某個隨機(jī)變量的數(shù)學(xué)期望，或者是與概率、數(shù)學(xué)期望有關(guān)的量時，通過某種試驗的方法，得出該事件發(fā)生的頻率，或者該隨機(jī)變量若干個具體觀察值的算術(shù)平均值，通過它得到問題的解。這就是蒙特卡羅方法的基本思想。 當(dāng)隨機(jī)變量的取值僅為

5、1或0時，它的數(shù)學(xué)期望就是某個事件的概率?；蛘哒f，某種事件的概率也是隨機(jī)變量（僅取值為1或0）的數(shù)學(xué)期望。,因此，可以通俗地說，蒙特卡羅方法是用隨機(jī)試驗的方法計算積分，即將所要計算的積分看作服從某種分布密度函數(shù)f(r)的隨機(jī)變量(r)的數(shù)學(xué)期望 通過某種試驗，得到個觀察值r1，r2，rN（用概率語言來說，從分布密度函數(shù)f(r)中抽取個子樣r1，r2，rN，），將相應(yīng)的個隨機(jī)變量的值g(r1)，g(r2)，g(rN)的算術(shù)平均值 作為積分的估計值（近似值）。,為了得到具有一定精確度的近似解，所需試驗的次數(shù)是很多的，通過人工方法作大量的試驗相當(dāng)困難，甚至是不可能的。因此，蒙特卡羅方法的基本思想雖然

6、早已被人們提出，卻很少被使用。本世紀(jì)四十年代以來，由于電子計算機(jī)的出現(xiàn)，使得人們可以通過電子計算機(jī)來模擬隨機(jī)試驗過程，把巨大數(shù)目的隨機(jī)試驗交由計算機(jī)完成，使得蒙特卡羅方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中發(fā)揮應(yīng)有的作用。,計算機(jī)模擬試驗過程,計算機(jī)模擬試驗過程，就是將試驗過程（如投針，射擊）化為數(shù)學(xué)問題，在計算機(jī)上實(shí)現(xiàn)。以上述兩個問題為例，分別加以說明。 例1. 蒲豐氏問題 例2. 射擊問題（打靶游戲） 由上面兩個例題看出，蒙特卡羅方法常以一個“概率模型”為基礎(chǔ)，按照它所描述的過程，使用由已知分布抽樣的方法，得到部分試驗結(jié)果的觀察值，求得問題的近似解。,例蒲豐氏問題,設(shè)針投到地面上的位置可以

7、用一組參數(shù)（x,）來描述，x為針中心的坐標(biāo)，為針與平行線的夾角，如圖所示。 任意投針，就是意味著x與都是任意取的，但x的范圍限于0，a，夾角的范圍限于0，。在此情況下，針與平行線相交的數(shù)學(xué)條件是,針在平行線間的位置,如何產(chǎn)生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均勻分布的，其分布密度函數(shù)為： 類似地，的分布密度函數(shù)為： 因此，產(chǎn)生任意的（x,）的過程就變成了由f1(x)抽樣x及由f2()抽樣的過程了。由此得到： 其中1，2均為（0,1）上均勻分布的隨機(jī)變量。,每次投針試驗，實(shí)際上變成在計算機(jī)上從兩個均勻分布的隨機(jī)變量中抽樣得到（x,），然后定義描述針與平行線相交狀況的隨機(jī)變量

8、s(x,)，為 如果投針次，則 是針與平行線相交概率的估計值。事實(shí)上， 于是有,例射擊問題,設(shè)射擊運(yùn)動員的彈著點(diǎn)分布為 用計算機(jī)作隨機(jī)試驗（射擊）的方法為，選取一個隨機(jī)數(shù)，按右邊所列方法判斷得到成績。 這樣，就進(jìn)行了一次隨機(jī)試驗（射擊），得到了一次成績 (r)，作次試驗后，得到該運(yùn)動員射擊成績的近似值,蒙特卡羅方法的收斂性，誤差,蒙特卡羅方法作為一種計算方法，其收斂性與誤差是普遍關(guān)心的一個重要問題。 收斂性 誤差 減小方差的各種技巧 效率,收斂性,由前面介紹可知，蒙特卡羅方法是由隨機(jī)變量X的簡單子樣X1，X2，XN的算術(shù)平均值： 作為所求解的近似值。由大數(shù)定律可知， 如X1，X2，XN獨(dú)立同分

9、布，且具有有限期望值（E(X)），則 即隨機(jī)變量X的簡單子樣的算術(shù)平均值 ，當(dāng)子樣數(shù)充分大時，以概率1收斂于它的期望值E(X)。,誤差,蒙特卡羅方法的近似值與真值的誤差問題，概率論的中心極限定理給出了答案。該定理指出，如果隨機(jī)變量序列X1，X2，XN獨(dú)立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函數(shù)。則,當(dāng)N充分大時，有如下的近似式 其中稱為置信度，1稱為置信水平。 這表明，不等式 近似地以概率 1成立，且誤差收斂速度的階為 。 通常，蒙特卡羅方法的誤差定義為 上式中 與置信度是一一對應(yīng)的，根據(jù)問題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出 。,下面給出幾個常用的

10、與的數(shù)值： 關(guān)于蒙特卡羅方法的誤差需說明兩點(diǎn)：第一，蒙特卡羅方法的誤差為概率誤差，這與其他數(shù)值計算方法是有區(qū)別的。第二，誤差中的均方差是未知的，必須使用其估計值 來代替，在計算所求量的同時，可計算出 。,減小方差的各種技巧,顯然，當(dāng)給定置信度后，誤差由和N決定。要減小，或者是增大N，或者是減小方差2。在固定的情況下，要把精度提高一個數(shù)量級，試驗次數(shù)N需增加兩個數(shù)量級。因此，單純增大N不是一個有效的辦法。 另一方面，如能減小估計的均方差，比如降低一半，那誤差就減小一半，這相當(dāng)于N增大四倍的效果。因此降低方差的各種技巧，引起了人們的普遍注意。后面課程將會介紹一些降低方差的技巧。,效率,一般來說，降

11、低方差的技巧，往往會使觀察一個子樣的時間增加。在固定時間內(nèi)，使觀察的樣本數(shù)減少。所以，一種方法的優(yōu)劣，需要由方差和觀察一個子樣的費(fèi)用（使用計算機(jī)的時間）兩者來衡量。這就 是蒙特卡羅方法中效率的概念。它定義為 ，其中c 是觀察一個子樣的平均費(fèi)用。顯然 越小，方法越有效。,蒙特卡羅方法的特點(diǎn),優(yōu)點(diǎn) 能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗過程。 受幾何條件限制小。 收斂速度與問題的維數(shù)無關(guān)。 具有同時計算多個方案與多個未知量的能力。 誤差容易確定。 程序結(jié)構(gòu)簡單，易于實(shí)現(xiàn)。,缺點(diǎn) 收斂速度慢。 誤差具有概率性。 在粒子輸運(yùn)問題中，計算結(jié)果與系統(tǒng)大小有關(guān)。,能夠比較逼真地描述具有隨機(jī)性質(zhì)

12、的事物的特點(diǎn)及物理實(shí)驗過程,從這個意義上講，蒙特卡羅方法可以部分代替物理實(shí)驗，甚至可以得到物理實(shí)驗難以得到的結(jié)果。用蒙特卡羅方法解決實(shí)際問題，可以直接從實(shí)際問題本身出發(fā)，而不從方程或數(shù)學(xué)表達(dá)式出發(fā)。它有直觀、形象的特點(diǎn)。,受幾何條件限制小,在計算s維空間中的任一區(qū)域Ds上的積分 時，無論區(qū)域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產(chǎn)生N個點(diǎn) ，得到積分的近似值。 其中Ds為區(qū)域Ds的體積。這是數(shù)值方法難以作到的。 另外，在具有隨機(jī)性質(zhì)的問題中，如考慮的系統(tǒng)形狀很復(fù)雜，難以用一般數(shù)值方法求解，而使用蒙特卡羅方法，不會有原則上的困難。,收斂速度與問題的維數(shù)無關(guān),由

13、誤差定義可知，在給定置信水平情況下，蒙特卡羅方法的收斂速度為，與問題本身的維數(shù)無關(guān)。維數(shù)的變化，只引起抽樣時間及估計量計算時間的變化，不影響誤差。也就是說，使用蒙特卡羅方法時，抽取的子樣總數(shù)N與維數(shù)s無關(guān)。維數(shù)的增加，除了增加相應(yīng)的計算量外，不影響問題的誤差。這一特點(diǎn)，決定了蒙特卡羅方法對多維問題的適應(yīng)性。而一般數(shù)值方法，比如計算定積分時，計算時間隨維數(shù)的冪次方而增加，而且，由于分點(diǎn)數(shù)與維數(shù)的冪次方成正比，需占用相當(dāng)數(shù)量的計算機(jī)內(nèi)存，這些都是一般數(shù)值方法計算高維積分時難以克服的問題。,具有同時計算多個方案與多個未知量的能力,對于那些需要計算多個方案的問題，使用蒙特卡羅方法有時不需要像常規(guī)方法那

14、樣逐個計算，而可以同時計算所有的方案，其全部計算量幾乎與計算一個方案的計算量相當(dāng)。例如，對于屏蔽層為均勻介質(zhì)的平板幾何，要計算若干種厚度的穿透概率時，只需計算最厚的一種情況，其他厚度的穿透概率在計算最厚一種情況時稍加處理便可同時得到。 另外，使用蒙特卡羅方法還可以同時得到若干個所求量。例如，在模擬粒子過程中，可以同時得到不同區(qū)域的通量、能譜、角分布等，而不像常規(guī)方法那樣，需要逐一計算所求量。,誤差容易確定,對于一般計算方法，要給出計算結(jié)果與真值的誤差并不是一件容易的事情，而蒙特卡羅方法則不然。根據(jù)蒙特卡羅方法的誤差公式，可以在計算所求量的同時計算出誤差。對干很復(fù)雜的蒙特卡羅方法計算問題，也是容

15、易確定的。 一般計算方法常存在著有效位數(shù)損失問題，而要解決這一問題有時相當(dāng)困難，蒙特卡羅方法則不存在這一問題。,程序結(jié)構(gòu)簡單，易于實(shí)現(xiàn),在計算機(jī)上進(jìn)行蒙特卡羅方法計算時，程序結(jié)構(gòu)簡單，分塊性強(qiáng)，易于實(shí)現(xiàn)。,收斂速度慢,如前所述，蒙特卡羅方法的收斂速度為 ，一般不容易得到精確度較高的近似結(jié)果。對于維數(shù)少（三維以下）的問題，不如其他方法好。,誤差具有概率性,由于蒙特卡羅方法的誤差是在一定置信水平下估計的，所以它的誤差具有概率性，而不是一般意義下的誤差。,在粒子輸運(yùn)問題中，計算結(jié)果與系統(tǒng)大小有關(guān),經(jīng)驗表明，只有當(dāng)系統(tǒng)的大小與粒子的平均自由程可以相比較時（一般在十個平均自由程左右），蒙特卡羅方法計算的

16、結(jié)果較為滿意。但對于大系統(tǒng)或小概率事件的計算問題，計算結(jié)果往往比真值偏低。而對于大系統(tǒng)，數(shù)值方法則是適用的。 因此，在使用蒙特卡羅方法時，可以考慮把蒙特卡羅方法與解析（或數(shù)值）方法相結(jié)合，取長補(bǔ)短，既能解決解析（或數(shù)值）方法難以解決的問題，也可以解決單純使用蒙特卡羅方法難以解決的問題。這樣，可以發(fā)揮蒙特卡羅方法的特長，使其應(yīng)用范圍更加廣泛。,蒙特卡羅方法的主要應(yīng)用范圍,蒙特卡羅方法所特有的優(yōu)點(diǎn)，使得它的應(yīng)用范圍越來越廣。它的主要應(yīng)用范圍包括：粒子輸運(yùn)問題，統(tǒng)計物理，典型數(shù)學(xué)問題，真空技術(shù)，激光技術(shù)以及醫(yī)學(xué)，生物，探礦等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。 蒙特卡羅方法在粒子輸運(yùn)問題

17、中的應(yīng)用范圍主要包括：實(shí)驗核物理，反應(yīng)堆物理，高能物理等方面。 蒙特卡羅方法在實(shí)驗核物理中的應(yīng)用范圍主要包括：通量及反應(yīng)率，中子探測效率，光子探測效率，光子能量沉積譜及響應(yīng)函數(shù)，氣體正比計數(shù)管反沖質(zhì)子譜，多次散射與通量衰減修正等方面。,作 業(yè),用蒲豐投針法在計算機(jī)上計算值，取a=4、l=3。 分別用理論計算和計算機(jī)模擬計算，求連續(xù)擲兩顆骰子，點(diǎn)數(shù)之和大于6且第一次擲出的點(diǎn)數(shù)大于第二次擲出點(diǎn)數(shù)的概率。,第二章 隨機(jī)數(shù),隨機(jī)數(shù)的定義及產(chǎn)生方法 偽隨機(jī)數(shù) 產(chǎn)生偽隨機(jī)數(shù)的乘同余方法 產(chǎn)生偽隨機(jī)數(shù)的乘加同余方法 產(chǎn)生偽隨機(jī)數(shù)的其他方法 偽隨機(jī)數(shù)序列的均勻性和獨(dú)立性 作 業(yè),第二章 隨機(jī)數(shù),由具有已知分布

18、的總體中抽取簡單子樣，在蒙特卡羅方法中占有非常重要的地位?？傮w和子樣的關(guān)系，屬于一般和個別的關(guān)系，或者說屬于共性和個性的關(guān)系。由具有已知分布的總體中產(chǎn)生簡單子樣，就是由簡單子樣中若干個性近似地反映總體的共性。 隨機(jī)數(shù)是實(shí)現(xiàn)由已知分布抽樣的基本量，在由已知分布的抽樣過程中，將隨機(jī)數(shù)作為已知量，用適當(dāng)?shù)臄?shù)學(xué)方法可以由它產(chǎn)生具有任意已知分布的簡單子樣。,隨機(jī)數(shù)的定義及產(chǎn)生方法,隨機(jī)數(shù)的定義及性質(zhì) 隨機(jī)數(shù)表 物理方法,隨機(jī)數(shù)的定義及性質(zhì),在連續(xù)型隨機(jī)變量的分布中，最簡單而且最基本的分布是單位均勻分布。由該分布抽取的簡單子樣稱，隨機(jī)數(shù)序列，其中每一個體稱為隨機(jī)數(shù)。 單位均勻分布也稱為0，1上的均勻分布，

19、其分布密度函數(shù)為： 分布函數(shù)為 ：,由于隨機(jī)數(shù)在蒙特卡羅方法中占有極其重要的位置，我們用專門的符號表示。由隨機(jī)數(shù)序列的定義可知，1，2，是相互獨(dú)立且具有相同單位均勻分布的隨機(jī)數(shù)序列。也就是說，獨(dú)立性、均勻性是隨機(jī)數(shù)必備的兩個特點(diǎn)。 隨機(jī)數(shù)具有非常重要的性質(zhì)：對于任意自然數(shù)s，由s個隨機(jī)數(shù)組成的s維空間上的點(diǎn)（n+1，n+2，n+s）在s維空間的單位立方體Gs上均勻分布，即對任意的ai， 如下等式成立：,其中P()表示事件發(fā)生的概率。反之，如果隨機(jī)變量序列1, 2對于任意自然數(shù)s，由s個元素所組成的s維空間上的點(diǎn)（n+1，n+s）在Gs上均勻分布，則它們是隨機(jī)數(shù)序列。 由于隨機(jī)數(shù)在蒙特卡羅方法中

20、所處的特殊地位，它們雖然也屬于由具有已知分布的總體中產(chǎn)生簡單子樣的問題，但就產(chǎn)生方法而言，卻有著本質(zhì)上的差別。,隨機(jī)數(shù)表,為了產(chǎn)生隨機(jī)數(shù)，可以使用隨機(jī)數(shù)表。隨機(jī)數(shù)表是由0，1，9十個數(shù)字組成，每個數(shù)字以0.1的等概率出現(xiàn)，數(shù)字之間相互獨(dú)立。這些數(shù)字序列叫作隨機(jī)數(shù)字序列。如果要得到n位有效數(shù)字的隨機(jī)數(shù)，只需將表中每n個相鄰的隨機(jī)數(shù)字合并在一起，且在最高位的前邊加上小數(shù)點(diǎn)即可。例如，某隨機(jī)數(shù)表的第一行數(shù)字為7634258910，要想得到三位有效數(shù)字的隨機(jī)數(shù)依次為0.763，0.425，0.891。 因為隨機(jī)數(shù)表需在計算機(jī)中占有很大內(nèi)存，而且也難以滿足蒙特卡羅方法對隨機(jī)數(shù)需要量非常大的要求，因此，該

21、方法不適于在計算機(jī)上使用。,物理方法,用物理方法產(chǎn)生隨機(jī)數(shù)的基本原理是：利用某些物理現(xiàn)象，在計算機(jī)上增加些特殊設(shè)備，可以在計算機(jī)上直接產(chǎn)生隨機(jī)數(shù)。這些特殊設(shè)備稱為隨機(jī)數(shù)發(fā)生器。用來作為隨機(jī)數(shù)發(fā)生器的物理源主要有兩種：一種是根據(jù)放射性物質(zhì)的放射性，另一種是利用計算機(jī)的固有噪聲。 一般情況下，任意一個隨機(jī)數(shù)在計算機(jī)內(nèi)總是用二進(jìn)制的數(shù)表示的： 其中i（i=1，2，m）或者為0，或者為1。,因此，利用物理方法在計算機(jī)上產(chǎn)生隨機(jī)數(shù)，就是要產(chǎn)生只取0或1的隨機(jī)數(shù)字序列，數(shù)字之間相互獨(dú)立，每個數(shù)字取0或1的概率均為0.5。 用物理方法產(chǎn)生的隨機(jī)數(shù)序列無法重復(fù)實(shí)現(xiàn)，不能進(jìn)行程序復(fù)算，給驗證結(jié)果帶來很大困難。而

22、且，需要增加隨機(jī)數(shù)發(fā)生器和電路聯(lián)系等附加設(shè)備，費(fèi)用昂貴。因此，該方法也不適合在計算機(jī)上使用。,偽隨機(jī)數(shù),偽隨機(jī)數(shù) 偽隨機(jī)數(shù)存在的兩個問題 偽隨機(jī)數(shù)的周期和最大容量,偽隨機(jī)數(shù),在計算機(jī)上產(chǎn)生隨機(jī)數(shù)最實(shí)用、最常見的方法是數(shù)學(xué)方法，即用如下遞推公式： 產(chǎn)生隨機(jī)數(shù)序列。對于給定的初始值1,2，k，確定n+k，=1，2，。經(jīng)常使用的是k=1的情況，其遞推公式為： 對于給定的初始值1，確定n+1，=,偽隨機(jī)數(shù)存在的兩個問題,用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)，存在兩個問題： 遞推公式和初始值1,2，k確定后，整個隨機(jī)數(shù)序列便被唯一確定。不滿足隨機(jī)數(shù)相互獨(dú)立的要求。 由于隨機(jī)數(shù)序列是由遞推公式確定的，而在計算機(jī)上所能表

23、示的0，1上的數(shù)又是有限的，因此，這種方法產(chǎn)生的隨機(jī)數(shù)序列就不可能不出現(xiàn)無限重復(fù)。一旦出現(xiàn)這樣的n，n (n n )，使得下面等式成立： 隨機(jī)數(shù)序列便出現(xiàn)了周期性的循環(huán)現(xiàn)象。對于k=1的情況，只要有一個隨機(jī)數(shù)重復(fù)，其后面的隨機(jī)數(shù)全部重復(fù)，這與隨機(jī)數(shù)的要求是不相符的。,由于這兩個問題的存在，常稱用數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)為偽隨機(jī)數(shù)。對于以上存在的兩個問題，作如下具體分析。 關(guān)于第一個問題，不能從本質(zhì)上加以改變，但只要遞推公式選得比較好，隨機(jī)數(shù)間的相互獨(dú)立性是可以近似滿足的。至于第二個問題，則不是本質(zhì)的。因為用蒙特卡羅方法解任何具體問題時，所使用的隨機(jī)數(shù)的個數(shù)總是有限的，只要所用隨機(jī)數(shù)的個數(shù)不超過偽隨

24、機(jī)數(shù)序列出現(xiàn)循環(huán)現(xiàn)象時的長度就可以了。 用數(shù)學(xué)方法產(chǎn)生的偽隨機(jī)數(shù)容易在計算機(jī)上得到，可以進(jìn)行復(fù)算，而且不受計算機(jī)型號的限制。因此，這種方法雖然存在著一些問題，但仍然被廣泛地在計算機(jī)上使用，是在計算機(jī)上產(chǎn)生偽隨機(jī)數(shù)的主要方法。,偽隨機(jī)數(shù)的周期和最大容量,發(fā)生周期性循環(huán)現(xiàn)象的偽隨機(jī)數(shù)的個數(shù)稱為偽隨機(jī)數(shù)的周期。對于前面介紹的情況，偽隨機(jī)數(shù)的周期為nn。 從偽隨機(jī)數(shù)序列的初始值開始，到出現(xiàn)循環(huán)現(xiàn)象為止，所產(chǎn)生的偽隨機(jī)數(shù)的個數(shù)稱為偽隨機(jī)數(shù)的最大容量。前面的例子中，偽隨機(jī)數(shù)的最大容量為n 。,產(chǎn)生偽隨機(jī)數(shù)的乘同余方法,乘同余方法是由Lehmer在1951年提出來的，它的一般形式是：對于任一初始值x1，偽隨

25、機(jī)數(shù)序列由下面遞推公式確定： 其中a為常數(shù)。,乘同余方法的最大容量的上界,對于任意正整數(shù)M，根據(jù)數(shù)論中的標(biāo)準(zhǔn)分解定理，總可以分解成如下形式： 其中P0=2，P1， Pr表示不同的奇素數(shù)，0表示非負(fù)整數(shù)，1，r表示正整數(shù)。a無論取什么值，乘同余方法的最大容量的上界為： 的最小公倍數(shù)。其中：,關(guān)于a與x1的取值,如果a與x1滿足如下條件： 對于 ， x1與M互素，則乘同余方法產(chǎn)生的偽隨機(jī)數(shù)序列的最大容量達(dá)到最大可能值(M)。,乘同余方法在計算機(jī)上的使用,為了便于在計算機(jī)上使用，通常取 ：=2s 其中s為計算機(jī)中二進(jìn)制數(shù)的最大可能有效位數(shù) x1= 奇數(shù) a = 52k+1 其中k為使52k+1在計算

26、機(jī)上所能容納的最大整數(shù)，即a為計算機(jī)上所能容納的5的最大奇次冪。一般地，s=32時，a=513；s=48，a=515等。偽隨機(jī)數(shù)序列的最大容量(M)=2s-2 。 乘同余方法是使用的最多、最廣的方法，在計算機(jī)上被廣泛地使用。,產(chǎn)生偽隨機(jī)數(shù)的乘加同余方法,產(chǎn)生偽隨機(jī)數(shù)的乘加同余方法是由Rotenberg于1960年提出來的，由于這個方法有很多優(yōu)點(diǎn)，已成為僅次于乘同余方法產(chǎn)生偽隨機(jī)數(shù)的另一主要方法。 乘加同余方法的一般形式是，對任意初始值x1，偽隨機(jī)數(shù)序列由下面遞推公式確定： 其中a和c為常數(shù)。,乘加同余方法的最大容量,關(guān)于乘加同余方法的最大容量問題，有如下結(jié)論：如果對于正整數(shù)M的所有素數(shù)因子P，

27、下式均成立： 當(dāng)M為4的倍數(shù)時，還有下式成立： c與M互素，則乘加同余方法所產(chǎn)生的偽隨機(jī)數(shù)序列的最大容量達(dá)到最大可能值M。,M，x1，a，c的取值,為了便于在計算機(jī)上使用，通常取 M = 2s 其中s為計算機(jī)中二進(jìn)制數(shù)的最大可能有效位數(shù)。 a = 2b + 1(b2) c = 1 這樣在計算中可以使用移位和指令加法，提高計算速度。,產(chǎn)生偽隨機(jī)數(shù)的其他方法,取中方法 加同余方法,偽隨機(jī)數(shù)序列的均勻性和獨(dú)立性,判斷偽隨機(jī)數(shù)序列是否滿足均勻和相互獨(dú)立的要求，要靠統(tǒng)計檢驗的方法實(shí)現(xiàn)。對于偽隨機(jī)數(shù)的統(tǒng)計檢驗，一般包括兩大類：均勻性檢驗和獨(dú)立性檢驗。 六十年代初，人們開始用定性的方法研究偽隨機(jī)數(shù)序列的均勻

28、性和獨(dú)立性問題，簡要敘述如下。,偽隨機(jī)數(shù)的均勻性,這里只考慮偽隨機(jī)數(shù)序列1,2，n全體作為子樣時的均勻性問題。其中n為偽隨機(jī)數(shù)序列的最大容量。 對于任意的0 x1，令Nn(x)表示偽隨機(jī)數(shù)序列1,2，n中適合不等式 i x i=1，2，n 的個數(shù)，則 標(biāo)志偽隨機(jī)數(shù)序列1,2，n的均勻程度，稱為均勻偏度。,將偽隨機(jī)數(shù)序列1,2，n從小至大重新排列 并令 ，則由(n)的定義，容易證明 很明顯，對于固定的，(n)的值越小越好。它是描述偽隨機(jī)數(shù)序列均勻程度的基本量。對于任意隨機(jī)數(shù)序列，均有如下不等式成立： 當(dāng) 時，所對應(yīng)的偽隨機(jī)數(shù)序列為最佳分布。,可以證明，偽隨機(jī)數(shù)序列為最佳分布的充要條件是它取遍序列

29、 的所有值。 對于計算機(jī)上使用的乘同余方法，按照前面介紹的方法選取a、x1時，所產(chǎn)生的偽隨機(jī)數(shù)序列的均勻偏度 對于乘加同余方法 對于部分偽隨機(jī)數(shù)的均勻性問題通常用統(tǒng)計檢驗方法檢驗。,偽隨機(jī)數(shù)的獨(dú)立性,對于任意 ，令 表示(1,2)， (2,3)， (n,n+1)中適合不等式 的個數(shù)，根據(jù)隨機(jī)變量間相互獨(dú)立的定義和頻率近似概率的方法，令 則(n)標(biāo)志偽隨機(jī)數(shù)序列1,2，n的獨(dú)立程度，簡稱為獨(dú)立偏度。對于固定的n,(n)的值越接近于零，偽隨機(jī)數(shù)序列的獨(dú)立性越好。,對于乘同余方法， 對于乘加同余方法， 因此，這兩種方法的獨(dú)立性都是很好的。 同偽隨機(jī)數(shù)的均勻性問題一樣，偽隨機(jī)數(shù)序列的獨(dú)立性問題也是對它

30、的全體討論的。若只考慮偽隨機(jī)數(shù)的一部分，在通常情況下給出(i)是相當(dāng)因難的。因此，偽隨機(jī)數(shù)序列的獨(dú)立性問題的統(tǒng)計檢驗方法同樣是非常重要的。,作 業(yè),證明1是隨機(jī)數(shù)。 證明 與 同分布 。,第三章 由已知分布的隨機(jī)抽樣,隨機(jī)抽樣及其特點(diǎn) 直接抽樣方法 挑選抽樣方法 復(fù)合抽樣方法 復(fù)合挑選抽樣方法 替換抽樣方法 隨機(jī)抽樣的一般方法 隨機(jī)抽樣的其它方法 作 業(yè),第三章 由已知分布的隨機(jī)抽樣,本章敘述由己知分布抽樣的各主要方法，并給出在粒子輸運(yùn)問題中經(jīng)常用到的具體實(shí)例。,隨機(jī)抽樣及其特點(diǎn),由巳知分布的隨機(jī)抽樣指的是由己知分布的總體中抽取簡單子樣。隨機(jī)數(shù)序列是由單位均勻分布的總體中抽取的簡單子樣，屬于一

31、種特殊的由已知分布的隨機(jī)抽樣問題。本章所敘述的由任意已知分布中抽取簡單子樣，是在假設(shè)隨機(jī)數(shù)為已知量的前提下，使用嚴(yán)格的數(shù)學(xué)方法產(chǎn)生的。 為方便起見，用XF表示由己知分布F(x)中產(chǎn)生的簡單子樣的個體。對于連續(xù)型分布，常用分布密度函數(shù)f(x)表示總體的己知分布，用Xf表示由己知分布密度函數(shù)f(x)產(chǎn)生的簡單子樣的個體。另外，在抽樣過程中用到的偽隨機(jī)數(shù)均稱隨機(jī)數(shù)。,直接抽樣方法,對于任意給定的分布函數(shù)F(x)，直接抽樣方法如下： 其中，1，2，N為隨機(jī)數(shù)序列。為方便起見，將上式簡化為： 若不加特殊說明，今后將總用這種類似的簡化形式表示，總表示隨機(jī)數(shù)。,證明,下面證明用前面介紹的方法所確定的隨機(jī)變量

32、序列X1，X2，XN具有相同分布F(x)。 對于任意的n成立，因此隨機(jī)變量序列X1，X2，XN具有相同分布F(x)。另外，由于隨機(jī)數(shù)序列1，2，N是相互獨(dú)立的，而直接抽樣公式所確定的函數(shù)是波雷爾（Borel）可測的，因此，由它所確定的X1，X2，XN也是相互獨(dú)立的（P.R.Halmos, Measure theory, N.Y.Von Nosrtand,195045定理2）。,離散型分布的直接抽樣方法,對于任意離散型分布： 其中x1，x2，為離散型分布函數(shù)的跳躍點(diǎn)，P1，P2，為相應(yīng)的概率，根據(jù)前述直接抽樣法，有離散型分布的直接抽樣方法如下： 該結(jié)果表明，為了實(shí)現(xiàn)由任意離散型分布的隨機(jī)抽樣，直

33、接抽樣方法是非常理想的。,例1. 二項分布的抽樣,二項分布為離散型分布，其概率函數(shù)為： 其中，P為概率。對該分布的直接抽樣方法如下：,例2. 泊松(Possion)分布的抽樣,泊松(Possion)分布為離散型分布，其概率函數(shù)為： 其中，0 。對該分布的直接抽樣方法如下：,例3. 擲骰子點(diǎn)數(shù)的抽樣,擲骰子點(diǎn)數(shù)X=n的概率為： 選取隨機(jī)數(shù)，如 則 在等概率的情況下，可使用如下更簡單的方法： 其中表示取整數(shù)。,例4. 碰撞核種類的確定,中子或光子在介質(zhì)中發(fā)生碰撞時，如介質(zhì)是由多種元素組成，需要確定碰撞核的種類。假定介質(zhì)中每種核的宏觀總截面分別為1，2，n，則中子或光子與每種核碰撞的概率分別為： 其

34、中t12n。碰撞核種類的確定方法為：產(chǎn)生一個隨機(jī)數(shù)，如果 則中子或光子與第I種核發(fā)生碰撞。,例5. 中子與核的反應(yīng)類型的確定,假設(shè)中子與核的反應(yīng)類型有如下幾種:彈性散射，非彈性散射，裂變，吸收，相應(yīng)的反應(yīng)截面分別為el，in，f，a。則發(fā)生每一種反應(yīng)類型的概率依次為 ： 其中反應(yīng)總截面telinfa。,反應(yīng)類型的確定方法為：產(chǎn)生一個隨機(jī)數(shù),連續(xù)型分布的直接抽樣方法,對于連續(xù)型分布，如果分布函數(shù)F(x) 的反函數(shù) F1(x)存在，則直接抽樣方法是 ：,例6. 在a，b上均勻分布的抽樣,在a，b上均勻分布的分布函數(shù)為： 則,例7. 分布,分布為連續(xù)型分布，作為它的一個特例是： 其分布函數(shù)為： 則,

35、例8. 指數(shù)分布,指數(shù)分布為連續(xù)型分布，其一般形式如下： 其分布函數(shù)為： 則 因為1也是隨機(jī)數(shù)，可將上式簡化為,連續(xù)性分布函數(shù)的直接抽樣方法對于分布函數(shù)的反函數(shù)存在且容易實(shí)現(xiàn)的情況，使用起來是很方便的。但是對于以下幾種情況，直接抽樣法是不合適的。 分布函數(shù)無法用解析形式給出，因而其反函數(shù)也無法給出。 分布函數(shù)可以給出其解析形式，但是反函數(shù)給不出來。 分布函數(shù)即使能夠給出反函數(shù)，但運(yùn)算量很大。 下面敘述的挑選抽樣方法是克服這些困難的比較好的方法。,挑選抽樣方法,為了實(shí)現(xiàn)從己知分布密度函數(shù)f(x)抽樣，選取與f(x)取值范圍相同的分布密度函數(shù)h(x)，如果 則挑選抽樣方法為：,即從h(x)中抽樣x

36、h，以 的概率接受它。 下面證明xf 服從分布密度函數(shù)f(x)。 證明：對于任意x,使用挑選抽樣方法時，要注意以下兩點(diǎn)：選取h(x)時要使得h(x)容易抽樣且M的值要盡量小。因為M小能提高抽樣效率。抽樣效率是指在挑選抽樣方法中進(jìn)行挑選時被選中的概率。按此定義，該方法的抽樣效率E為： 所以，M越小，抽樣效率越高。,當(dāng) f(x) 在0，1上定義時，取 h(x)=1，Xh=， 此時挑選抽樣方法為,例9. 圓內(nèi)均勻分布抽樣,令圓半徑為R0，點(diǎn)到圓心的距離為r，則r的分布密度函數(shù)為 分布函數(shù)為 容易知道，該分布的直接抽樣方法是,由于開方運(yùn)算在計算機(jī)上很費(fèi)時間，該方法不是好方法。下面使用挑選抽樣方法：取

37、則抽樣框圖為,顯然，沒有必要舍棄12的情況，此時，只需取 就可以了，亦即 另一方面，也可證明 與 具有相同的分布 。,復(fù)合抽樣方法,在實(shí)際問題中，經(jīng)常有這樣的隨機(jī)變量，它服從的分布與一個參數(shù)有關(guān)，而該參數(shù)也是一個服從確定分布的隨機(jī)變量，稱這樣的隨機(jī)變量服從復(fù)合分布。例如，分布密度函數(shù) 是一個復(fù)合分布。其中Pn0，n=1，2，且 fn(x)為與參數(shù)n有關(guān)的分布密度函數(shù)，n=1，2， 參數(shù)n服從如下分布,復(fù)合分布的一般形式為： 其中f2(x/y)表示與參數(shù)y有關(guān)的條件分布密度函數(shù)， F1(y)表示分布函數(shù)。 復(fù)合分布的抽樣方法為:首先由分布函數(shù)F1(y) 或分布密度函數(shù)f1(y)中抽樣YF1或Yf

38、1，然后再由分布密度函數(shù)f2(x/ YF1)中抽樣確定Xf2 (x/YF) 證明： 所以，Xf所服從的分布為f (x)。,例10. 指數(shù)函數(shù)分布的抽樣,指數(shù)函數(shù)分布的一般形式為： 引入如下兩個分布密度函數(shù)：,則 使用復(fù)合抽樣方法，首先從f1(y)中抽取y 再由f2(x/ YF1)中抽取x,復(fù)合挑選抽樣方法,考慮另一種形式的復(fù)合分布如下： 其中0H(x,y)M，f2(x/y)表示與參數(shù)y有關(guān)的條件分布密度函數(shù)，F(xiàn)1(y)表示分布函數(shù)。抽樣方法如下：,證明： 抽樣效率為：E=1/M,為了實(shí)現(xiàn)某個復(fù)雜的隨機(jī)變量 y 的抽樣，將其表示成若干個簡單的隨機(jī)變量 x1，x2，xn 的函數(shù) 得到 x1，x2，

39、xn 的抽樣后，即可確定 y 的抽樣，這種方法叫作替換法抽樣。即,替換抽樣方法,例11. 散射方位角余弦分布的抽樣,散射方位角在0,2上均勻分布，則其正弦和余弦sin和cos服從如下分布： 直接抽樣方法為：,令=2，則在0,上均勻分布，作變換 其中01，0，則 (x,y) 表示上半個單位圓內(nèi)的點(diǎn)。如果 (x,y) 在上半個單位圓內(nèi)均勻分布，則在0,上均勻分布，由于,因此抽樣sin和cos的問題就變成在上半個單位圓內(nèi)均勻抽樣 (x,y) 的問題。 為獲得上半個單位圓內(nèi) 的均勻點(diǎn)，采用挑選法，在 上半個單位圓的外切矩形內(nèi) 均勻投點(diǎn)（如圖）。 舍棄圓外的點(diǎn)，余下的就是所要求的點(diǎn)。 抽樣方法為： 抽樣

40、效率 E=/40.785,為實(shí)現(xiàn)散射方位角余弦分布抽樣，最重要的是在上半個單位圓內(nèi)產(chǎn)生均勻分布點(diǎn)。下面這種方法，首先在單位圓的半個外切正六邊形內(nèi)產(chǎn)生均勻分布點(diǎn)，如圖所示。,于是便有了抽樣效率更高的抽樣方法： 抽樣效率,例12. 正態(tài)分布的抽樣,標(biāo)準(zhǔn)正態(tài)分布密度函數(shù)為： 引入一個與標(biāo)準(zhǔn)正態(tài)隨機(jī)變量X獨(dú)立同分布的隨機(jī)變量Y，則（X,Y）的聯(lián)合分布密度為： 作變換,則（,）的聯(lián)合分布密度函數(shù)為： 由此可知，與相互獨(dú)立，其分布密度函數(shù)分別為 分別抽取， ：,從而得到一對服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量X和Y： 對于一般的正態(tài)分布密度函數(shù) N(,2) 的抽樣，其抽樣結(jié)果為：,例13. 分布的抽樣,分布密度函數(shù)

41、的一般形式為： 其中n，k為整數(shù)。為了實(shí)現(xiàn)分布的抽樣，將其看作一組簡單的相互獨(dú)立隨機(jī)變量的函數(shù)，通過這些簡單隨機(jī)變量的抽樣，實(shí)現(xiàn)分布的抽樣。設(shè) x1，x2，xn 為一組相互獨(dú)立、具有相同分布 F(x) 的隨機(jī)變量，k為 x1，x2，xn 按大小順序排列后的第k個，記為：,則k的分布函數(shù)為： 當(dāng) F(x)=x 時， 不難驗證，k的分布密度函數(shù)為分布。因此， 分布的抽樣可用如下方法實(shí)現(xiàn)： 選取n個隨機(jī)數(shù)，按大小順序排列后取第k個，即,隨機(jī)抽樣的一般方法,加抽樣方法 減抽樣方法 乘抽樣方法 乘加抽樣方法 乘減抽樣方法 對稱抽樣方法 積分抽樣方法,加抽樣方法,加抽樣方法是對如下加分布給出的一種抽樣方法

42、： 其中Pn0， ，且 fn(x)為與參數(shù)n有關(guān)的分布密度函數(shù)，n=1，2，。 由復(fù)合分布抽樣方法可知，加分布的抽樣方法為:首先抽樣確定n，然后由 fn(x)中抽樣x，即：,例14. 多項式分布抽樣,多項式分布密度函數(shù)的一般形式為： 將 f(x) 改寫成如下形式： 則該分布的抽樣方法為：,例15. 球殼內(nèi)均勻分布抽樣,設(shè)球殼內(nèi)半徑為R0，外半徑為R1，點(diǎn)到球心的距離為r，則r的分布密度函數(shù)為 分布函數(shù)為 該分布的直接抽樣方法是,為避免開立方根運(yùn)算，作變換： 則 x0,1，其分布密度函數(shù)為： 其中,則x及r的抽樣方法為：,減抽樣方法,減抽樣方法是對如下形式的分布密度所給出的一種抽樣方法： 其中A

43、1、A2為非負(fù)實(shí)數(shù)，f1(x) 、f2(x)均為分布密度函數(shù)。 減抽樣方法分為以下兩種形式：,以上兩種形式的抽樣方法，究竟選擇哪種好，要看f1(x) 、f2(x)哪一個容易抽樣，如相差不多，選用第一種方法抽樣效率高。,（1）將f (x)表示為 令m表示f2(x)f1(x)的下界，使用挑選法，從f1(x)中抽取Xf1 抽樣效率為：,（2）將f (x)表示為 使用挑選法，從f2(x)中抽取Xf2 抽樣效率為：,例16. 分布抽樣,分布的一個特例： 取A12，A21，f1(x)1，f2(x)2x，此時m0，則根據(jù)第一種形式的減抽樣方法，有 或,由于11可用1代替，該抽樣方法可簡化為： 對于21的情況

44、，可取 Xf1 ，因此 與分布的推論相同。,如下形式的分布稱為乘分布： 其中H(x)為非負(fù)函數(shù)， f1(x)為任意分布密度函數(shù)。 令M為H(x)的上界，乘抽樣方法如下： 抽樣效率為：,乘抽樣方法,例17. 倒數(shù)分布抽樣,倒數(shù)分布密度函數(shù)為： 其直接抽樣方法為： 下面采用乘抽樣方法，考慮如下分布族： 其中 i = 1，2，該分布的直接抽樣方法為：,利用這一分布族，將倒數(shù)分布 f(x) 表示成: 其中， 乘法分布的抽樣方法如下: 該分布的抽樣效率為：,例18. 麥克斯韋(Maxwell)分布抽樣,麥克斯韋分布密度函數(shù)的一般形式為： 使用乘抽樣方法，令 該分布的直接抽樣方法為：,此時 則麥克斯韋分布

45、的抽樣方法為: 該分布的抽樣效率為：,在實(shí)際問題中，經(jīng)常會遇到如下形式的分布： 其中Hn(x)為非負(fù)函數(shù)，fn(x) 為任意分布密度函數(shù)，n=1，2，。不失一般性，只考慮n=2的情況： 將 f(x) 改寫成如下的加分布形式：,乘加抽樣方法,其中,乘加抽樣方法為： 該方法的抽樣效率為：,這種方法需要知道P1的值(P2=1P1），這對有些分布是很困難的。下面的方法可以不用計算P1 ： 對于任意小于1的正數(shù)P1 ，令P2=1P1 ；,則采用復(fù)合挑選抽樣方法，有：,當(dāng)取 時，抽樣效率最高 這時，乘加抽樣方法為：,由于 可知第一種方法比第二種方法的抽樣效率高。,例19. 光子散射后能量分布的抽樣,令光子

46、散射前后的能量分別為 和 （以 m0c2 為單位，m0為電子靜止質(zhì)量，c 為光速）， ， 則 x 的分布密度函數(shù)為： 該分布即為光子散射能量分布，它是由著名的KlinNishina 公式確定的。其中 K() 為歸一因子：,把光子散射能量分布改寫成如下形式： 在1, 1+2上定義如下函數(shù):,則有 使用乘加抽樣方法：,光子散射能量分布的抽樣方法為： 該方法的抽樣效率為：,乘減分布的形式為： 其中H1(x) 、H2(x)為非負(fù)函數(shù)，f1(x)、f2(x) 為任意分布密度函數(shù)。 與減抽樣方法類似，乘減分布的抽樣方法也分為兩種。,乘減抽樣方法,（1）將 f (x) 表示為 令H1(x)的上界為M1， 的

47、下界為m，使用乘抽 樣方法得到如下乘減抽樣方法：,（2）將 f (x) 表示為 令H2(x)的上界為M2，使用乘抽樣方法，得到另一種乘減抽樣方法：,例20. 裂變中子譜分布抽樣,裂變中子譜分布的一般形式為： 其中A，B，C，Emin，Emax 均為與元素有關(guān)的量。令 其中為歸一因子，為任意參數(shù)。,相應(yīng)的 H1(E)，H2(E) 為： 于是裂變中子譜分布可以表示成乘減分布形式: 容易確定 H1(E) 的上界為： 為提高抽樣效率，應(yīng)取使得M1 達(dá)到最小，此時,取 m0，令 則裂變中子譜分布的抽樣方法為: 抽樣效率,對稱分布的一般形式為： 其中 f1(x) 為任意分布密度函數(shù)，滿足偶函數(shù)對稱條件，H

48、(x) 為任意奇函數(shù)，即對任意x滿足： 對稱分布的抽樣方法如下：取=21,對稱抽樣方法,證明： 因為=21，x 相當(dāng)于 ，因此,例21. 質(zhì)心系各向同性散射角余弦分布抽樣,在質(zhì)心系各向同性散射的假設(shè)下，為得到實(shí)驗室系散射角余弦，需首先抽樣確定質(zhì)心條散射角余弦： 再利用下面轉(zhuǎn)換公式: 得到實(shí)驗室系散射角余弦L。其中A為碰撞核質(zhì)量，C、L 分別為質(zhì)心系和實(shí)驗室系散射角。,為避免開方運(yùn)算，可以使用對稱分布抽樣。 根據(jù)轉(zhuǎn)換公式可得： 依照質(zhì)心系散射各向同性的假定，可得到實(shí)驗室系散射角余弦L 的分布如下： 該密度函數(shù)中的第一項為偶函數(shù)，第二項為奇函數(shù)，因而是對稱分布。其中,從 f1(L) 的抽樣可使用挑

49、選法 然后再以 的概率決定接受或取負(fù)值。 上述公式涉及開方運(yùn)算，需要進(jìn)一步簡化。,注意以下事實(shí)：對于任意0a1 令 則上述挑選抽樣中的挑選條件簡化為： 另一方面，在 即 的條件下，2/a 在1, 1上均勻分布，故可令2/a，則最終決定取正負(fù)值的條件簡化為：,于是，得到質(zhì)心系各向同性散射角余弦分布的抽樣方法為：,如下形式的分布密度函數(shù) 稱為積分分布密度函數(shù)，其中 f0(x，y) 為任意二維分布密度函數(shù)，H(x)為任意函數(shù)。該分布密度函數(shù)的抽樣方法為：,積分抽樣方法,證明：對于任意x,例22. 各向同性散射方向的抽樣,為了確定各向同性散射方向 ，根據(jù)公式： 對于各向同性散射，cos在1, 1上均勻

50、分布，在0, 2上均勻分布。由于 直接抽樣需要計算三角函數(shù)和開方。,定義兩個隨機(jī)變量： 可以證明，當(dāng) 時，隨機(jī)變量 x 和 y 服從如下分布: 定義區(qū)域為：,則 wcos 的分布可以用上述分布表示成積分分布的形式： 令 ，則屬于上述積分限內(nèi)的 y 一定滿足 條件 。,各向同性散射方向的抽樣方法為： 抽樣效率為：,隨機(jī)抽樣的其它方法,偏倚抽樣方法 近似抽樣方法 近似-修正抽樣方法 多維分布抽樣方法 指數(shù)分布的抽樣,使用蒙特卡羅方法計算積分 時，可考慮將積分I改寫為 其中 f *(x) 為一個與 f (x) 有相同定義域的新的分布密度函數(shù)。于是可以這樣計算積分I： 這里 Xi 是從 f *(x)

51、中抽取的第 i 個子樣。,偏移抽樣方法,由此可以看出，原來由 f (x) 抽樣，現(xiàn)改為由另一個分布密度函數(shù) f *(x) 抽樣，并附帶一個權(quán)重糾偏因子 這種方法稱為偏倚抽樣方法。 從 f (x) 中抽取的 Xf ，滿足 而對于偏倚抽樣，有 一般情況下，Xf 是具有分布 f (x) 總體的簡單子樣的個體，只代表一個。Xf* 是具有分布 f *(x) 總體的簡單子樣的個體，但不代表一個，而是代表 W(Xf*) 個，這時Xf*是帶權(quán)W(Xf*)服從分布 f (x) 。,在實(shí)際問題中，分布密度函數(shù)的形式有時是非常復(fù)雜的，有些甚至不能用解析形式給出，只能用數(shù)據(jù)或曲線形式給出。如中子散射角余弦分布多數(shù)是以

52、曲線形式給出的。對于這樣的分布，需要用近似分布密度函數(shù)代替原來的分布密度函數(shù)，用近似分布密度函數(shù)的抽樣代替原分布密度函數(shù)的抽樣，這種方法稱為近似抽樣方法。,近似抽樣方法,設(shè) fa(x) f (x)，即 fa(x) 是 f (x) 的一個近似分布密度函數(shù)。對于階梯近似，有 其中，x0，x1， ，xn為任意分點(diǎn)。在此情況下，近似抽樣方法為： 或,階梯近似,對于梯形近似，有 其中，c 為歸一因子， fi f (xi) ，x0，x1， ，xn為任意分點(diǎn)。根據(jù)對稱抽樣方法，梯形近似抽樣方法為：,梯形近似,除了上述這種近似外，近似抽樣方法還包括對直接抽樣方法中分布函數(shù)反函數(shù)的近似處理，以及用具有近似分布的

53、隨機(jī)變量代替原分布的隨機(jī)變量。,例23. 正態(tài)分布的近似抽樣,我們知道，隨機(jī)數(shù)的期望值為 1/2，方差為 1/12，則隨機(jī)變量 漸近正態(tài)分布，因此，當(dāng) n 足夠大時便可用 Xn 作為正態(tài)分布的近似抽樣。特別是 n12 時，有,對于任意分布密度函數(shù) f (x) ，設(shè) fa(x) 是 f (x) 的一個近似分布密度函數(shù)，它的特點(diǎn)是抽樣簡單，運(yùn)算量小。令 則分布密度函數(shù) f(x) 可以表示為乘加分布形式： 其中 H1(x) 為非負(fù)函數(shù)，f1(x) 為一分布密度函數(shù)。 對 f(x) 而言，fa(x) 是它的近似分布密度函數(shù)，而H1(x) f1(x)正好是這種近似的修正。,近似-修正抽樣方法,近似-修正

54、抽樣方法如下： 抽樣效率 由上述近似-修正抽樣方法可以看出，如果近似分布密度函數(shù) fa(x) 選得好，m 接近 1，這時有很大可能直接從 fa(x) 中抽取 Xfa ，而只有很少的情況需要計算與f (x) 有關(guān)的函數(shù) H1(Xf1)。在乘抽樣方法中，每一次都要計算 H(Xfa)f (Xfa)fa(Xfa)。因此，當(dāng) f (x) 比較復(fù)雜時，近似-修正抽樣方法有很大好處。,例24. 裂變中子譜分布的近似-修正抽樣,裂變中子譜分布的一般形式為： 其中A，B，C，Emin，Emax 均為與元素有關(guān)的量。 對于鈾-235， A=0.965，B=2.29，C=0.453，Emin=0，Emax=。 若采

55、用乘減抽樣方法，其抽樣效率約為0.5。,令 相應(yīng)的 則 從 fa(x) 的抽樣為 從 f1(x) 的抽樣為,參數(shù)的確定，使1A0，且使 H1(E) 的上界M1 最小。裂變中子譜的近似修正抽樣方法為 對于鈾-235，m0.8746，M0.2678，0.5543，抽樣效率 E0.9333。而且近似修正抽樣方法有0.8746的概率直接用近似分布抽樣，只計算一次對數(shù)。因此，較之乘減抽樣方法大大節(jié)省了計算時間，提高了抽樣效率。,為方便起見，這里僅討論二維分布的情況，對于更高維數(shù)的分布，可用類似的方法處理。 對于任意二維分布密度函數(shù)，總可以用其邊緣分布密度函數(shù)和條件分布密度函數(shù)的乘積表示: 其中 fl(x

56、)，f2(y|x) 分別為分布 f (x,y) 的邊緣分布密度函數(shù)和條件分布密度函數(shù)，即,多維分布抽樣方法,二維分布密度函數(shù)的抽樣方法是: 首先由 fl(x) 中抽取 Xf1，再由 f2(y|Xf1) 中抽樣確定 Yf2 。 對于多維分布密度函數(shù)，也可直接采用類似于一維分布密度函數(shù)的抽樣方法。例如，對如下形式的二維分布密度函數(shù)： 其中 H(x,y) 為非負(fù)函數(shù)，f1(x,y) 為任意二維分布密度函數(shù)。設(shè) M 為 H(x,y) 的上界，則有二維分布的乘抽樣方法如下：,例25. 下面二維分布密度函數(shù)的抽樣,將 f (x,y) 寫為 其中 用直接抽樣方法分別從 fl(x) 和 f2(y|Xf1) 中抽樣，得到,前面已經(jīng)介紹了，指數(shù)分布 的直接抽樣為： 這不僅需要計算對數(shù)，而且由于要使用偽隨機(jī)數(shù)，受精度的限制，該抽樣值在小概率處即數(shù)值較大處呈現(xiàn)明顯得離散性。 下面介紹兩種抽樣方法可以避免這些問題。,指數(shù)分布的抽樣,所用隨機(jī)數(shù)的平均個數(shù) Ne2 / ( e1)4.3,方法一,