蒙特卡洛應(yīng)用講義 完整版

蒙特卡羅方法 在核技術(shù)中的應(yīng)用 林謙 目 錄 第一章 蒙特卡羅方法概述 第二章 隨機(jī)數(shù) 第三章 由已知分布的隨機(jī)抽樣 第四章 蒙特卡羅方法解粒子輸運(yùn)問題 教材 蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用 許淑艷 編著 原子能出版社 蒙特卡羅方法 清華大學(xué) 參考書 蒙特卡羅方法及其在粒子輸運(yùn)問題中的應(yīng)用 裴鹿成 張孝澤 編著 科學(xué)出版社 蒙特卡羅方法 徐鐘濟(jì) 編著 上海科學(xué)技術(shù)出版社 聯(lián)系方式 電話 83918 電子郵件 第一章 蒙特卡羅方法概述 蒙特卡羅方法的基本思想 蒙特卡羅方法的收斂性，誤差 蒙特卡羅方法的特點(diǎn) 蒙特卡羅方法的主要應(yīng)用范圍 作 業(yè) 第一章 蒙特卡羅方法概述 蒙特卡羅方法 又稱 隨機(jī)抽樣技巧 或 統(tǒng)計(jì)試驗(yàn)方法 。半個(gè)多世紀(jì)以來(lái)，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明 ，這種方法作為一種獨(dú)立的方法被提出來(lái)，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。蒙特卡羅方法是一種計(jì)算方法，但與一般數(shù)值計(jì)算方法有很大區(qū)別。它是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。由于蒙特卡羅方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過程，解決一些數(shù)值方法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。 1. 蒙特卡羅方法的基本思想 二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨(dú)立的方法被提出來(lái)，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用。但其基本思想并非新穎，人們?cè)谏a(chǎn)實(shí)踐和科學(xué)試驗(yàn)中就已發(fā)現(xiàn)，并加以利用。 兩個(gè)例子 例 1. 蒲豐氏問題 例 2. 射擊問題（打靶游戲） 基本思想 計(jì)算機(jī)模擬試驗(yàn)過程 例 1. 蒲豐氏問題 為了求得圓周率 值，在十九世紀(jì)后期，有很多人作了這樣的試驗(yàn)：將長(zhǎng)為 2l的一根針任意投到地面上，用針與一組相間距離為 2a（ l a）的平行線相交的頻率代替概率 P，再利用準(zhǔn)確的關(guān)系式： 求出 值 其中 為投計(jì)次數(shù)， n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。 alP2)(22 nNa laP l 一些人進(jìn)行了實(shí)驗(yàn)，其結(jié)果列于下表 ： 實(shí)驗(yàn)者 年份 投計(jì)次數(shù) 的實(shí)驗(yàn)值 沃爾弗 (Wolf) 1850 5000 3.1596 斯密思 (Smith) 1855 3204 3.1553 福克斯 (Fox) 1894 1120 3.1419 拉查里尼(Lazzarini) 1901 3408 3.1415929 例 2. 射擊問題（打靶游戲） 設(shè) r表示射擊運(yùn)動(dòng)員的彈著點(diǎn)到靶心的距離， (r)表示擊中 r處相應(yīng)的得分?jǐn)?shù)（環(huán)數(shù)）， f(r)為該運(yùn)動(dòng)員的彈著點(diǎn)的分布密度函數(shù)，它反映運(yùn)動(dòng)員的射擊水平。該運(yùn)動(dòng)員的射擊成績(jī)?yōu)?用概率語(yǔ)言來(lái)說， 是隨機(jī)變量 (r)的數(shù)學(xué)期望，即 )( rgEg 0 )()( drrfrgg 現(xiàn)假設(shè)該運(yùn)動(dòng)員進(jìn)行了 次射擊 ， 每次射擊的彈著點(diǎn)依次為 r1 ， r2 ， ， rN ， 則 次得分 g(r1) ，g(r2)， ， g(rN)的算術(shù)平均值 代表了該運(yùn)動(dòng)員的成績(jī) 。 換言之 ， 為積分 的估計(jì)值 ， 或近似值 。 在該例中，用 次試驗(yàn)所得成績(jī)的算術(shù)平均值作為數(shù)學(xué)期望 的估計(jì)值（積分近似值）。 NiiN rgNg1)(1 基本思想 由以上兩個(gè)例子可以看出 ， 當(dāng)所求問題的解是某個(gè)事件的概率 ， 或者是某個(gè)隨機(jī)變量的數(shù)學(xué)期望 ， 或者是與概率 、 數(shù)學(xué)期望有關(guān)的量時(shí) ， 通過某種試驗(yàn)的方法 ， 得出該事件發(fā)生的頻率 ， 或者該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值 ， 通過它得到問題的解 。這就是蒙特卡羅方法的基本思想 。 當(dāng)隨機(jī)變量的取值僅為 1或 0時(shí) ， 它的數(shù)學(xué)期望就是某個(gè)事件的概率 。 或者說 ， 某種事件的概率也是隨機(jī)變量 （ 僅取值為 1或 0） 的數(shù)學(xué)期望 。 因此，可以通俗地說，蒙特卡羅方法是用隨機(jī)試驗(yàn)的方法計(jì)算積分，即將所要計(jì)算的積分看作服從某種分布密度函數(shù) f(r)的隨機(jī)變量 (r)的數(shù)學(xué)期望 通過某種試驗(yàn) ， 得到 個(gè)觀察值 r1， r2， ， rN（ 用概率語(yǔ)言來(lái)說 ， 從分布密度函數(shù) f(r)中抽取 個(gè)子樣 r1，r2， ， rN， ） ， 將相應(yīng)的 個(gè)隨機(jī)變量的值 g(r1)，g(r2)， ， g(rN)的算術(shù)平均值 作為積分的估計(jì)值（近似值）。 NiiN rgNg1)(1 0 )()( drrfrgg 為了得到具有一定精確度的近似解，所需試驗(yàn)的次數(shù)是很多的，通過人工方法作大量的試驗(yàn)相當(dāng)困難，甚至是不可能的。因此，蒙特卡羅方法的基本思想雖然早已被人們提出，卻很少被使用。本世紀(jì)四十年代以來(lái)，由于電子計(jì)算機(jī)的出現(xiàn)，使得人們可以通過電子計(jì)算機(jī)來(lái)模擬隨機(jī)試驗(yàn)過程，把巨大數(shù)目的隨機(jī)試驗(yàn)交由計(jì)算機(jī)完成，使得蒙特卡羅方法得以廣泛地應(yīng)用，在現(xiàn)代化的科學(xué)技術(shù)中發(fā)揮應(yīng)有的作用。 計(jì)算機(jī)模擬試驗(yàn)過程 計(jì)算機(jī)模擬試驗(yàn)過程 ， 就是將試驗(yàn)過程 （ 如投針 ，射擊 ） 化為數(shù)學(xué)問題 ， 在計(jì)算機(jī)上實(shí)現(xiàn) 。 以上述兩個(gè)問題為例 ， 分別加以說明 。 例 1. 蒲豐氏問題 例 2. 射擊問題（打靶游戲） 由上面兩個(gè)例題看出 ， 蒙特卡羅方法常以一個(gè)“ 概率模型 ” 為基礎(chǔ) ， 按照它所描述的過程 ， 使用由已知分布抽樣的方法 ， 得到部分試驗(yàn)結(jié)果的觀察值 ，求得問題的近似解 。 例蒲豐氏問題 設(shè)針投到地面上的位置可以用一組參數(shù)（ x,）來(lái)描述， x為針中心的坐標(biāo)， 為針與平行線的夾角，如圖所示。 任意投針，就是意味著 x與都是任意取的，但 x的范圍限于 0， a，夾角 的范圍限于 0， 。在此情況下，針與平行線相交的數(shù)學(xué)條件是 針在平行線間的位置 s in lx 如何產(chǎn)生任意的（ x,）？x在 0， a上任意取值，表示x在 0， a上是均勻分布的，其分布密度函數(shù)為： 類似地， 的分布密度函數(shù)為： 因此，產(chǎn)生任意的（ x,）的過程就變成了由 f1(x)抽樣 x及由 f2()抽樣 的過程了。由此得到： 其中 1， 2均為（ 0,1）上均勻分布的隨機(jī)變量。 其他,00,/1)(1axaxf 其他,00,/1)(2f21 ax 每次投針試驗(yàn) ， 實(shí)際上變成在計(jì)算機(jī)上從兩個(gè)均勻分布的隨機(jī)變量中抽樣得到 （ x,） ， 然后定義描述針與平行線相交狀況的隨機(jī)變量 s(x,)， 為 如果投針 次 ， 則 是針與平行線相交概率 的估計(jì)值。事實(shí)上， 于是有 其他當(dāng),0s i n,1),(lxxsNiiiN xsNs1),(1 aladxdd xdfxfxsPl 2)()(),(s i n0021NsalaPl 22 例射擊問題 設(shè)射擊運(yùn)動(dòng)員的彈著點(diǎn)分布為 用計(jì)算機(jī)作隨機(jī)試驗(yàn)（射擊）的方法為，選取一個(gè)隨機(jī)數(shù) ，按右邊所列方法判斷得到成績(jī)。 這樣，就進(jìn)行了一次隨機(jī)試驗(yàn)（射擊），得到了一次成績(jī) (r)，作 次試驗(yàn)后，得到該運(yùn)動(dòng)員射擊成績(jī)的近似值 環(huán)數(shù) 7 8 9 10 概率 0.1 0.1 0.3 0.5 環(huán)中命環(huán)命中環(huán)命中環(huán)命中1095.082.071.0 NiiN rgNg1)(12. 蒙特卡羅方法的收斂性，誤差 蒙特卡羅方法作為一種計(jì)算方法 ， 其收斂性與誤差是普遍關(guān)心的一個(gè)重要問題 。 收斂性 誤差 減小方差的各種技巧 效率 收斂性 由前面介紹可知，蒙特卡羅方法是由隨機(jī)變量 X的簡(jiǎn)單子樣 X1， X2， ， XN的算術(shù)平均值： 作為所求解的近似值。由大數(shù)定律可知， 如 X1， X2， ， XN獨(dú)立同分布，且具有有限期望值（ E(X)0 。對(duì)該分布的直接抽樣方法如下： !)(nePnxPnn n0ii1n0iii!i!,當(dāng) enX F例 3. 擲骰子點(diǎn)數(shù)的抽樣 擲骰子點(diǎn)數(shù) X=n的概率為： 選取隨機(jī)數(shù) ，如 則 在等概率的情況下，可使用如下更簡(jiǎn)單的方法： 其中表示取整數(shù)。 61)( nXP661 nn nX F 16 FX例 4. 碰撞核種類的確定 中子或光子在介質(zhì)中發(fā)生碰撞時(shí)，如介質(zhì)是由多種元素組成，需要確定碰撞核的種類。假定介質(zhì)中每種核的宏觀總截面分別為 1， 2， ， n，則中子或光子與每種核碰撞的概率分別為： 其中 t 1 2 n。碰撞核種類的確定方法為：產(chǎn)生一個(gè)隨機(jī)數(shù) ，如果 則中子或光子與第 I種核發(fā)生碰撞。 niPtii ,2,1 I1ii1I1ii PP例 5. 中子與核的反應(yīng)類型的確定 假設(shè)中子與核的反應(yīng)類型有如下幾種 :彈性散射，非彈性散射，裂變，吸收，相應(yīng)的反應(yīng)截面分別為 el，in， f， a。則發(fā)生每一種反應(yīng)類型的概率依次為 ： 其中反應(yīng)總截面 t el in f a。 taatfftinintelel PPPP 反應(yīng)類型的確定方法為：產(chǎn)生一個(gè)隨機(jī)數(shù) 收吸裂變非彈性散射彈性散射 finelinelelPPPPPP2) 連續(xù)型分布的直接抽樣方法 對(duì)于連續(xù)型分布 ， 如果分布函數(shù) F(x) 的反函數(shù) F 1(x)存在 ， 則直接抽樣方法是 ： )(1 FX F例 6. 在 a， b上均勻分布的抽樣 在 a， b上均勻分布的分布函數(shù)為： 則 bxbxaabaxaxxF當(dāng)當(dāng)當(dāng)10)( )( abaX F例 7. 分布 分布為連續(xù)型分布，作為它的一個(gè)特例是： 其分布函數(shù)為： 則 FX10,2)()( 20 xxtdtdttfxF xx10,2)( xxxf例 8. 指數(shù)分布 指數(shù)分布為連續(xù)型分布，其一般形式如下： 其分布函數(shù)為： 則 因?yàn)?1 也是隨機(jī)數(shù)，可將上式簡(jiǎn)化為 0,1)()( 0 xedteadttfxF axx atx0,)( xeaxf ax)1ln (1 aX Fln1aX F 連續(xù)性分布函數(shù)的直接抽樣方法對(duì)于分布函數(shù)的反函數(shù)存在且容易實(shí)現(xiàn)的情況 ， 使用起來(lái)是很方便的 。但是對(duì)于以下幾種情況 ， 直接抽樣法是不合適的 。 分布函數(shù)無(wú)法用解析形式給出 ， 因而其反函數(shù)也無(wú)法給出 。 分布函數(shù)可以給出其解析形式 ， 但是反函數(shù)給不出來(lái) 。 分布函數(shù)即使能夠給出反函數(shù) ， 但運(yùn)算量很大 。 下面敘述的挑選抽樣方法是克服這些困難的比較好的方法 。 3. 挑選抽樣方法 為了實(shí)現(xiàn)從己知分布密度函數(shù) f(x)抽樣 ， 選取與f(x)取值范圍相同的分布密度函數(shù) h(x)， 如果 則挑選抽樣方法為： )()(s u pxhxfMxhfhhXXXhMXf)()( 即從 h(x)中抽樣 xh， 以 的概率接受它 。 下面證明 xf 服從分布密度函數(shù) f(x)。 證明：對(duì)于任意 x )()(hhxhM xf)()()()(,)()()(hhhhhhhhfXhMXfPXhMXfdxxXxPXhMXfdxxXxPdxxXxP dxxfdXXfdXXfdXXhXhMXfdXXhXhMXfddXXhddXXhhhdxxxhhhhhhdxxxhhhhXhMXfhhdxxxXhMXfhhhhhh)()()()()()()()()()()()()(0)()(0 使用挑選抽樣方法時(shí) ， 要注意以下兩點(diǎn)：選取 h(x)時(shí)要使得 h(x)容易抽樣且 M的值要盡量小 。 因?yàn)?M小能提高抽樣效率 。 抽樣效率是指在挑選抽樣方法中進(jìn)行挑選時(shí)被選中的概率 。 按此定義 ， 該方法的抽樣效率 E為： 所以 ， M越小 ， 抽樣效率越高 。 MdXXhXhMXfXhMXfPEhhhhhh1)()()()()( 當(dāng) f(x) 在 0， 1 上定義時(shí) ， 取 h(x)=1， Xh= ， 此時(shí)挑選抽樣方法為 )(s u p10xfMx fXMf )(例 9. 圓內(nèi)均勻分布抽 樣 令圓半徑為 R0，點(diǎn)到圓心的距離為 r，則 r的分布密度函數(shù)為 分布函數(shù)為 容易知道，該分布的直接抽樣方法是 其它當(dāng)002)( 020RrRrrf202)(RrrF 0Rr f 由于開方運(yùn)算在計(jì)算機(jī)上很費(fèi)時(shí)間 ， 該方法不是好方法 。 下面使用挑選抽樣方法：取 則抽樣框圖為 00022)()(1)( RrMRrrhrfRrh h， 2021Rr f 顯然 ， 沒有必要舍棄 1 2的情況 ， 此時(shí) ， 只需取 就可以了 ， 亦即 另一方面，也可證明 與 具有相同的分布 。 10 Rr f),m ax ( 210 Rr f ),m a x( 21 2)( rrF 4. 復(fù)合抽樣方法 在實(shí)際問題中 ， 經(jīng)常有這樣的隨機(jī)變量 ， 它服從的分布與一個(gè)參數(shù)有關(guān) ， 而該參數(shù)也是一個(gè)服從確定分布的隨機(jī)變量 ， 稱這樣的隨機(jī)變量服從復(fù)合分布 。例如 ， 分布密度函數(shù) 是一個(gè)復(fù)合分布 。 其中 Pn0， n=1， 2， ， 且 fn(x)為與參數(shù) n有關(guān)的分布密度函數(shù) ， n=1， 2， ， 參數(shù) n服從如下分布 1)()(nnn xfPxf11nnPynnPyF )( 復(fù)合分布的一般形式為： 其中 f2(x/y)表示與參數(shù) y有關(guān)的條件分布密度函數(shù) ， F1(y)表示分布函數(shù) 。 復(fù)合分布的抽樣方法為 :首先由分布函數(shù) F1(y) 或分布密度函數(shù) f1(y)中抽樣 YF1或 Yf1，然后再由分布密度函數(shù) f2(x/ YF1)中抽樣確定 Xf2 (x/YF) 證明： 所以 ， Xf所服從的分布為 f (x)。 )()()( 12 ydFyxfxf )/( 12 FYxff XX dxxfYd x dFYxfdxxXxpdxxXxpFYxff)()()()()(12)/( 12例 10. 指數(shù)函數(shù)分布的抽 樣 指數(shù)函數(shù)分布的一般形式為： 引入如下兩個(gè)分布密度函數(shù)： 其它當(dāng)00)( 1xdyyenxE nxyn 其它當(dāng)其它當(dāng)00)(01)(211xeyyxfyynyfxyn 則 使用復(fù)合抽樣方法 ， 首先從 f1(y)中抽取 y 再由 f2(x/ YF1)中抽取 x 1 12 )()()( dyyfyxfxE n),m a x (11211nnfY 1211ln),m a x (ln1nnfnfYX5. 復(fù)合挑選抽樣方法 考慮另一種形式的復(fù)合分布如下： 其中 0H(x,y)M， f2(x/y)表示與參數(shù) y有關(guān)的條件分布密度函數(shù) ， F1(y)表示分布函數(shù) 。 抽樣方法如下： )()(),()( 12 ydFyxfyxHxf )/()/(12112),(FFYxffFYxfXXMYXH 證明： 抽樣效率為： E=1/M dxxfdxydFyxfyxHyd x d FyxfMyxHyd x d FyxfMyxHdyd x d Fyxfdyd x d FyxfMYXHPMYXHdxxXxPMYXHdxxXxPdxxXxPdxxxMyxHdxxxMyxHFfFffFfff)()()(),()()(),()()(),()()()()(),(),(,),()(121212),(012),(01212122122 為了實(shí)現(xiàn)某個(gè)復(fù)雜的隨機(jī)變量 y 的抽樣，將其表示成若干個(gè)簡(jiǎn)單的隨機(jī)變量 x1， x2， ， xn 的函數(shù) 得到 x1， x2， ， xn 的抽樣后，即可確定 y 的抽樣，這種方法叫作替換法抽樣。即 6. 替換抽樣方法 ),( 21 nxxxgy ),( 21 nf XXXgY 例 11. 散射方位角余弦分布的抽 樣 散射方位角 在 0,2上均勻分布，則其正弦和余弦 sin和 cos服從如下分布： 直接抽樣方法為： 其它當(dāng)011111)( 2xxxf 2co sco s2s ins in 令 =2， 則 在 0, 上均勻分布 ， 作變換 其中 01， 0， 則 (x,y) 表示上半個(gè)單位圓內(nèi)的點(diǎn) 。 如果 (x,y) 在上半個(gè)單位圓內(nèi)均勻分布 ， 則 在 0, 上均勻分布 ， 由于 s inc o syx2222s inc o syxyyxx222222222co ss i n22s i ns i ns i nco s2co sco syxxyyxyx 因此抽樣 sin和 cos的問題就變成在上半個(gè)單位圓內(nèi)均勻抽樣 (x,y) 的問題 。 為獲得上半個(gè)單位圓內(nèi) 的均勻點(diǎn) ， 采用挑選法 ， 在 上半個(gè)單位圓的外切矩形內(nèi) 均勻投點(diǎn) （ 如圖 ） 。 舍棄圓外的點(diǎn) ， 余下的就是所要求的點(diǎn) 。 抽樣方法為： 抽樣效率 E=/40.785 21 yx2221212221222122212s i n,c o s1 為實(shí)現(xiàn)散射方位角余弦分布抽樣 ， 最重要的是在上半個(gè)單位圓內(nèi)產(chǎn)生均勻分布點(diǎn) 。 下面這種方法 ， 首先在單位圓的半個(gè)外切正六邊形內(nèi)產(chǎn)生均勻分布點(diǎn) ，如圖所示 。 于是便有了抽樣效率更高的抽樣方法： 抽樣效率 222121222122212221221121332s i n,33c o s131,123 906.032 E例 12. 正態(tài)分布的抽 樣 標(biāo)準(zhǔn)正態(tài)分布密度函數(shù)為： 引入一個(gè)與標(biāo)準(zhǔn)正態(tài)隨機(jī)變量 X獨(dú)立同分布的隨機(jī)變量 Y，則（ X,Y）的聯(lián)合分布密度為： 作變換 2221)( xexf 2)( 2221),( yxeyxf s inco syx 則 （ ,） 的聯(lián)合分布密度函數(shù)為： 由此可知 ， 與 相互獨(dú)立 ， 其分布密度函數(shù)分別為 分別抽取 ， ： 222),( ef 21)()(2212 ff e212ln2 從而得到一對(duì)服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量 X和 Y： 對(duì)于一般的正態(tài)分布密度函數(shù) N(,2) 的抽樣 ， 其抽樣結(jié)果為： ffffYYXX)2s i n(ln2)2c os (ln22121ffYX例 13. 分布的抽 樣 分布密度函數(shù)的一般形式為： 其中 n， k為整數(shù)。為了實(shí)現(xiàn) 分布的抽樣，將其看作一組簡(jiǎn)單的相互獨(dú)立隨機(jī)變量的函數(shù)，通過這些簡(jiǎn)單隨機(jī)變量的抽樣，實(shí)現(xiàn) 分布的抽樣。設(shè) x1， x2， ，xn 為一組相互獨(dú)立、具有相同分布 F(x) 的隨機(jī)變量，k為 x1， x2， ， xn 按大小順序排列后的第 k個(gè)，記為： 10)1()!()!1(!)( 1 xxxknknxf knk),( 21 nkk xxxR 則 k的分布函數(shù)為： 當(dāng) F(x)=x 時(shí) ， 不難驗(yàn)證 ， k的分布密度函數(shù)為 分布 。 因此 ， 分布的抽樣可用如下方法實(shí)現(xiàn)： 選取 n個(gè)隨機(jī)數(shù) ， 按大小順序排列后取第 k個(gè) ， 即 ininkiin xFxFCxF k )(1)()(ininkiin xxCxF k )1()(),( 21 nkf RX 7. 隨機(jī)抽樣的一般方法 加抽樣方法 減抽樣方法 乘抽樣方法 乘加抽樣方法 乘減抽樣方法 對(duì)稱抽樣方法 積分抽樣方法 1) 加抽樣方法 加抽樣方法是對(duì)如下加分布給出的一種抽樣方法： 其中 Pn0， ，且 fn(x)為與參數(shù) n有關(guān)的分布密度函數(shù)， n=1， 2， 。 由復(fù)合分布抽樣方法可知，加分布的抽樣方法為 :首先抽樣確定 n，然后由 fn(x)中抽樣 x，即： 1)()(nnn xfPxf11nnPn1nn1n1nn PP,當(dāng) nffXX例 14. 多項(xiàng)式分布抽 樣 多項(xiàng)式分布密度函數(shù)的一般形式為： 將 f(x) 改寫成如下形式： 則該分布的抽樣方法為： 0)(iii xaxf00)()1(1)(iiiiii xfPxiiaxf n0ii1n0ii11 PP),m ax (當(dāng) nfX 例 15. 球殼內(nèi)均勻分布抽 樣 設(shè)球殼內(nèi)半徑為 R0，外半徑為 R1，點(diǎn)到球心的距離為 r，則 r的分布密度函數(shù)為 分布函數(shù)為 該分布的直接抽樣方法是 其它當(dāng)03)( 1030312RrRRRrrf3031303)(RRRrrF 31303031 )( RRRr f 為避免開立方根運(yùn)算 ， 作變換： 則 x 0,1， 其分布密度函數(shù)為： 其中 001 )( RxRRr 132)(33)()(200102201 RxRRRxRRxf211020 RRRR 則 x及 r的抽樣方法為： 001432322101201)(),m a x (),m a x (33RXRRrXXXRRRfffff 2) 減抽樣方法 減抽樣方法是對(duì)如下形式的分布密度所給出的一種抽樣方法： 其中 A1、 A2為非負(fù)實(shí)數(shù)， f1(x) 、 f2(x)均為分布密度函數(shù)。 減抽樣方法分為以下兩種形式： )()()( 2211 xfAxfAxf 以上兩種形式的抽樣方法，究竟選擇哪種好，要看 f1(x) 、 f2(x)哪一個(gè)容易抽樣，如相差不多，選用第一種方法抽樣效率高。 （ 1） 將 f (x)表示為 令 m表示 f2(x) f1(x)的下界 ， 使用挑選法 ， 從 f1(x)中抽取 Xf1 抽樣效率為： )()()()(12211 xfxfAAxfxf 111)()(12212211ffffXXXfXfmAAAmAAA211mAAE （ 2） 將 f (x)表示為 使用挑選法 ， 從 f2(x)中抽取 Xf2 抽樣效率為： 22112 )()()()( AxfxfAxfxfmEmAAmE 21 22221221211)()(ffffXXmAAmAXfXfmAAmA例 16. 分布抽樣 分布的一個(gè)特例： 取 A1 2， A2 1， f1(x) 1， f2(x) 2x，此時(shí) m 0，則根據(jù)第一種形式的減抽樣方法，有 或 10),1(2)( xxxf 221 1fX 212 1fX 由于 1 1可用 1代替 ， 該抽樣方法可簡(jiǎn)化為： 對(duì)于 2 1的情況 ， 可取 Xf 1 ， 因此 與 分布的推論相同 。 212fX ),m in ( 21 fX 如下形式的分布稱為乘分布： 其中 H(x)為非負(fù)函數(shù)， f1(x)為任意分布密度函數(shù)。 令 M為 H(x)的上界，乘抽樣方法如下： 抽樣效率為： 3) 乘抽樣方法 )()()( 1 xfxHxf ME111)(fffXXMXH 例 17. 倒數(shù)分布抽樣 倒數(shù)分布密度函數(shù)為： 其直接抽樣方法為： 下面采用乘抽樣方法，考慮如下分布族： 其中 i = 1， 2， ，該分布的直接抽樣方法為： axxaxf 1,1ln1)( af eaX lniif aXi 1)1(1 axxxaixfiii 1,)1(1)( 11 利用這一分布族 ， 將倒數(shù)分布 f(x) 表示成 : 其中 ， 乘法分布的抽樣方法如下 : 該分布的抽樣效率為： )()()( xfxHxf i,1)(,ln)1(,ln)1()(1111iiiixMxHaaiMxaaixH iifiaXa1)1(11)1(21211 )1(ln1 iaiaE例 18. 麥克斯韋 (Maxwell)分布抽樣 麥克斯韋分布密度函數(shù)的一般形式為： 使用乘抽樣方法，令 該分布的直接抽樣方法為： 0,2)(23 xexxf x2ln231fX0,32)( 321 xexf x 此時(shí) 則麥克斯韋分布的抽樣方法為 : 該分布的抽樣效率為： eMxexxHx 2270,3)( 3121 22221ln23lnfXe 7 9 5.0272 eE 在實(shí)際問題中，經(jīng)常會(huì)遇到如下形式的分布： 其中 Hn(x)為非負(fù)函數(shù)， fn(x) 為任意分布密度函數(shù)，n=1， 2， 。不失一般性，只考慮 n=2的情況： 將 f(x) 改寫成如下的加分布形式： 4) 乘加抽樣方法 1)()()(nnn xfxHxf)()()()()( 2211 xfxHxfxHxf )()()()()()()(*22*1122221111xfPxfPxfPxHPxfPxHPxf 其中 )()()()()()()()()()(222*2111*1222111xfPxHxfxfPxHxfdxxfxHPdxxfxHP 乘加抽樣方法為： 該方法的抽樣效率為： 22222)(fffXXMXH 11112)(fffXXMXH 11 P 2221212221111 MPMPMPPMPPE 這種方法需要知道 P1的值 (P2=1 P1） ， 這對(duì)有些分布是很困難的 。 下面的方法可以不用計(jì)算 P1 ： 對(duì)于任意小于 1的正數(shù) P1 ， 令 P2=1 P1 ； ynnPyFyxfyxfyxfyPxHyPxHyxH)(2),(21),()(2)(21,)(),(12122211),m in (),m a x (22112211MPMPEPMPMM則采用復(fù)合挑選抽樣方法，有： 當(dāng)取 時(shí) ， 抽樣效率最高 這時(shí) ， 乘加抽樣方法為： 22222)(fffXXMXH 11112)(fffXXMXH 2111 MMM 2121MME21222111 MMMPMMMP 由于 可知第一種方法比第二種方法的抽樣效率高 。 0)()()(2)()()()()(121212211221212121222121222121212221212221212221212122211212122122212121MMMMPMPMMMMMPPMMPMPMMMMMMMPPMMPMPMMMMMMMPMMMPMMMMMMPMPEE例 19. 光子散射后能量分布的抽樣 令光子散射前后的能量分別為 和 （以 m0c2 為單位， m0為電子靜止質(zhì)量， c 為光速）， ， 則 x 的分布密度函數(shù)為： 該分布即為光子散射能量分布，它是由著名的 KlinNishina 公式確定的。其中 K() 為歸一因子： 211,1111)(1)(322 xxxxxxKxf x22 )21(21421)21l n ()1(21)( K 把光子散射能量分布改寫成如下形式： 在 1, 1+2 上定義如下函數(shù) : 3222)1(111)(1)(xxxxKxf32221221)1()(2)(11)21)(2)(21)(1221)(xxKxHxKxHxfxxf 則有 使用乘加抽樣方法： )()()()()( 2211 xfxHxfxHxf 221122)(278)()21)(4)(21212121MKxHMKxHXXff 光子散射能量分布的抽樣方法為： 該方法的抽樣效率為： 222211132321232)1(427212942711212121fffffffffXXXXXXXXX )294(4)()21(27121 KMME 乘減分布的形式為： 其中 H1(x) 、 H2(x)為非負(fù)函數(shù)， f1(x)、 f2(x) 為任意分布密度函數(shù)。 與減抽樣方法類似，乘減分布的抽樣方法也分為兩種。 5) 乘減抽樣方法 )()()()()( 2211 xfxHxfxHxf （ 1） 將 f (x) 表示為 令 H1(x)的上界為 M1， 的下界為 m， 使用乘抽 樣方法得到如下乘減抽樣方法： )()()()(1)()()(112211 xfxHxfxHxHxfxf 11111)()()()()1(112211ffffffXXXfXfXHXHmM)()()()(1122xfxHxfxH （ 2） 將 f (x) 表示為 令 H2(x)的上界為 M2， 使用乘抽樣方法 ， 得到另一種乘減抽樣方法： 1)()()()()()()(221122 xfxHxfxHxHxfxf 22222)()()()()1(22112ffffffXXXHXfXfXHmMm例 20. 裂變中子譜分布抽樣 裂變中子譜分布的一般形式為： 其中 A， B， C， Emin， Emax 均為與元素有關(guān)的量。令 其中 為歸一因子， 為任意參數(shù)。 m a xm i n,sh)( EEEBEeCEfAE m a xm i n21 ,)()( EEEeEfEfE 相應(yīng)的 H1(E)， H2(E) 為： 于是裂變中子譜分布可以表示成乘減分布形式 : 容易確定 H1(E) 的上界為： 為提高抽樣效率 ， 應(yīng)取 使得 M1 達(dá)到最小 ， 此時(shí) 1141e xp21 ABCMBEEACEHBEEACEH11e x p2)(11e x p2)(21)()()()()( 2211 EfEHEfEHEf 116181ABABA 取 m 0， 令 則裂變中子譜分布的抽樣方法為 : 抽樣效率 4,11 BA1111221m i n)( l nex plnfffffEEBEEEE eCME 211 對(duì)稱分布的一般形式為： 其中 f1(x) 為任意分布密度函數(shù)，滿足偶函數(shù)對(duì)稱條件，H(x) 為任意奇函數(shù)，即對(duì)任意 x滿足： 對(duì)稱分布的抽樣方法如下：取 =2 1 6) 對(duì)稱抽樣方法 )()()( 1 xHxfxf )()()()( 11xHxHxfxf1ff XX 1ff XX )()(111 ffXfXH 證明： 因?yàn)?=2 1， x 相當(dāng)于 ， 因此 xffxfffxfffxfffxfffxfffxffffxfffffffffffffffffdXXfdXXHXfdXXHXfdXXHXfdXXHXfdXXHXfdXXfXfXHdXXfXfXHXfXHxXPXfXHxXPXfXHxXPXfXHxXPxXP11111111111111111111111111111111111)()()()()(21)()(21)()(21)()(21)()()(121)()()(121)()(121,)()(121,)()(121,)()(121,)(1111111111111 x121例 21. 質(zhì)心系各向同性散射角余弦分布抽樣 在質(zhì)心系各向同性散射的假設(shè)下，為得到實(shí)驗(yàn)室系散射角余弦，需首先抽樣確定質(zhì)心條散射角余弦： 再利用下面轉(zhuǎn)換公式 : 得到實(shí)驗(yàn)室系散射角余弦 L。其中 A為碰撞核質(zhì)量，C、 L 分別為質(zhì)心系和實(shí)驗(yàn)室系散射角。 12co s11,21)(CCCCfCCLLAAA211c os2 為避免開方運(yùn)算 ， 可以使用對(duì)稱分布抽樣 。 根據(jù)轉(zhuǎn)換公式可得： 依照質(zhì)心系散射各向同性的假定 ， 可得到實(shí)驗(yàn)室系散射角余弦 L 的分布如下： 該密度函數(shù)中的第一項(xiàng)為偶函數(shù) ， 第二項(xiàng)為奇函數(shù) ，因而是對(duì)稱分布 。 其中 11,212121)(2222 LLLLLAAAf 11,12121)(22221 LLLLAAAf 111 222 LLLC AA 從 f1(L) 的抽樣可使用挑選法 然后再以 的概率決定接受或取負(fù)值 。 上述公式涉及開方運(yùn)算 ， 需要進(jìn)一步簡(jiǎn)化 。 22221212LLLAA1221221221121LAAAA 注意以下事實(shí)：對(duì)于任意 0a1 令 則上述挑選抽樣中的挑選條件簡(jiǎn)化為： 另一方面 ， 在 即 的條件下 ， 2/a 在 1, 1 上均勻分布 ， 故可令 2/a， 則最終決定取正負(fù)值的條件簡(jiǎn)化為： aaPaP )()( 2222 AAAAa 1121 2212212 21222122221)21(1 AAAA222 a 11 2 a12221 AA 于是 ， 得到質(zhì)心系各向同性散射角余弦分布的抽樣方法為： 112221222122221211)21(1 LLAAAAAA 如下形式的分布密度函數(shù) 稱為積分分布密度函數(shù)，其中 f0(x， y) 為任意二維分布密度函數(shù)， H(x)為任意函數(shù)。該分布密度函數(shù)的抽樣方法為： 7) 積分抽樣方法 )(0)(0),(),()(xHxHd xd yyxfdyyxfxf000)(ffffXXXHY 證明：對(duì)于任意 x dxxfd x d yyxfd x d yyxfXHYPXHYdxxXxPXHYdxxXxPdxxXxPxHxHfffffffff)(),(),()()(,)()()(0)(000000000 例 22. 各向同性散射方向的抽樣 為了確定各向同性散射方向 ，根據(jù)公式： 對(duì)于各向同性散射， cos在 1, 1上均勻分布， 在 0, 2上均勻分布。由于 直接抽樣需要計(jì)算三角函數(shù)和開方。 kji wvu co ss ins inco ss inwvu22 1c o s1s i n w 定義兩個(gè)隨機(jī)變量： 可以證明 ， 當(dāng) 時(shí) ， 隨機(jī)變量 x 和 y 服從如下分布 : 定義區(qū)域?yàn)椋?xxAyxxyAyxf12,12m i n0,111221),(22232222212322222123222221AAyAAAAx12322 則 w cos 的分布可以用上述分布表示成積分分布的形式： 令 ， 則屬于上述積分限內(nèi)的 y 一定滿足 條件 。 3 163A12322 21),(),()(31213121)()( xxd x d yyxfdyyxfxf 各向同性散射方向的抽樣方法為： 抽樣效率為： 2322222123222221232222213123222221211223222221c o s2s i ns i n2c o ss i n)(AAAAwAAAvAAAuAA555.012 2 AE 8. 隨機(jī)抽樣的其它方法 偏倚抽樣方法 近似抽樣方法 近似 -修正抽樣方法 多維分布抽樣方法 指數(shù)分布的抽樣 使用蒙特卡羅方法計(jì)算積分 時(shí)，可考慮將積分 I改寫為 其中 f *(x) 為一個(gè)與 f (x) 有相同定義域的新的分布密度函數(shù)。于是可以這樣計(jì)算積分 I： 這里 Xi 是從 f *(x) 中抽取的第 i 個(gè)子樣。 1) 偏移抽樣方法 dxxfxgdxxfxf xfxgI )()()()( )()( * dxxfxgI )()(NiiiiNi iiNiiN XgXWNXgXfXfNXgNI11*1* )()(1)()()(1)(1 由此可以看出 ， 原來(lái)由 f (x) 抽樣 ， 現(xiàn)改為由另一個(gè)分布密度函數(shù) f *(x) 抽樣 ， 并附帶一個(gè)權(quán)重糾偏因子 這種方法稱為偏倚抽樣方法 。 從 f (x) 中抽取的 Xf ， 滿足 而對(duì)于偏倚抽樣 ， 有 一般情況下 ， Xf 是具有分布 f (x) 總體的簡(jiǎn)單子樣的個(gè)體 ， 只代表一個(gè) 。 Xf* 是具有分布 f *(x) 總體的簡(jiǎn)單子樣的個(gè)體 ， 但不代表一個(gè) ， 而是代表 W(Xf*) 個(gè) ，這時(shí) Xf*是帶權(quán) W(Xf*)服從分布 f (x) 。 )()()( * xfxfxW dxxfdxdXxP f )()( dxxfdxxfxWdxdXxPXW ff )()()()()( * 在實(shí)際問題中，分布密度函數(shù)的形式有時(shí)是非常復(fù)雜的，有些甚至不能用解析形式給出，只能用數(shù)據(jù)或曲線形式給出。如中子散射角余弦分布多數(shù)是以曲線形式給出的。對(duì)于這樣的分布，需要用近似分布密度函數(shù)代替原來(lái)的分布密度函數(shù)， 用近似分布密度函數(shù)的抽樣代替原分布密度函數(shù)