計(jì)算物理第二章

1、第二章蒙特卡洛方法1§2.1MC方法的理論基礎(chǔ)概述 蒙特卡羅方法又稱隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法。 蒙特卡羅方法是數(shù)值計(jì)算方法之一。但是它與一般數(shù)值計(jì)算方法有很大區(qū)別。它以概率統(tǒng)計(jì)理論為基礎(chǔ)。由于蒙特卡羅方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過程，解決一些一般數(shù)值方法難以解決的問 題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。 基本思想：當(dāng)所求問題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的期望值，或與概率、數(shù)學(xué)期望有關(guān)的量時(shí)，通過某種試驗(yàn)的方法，得 出該事件發(fā)生的頻率，或該隨機(jī)變量若干個(gè)觀察值的算術(shù)平均值，根據(jù) 大數(shù)法則得到問題的解；由中心極限定理可以得到該解的概率誤差。2概率理論的基礎(chǔ)知識(shí)在

2、相同的實(shí)驗(yàn)條件下，隨機(jī)事件A，B按各自確定的概率發(fā)生-或概率 A.OR.B：P(A+B) = P(A) + P(B) P(AB)-與概率 A.AND.B：P(A*B) = P(A B)*P(B) = P(B A)*P(A)-條件概率 P(A|B) = 在隨機(jī)事件B發(fā)生的條件下，A發(fā)生的概率-互斥 P(A*B) = 0，i.e. 隨機(jī)事件AB不能在同一實(shí)驗(yàn)中同時(shí)發(fā)生-相互獨(dú)立 P(A*B) = P(A)*P(B)3古典概率：隨機(jī)事件A構(gòu)成互斥完備集合Ai，則任意隨機(jī)事件B可表述為P(B)=åP(B|Ai )P(Ai )i= P(Ai B) =P(B|Ai )P(Ai )P AB

3、9;P(B|A )P(A )iP(B)jjj4貝葉斯Bayes公式：全概率公式：）在相同的確定實(shí)件下，對(duì)的觀測(cè)無法給出單一固定值;）必須依據(jù)遍舉測(cè)量原則，對(duì)所有可能取值給出發(fā)生概率離散變量舉例:3MeV光子入射屏蔽鉛板的全吸收反應(yīng)過程反應(yīng)類型X :光電效應(yīng)Compton散射電子對(duì)產(chǎn)生23反應(yīng)幾率:其中e1 + e2 + e3 = 100%i.e.nX = xi , å pi = 1i =15隨機(jī)變量:連續(xù)型隨機(jī)變量:在連續(xù)區(qū)間取值，其取某確定值概率由分布密度函數(shù)給出f ( x) dx = Px < X < x + dx分布函數(shù)xF ()f ( x)dxò-

4、65;則有+¥F (+¥) = ò-¥f ( x)dx º 16聯(lián)合分布密度:描述兩個(gè)(i.e.多維)隨機(jī)變量與的相互關(guān)聯(lián)+¥òx=f1f ( x, y )dy-¥相互獨(dú)立:f ( x, y ) = f1 (x) × f2 ( y)7隨機(jī)變量的特征值1)期望值(mean):出現(xiàn)幾率最大或概率中心的觀測(cè)值+¥E y(x)y = ò y(x) × f (x)dx-¥2)標(biāo)準(zhǔn)方差(standard deviation):隨機(jī)變量x分布對(duì)期望值的離散程度+¥=(y

5、 - E y(x)2 × f (x)dx = E(y - y)2òDy(x)s2y-¥) 特征運(yùn)算:Ec = c, D= 0, Ecx = cEx, Dcx = c2DxEx1 + x2 = Ex1 + Ex2 Dx1 + x2 = Dx1 + Dx2 + 2 × Covx1 , x2 8幾種著名分布1) 二項(xiàng)式分布(Binominal): 現(xiàn)k次的幾率發(fā)生幾率為p, 不發(fā)生為q=(1-p),則N次試驗(yàn)中出N !(1- p)N -kPk =pkk !( N - k )!其中k=0,1,2,3, 0p1, p+q1Ek = Np, Dk = Np(1- p

6、)92) 泊松分布(Poisson): 在相同驗(yàn)條件下，相同時(shí)間內(nèi)，隨機(jī)過程發(fā)生k次的幾率lk- lP k =ek ³ 0, l > 0k !E k = Dk = l其中關(guān)于分布參數(shù)k有當(dāng)à時(shí)，Poisson分布過渡到Gaussian分布103) 高斯分布(Gaussian):- 1 ( x- m )21-¥< x < +¥E(x) = mD(x) = s 2N ( X ; m,sf (x) =s 2e22ps 標(biāo)準(zhǔn)化xà(x-m)/s, 則正態(tài)分布- 1 x221f (x) =eN (0,1)2p11分布函數(shù)yF(y) =

7、ò-1 x21edx22p-¥對(duì)稱分布F(- y) = 1- F( y)124)指數(shù)分布(Exponential):描述自由粒子壽命，或粒子平均自由程分函1x > 0, l > 0e - x / lf ( x ) =分布函數(shù)zò= 1 - e - z / lF ( z ) =f ( x ) dx-¥ 5)均勻分布(uniform):1f ( x ) =, a £x £ bb - aa + b( b - a ) 2E ( x ) =, D ( x ) =其中21 21313大數(shù)法則:在函數(shù)f(x)定義域a,b內(nèi)，以均勻概率分

8、布隨機(jī)地取N個(gè)數(shù)xi，函數(shù)值之和的算術(shù)平均收斂于函數(shù)的期望值1N1åb()f ( x)dxòlim IN ºlimfxIib - aNN ®¥N ®¥ai=1這一法則指出, 不論分布函數(shù)f(x)的分布如何, 只要N足夠大, 則算術(shù)平均與數(shù)學(xué)期望值可無限接近, 也就是說, 算術(shù)平均以幾率收斂于其數(shù)學(xué)期望值。14隨機(jī)變量大量抽樣值累加而成的隨機(jī)變量服從高斯分布。高斯分布可以由給定的期望值和方差完全確定下來，1N (m ,s2 ) =exp é- ( x - m )22 ù/ 2sëûs2p

9、高斯分布隨機(jī)測(cè)量表示x = m ±置信水平 (1-a)×sìs ( f )lim Pr ob-l£ I( )s f ü 1lòl2e-t /2dt = 1-a- I <=íýþ(1)n2p-ln®¥nnîsf< lIn - In1515中心極限定理:a,1-a和的數(shù)值關(guān)系由表可見,當(dāng)=2.5758時(shí)1立的幾率已經(jīng)為99%,也就是說,該式的可靠性已相當(dāng)高.16a (顯著水平)0.50.050.020.011-a (置信水平)0.50.950.980.990.6745

10、1.96002.32632.5758y = y(xi )統(tǒng)計(jì)量:隨機(jī)變量X的一組測(cè)量樣本xi的函數(shù)例，樣本平均值Nx = 1 åyxiNi=1作為維相互獨(dú)立量)隨機(jī)變量的函數(shù)，y亦為隨機(jī)變量，亦存在分布E(x ) = E(x) = m收斂于期望值(大數(shù)定理)s 2D(x)D(x ) =期望值測(cè)量誤差(中心極限)NN17蒙特卡羅方法基礎(chǔ)知識(shí)針對(duì)待求問題，根據(jù)物理現(xiàn)象本身的統(tǒng)計(jì)規(guī)律，或人為構(gòu)造一合適的依賴隨機(jī)變量的概率模型，使某些隨機(jī)變量的統(tǒng)計(jì)量為待求問題的解， 進(jìn)行大統(tǒng)計(jì)量Nà的統(tǒng)計(jì)實(shí)驗(yàn)方法或計(jì)算機(jī)隨機(jī)模擬方法。理論依據(jù):- 均勻分布的算術(shù)平均收斂于真值 (大數(shù)法則)- 置信

11、水平下的統(tǒng)計(jì)誤差 (中心極限)待求問題：1）自然界中本身存在的隨機(jī)過程，如粒子衰變過程、粒子在介質(zhì)中的輸運(yùn)過程等; 2）以慨率模型來解決不直接具有隨機(jī)性的確定性問題，如求, 求積分。18例1.Buffon投針實(shí)驗(yàn)(1777年)1）平行線間距=針長=s；針與平行線垂線方向夾角 則相交概率為P = s cos a= cos as2）各項(xiàng)同性均勻拋針，i.e.夾角在0, 均勻分布p= ò cos a × f2a da = pEcos a03）設(shè)N拋，M相，相概率的望值E( cos a ) = 2 º Mp = 2N（Nà大數(shù)定理）pNM19問: 的 量?p =

12、 M= 2服從二項(xiàng)式分布，單次相交概率pNpM (1- p)N -MCMNE(M ) = NpD(M ) = Np(1- p)s 2 (M )則p(1- p)s 2 ( p) = s 2 ( M ) =1s 2 (M ) =N 2NN202Qp =pdpdpp 2 (p - 2)p(1- p)4s 2 (p ) = () ×s( p) =22p4N2Nsp= 2.37Np = 2N ± 2.37NM對(duì)真值的測(cè)量精度與測(cè)量次數(shù)平方根反比，即106次實(shí)驗(yàn)才精確到10-321的投針實(shí)驗(yàn)結(jié)果實(shí)驗(yàn)者年份投針次數(shù)的實(shí)驗(yàn)值沃爾弗(Wolf)185050003.1596史密斯(Smith)

13、185532043.1553?？怂?Fox)189411203.1419190134083.1415929拉查里尼(Lazzarini)221I = òf (x)dx例2.投點(diǎn)法求定積分0隨機(jī)地向x0,1, y0,1正方形內(nèi)投點(diǎn)N， 統(tǒng)計(jì)落在曲線下的點(diǎn)數(shù)M，當(dāng)總擲點(diǎn)數(shù)Nà時(shí)N-MM=INM建立恰當(dāng)?shù)母怕誓Ｐ?，即確定某個(gè)隨機(jī)事件A或隨機(jī)變量X，使得待求問題的解等 于隨機(jī)事件出現(xiàn)的概率或隨機(jī)變量的數(shù)學(xué)期望值。然后重復(fù)進(jìn)行多次的隨機(jī)實(shí)驗(yàn)，對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)平均，求出A出現(xiàn)的頻數(shù)或X的平均值作為問題的近似解。這種方法也叫做蒙特卡洛模擬。23蒙特卡羅方法的優(yōu)缺點(diǎn)缺點(diǎn)優(yōu)點(diǎn)Ø

14、16;收斂速度慢。誤差具有概率性。在粒子輸運(yùn)問題中，計(jì)算結(jié)果能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過程。受幾何條件限制小。與系統(tǒng)大小有關(guān)。收斂速度與問題的維數(shù)無關(guān)。誤差容易確定。程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)。具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力。24§.隨機(jī)數(shù)與偽隨機(jī)數(shù)進(jìn)行計(jì)算機(jī)模擬需要大樣本的均勻分布隨機(jī)數(shù)數(shù)列，如何獲得？- 真隨機(jī)數(shù)：由隨機(jī)物理過程來產(chǎn)生，例如：放射性衰變、電子設(shè)備的熱噪音、宇宙射線的觸發(fā)時(shí)間等等- 偽隨機(jī)數(shù)：由計(jì)算機(jī)按遞推公式大量產(chǎn)生馮.諾曼平方取中法()= éë2-r xù22rxmod 2ûn+1nx= x/

15、 22rnnxn Î 0,1加同= (axn + c)(mod m)xn = xn / mxn +125偽隨機(jī)數(shù)的統(tǒng)計(jì)檢驗(yàn) 均勻性檢驗(yàn):-0,1分成k個(gè)相等子區(qū)間，進(jìn)行N次抽樣，投入各子區(qū)間-如均勻,則各區(qū)間落入數(shù)Ni應(yīng)為N = å Niki=1Ni = m = s2= N k計(jì)數(shù)PoissonàGaussian-Ni可視為(,)的一組無關(guān)樣本測(cè)量，服從26獨(dú)立性檢驗(yàn): 即xi與xi+1的前后無關(guān)性-0,1上進(jìn)行2N次抽樣，分成兩個(gè)序列X : x1,x3 ,.x2i-1,.,x2 N -1Y : x2 ,x4 ,.x2i ,.,x2 N在XY平面內(nèi)劃分k×

16、;k方格，如獨(dú)立, 則各格內(nèi)落入數(shù)應(yīng)為-kN = å niji, j =1Nm = s 2=n2ijk則服從c2分布-滿足以上統(tǒng)計(jì)性檢驗(yàn)的遞推抽樣序列，可視為0,1均勻分布偽隨機(jī)數(shù)2727§2.3任意分布的偽隨機(jī)變量的抽樣A. 離散型分布X :x1 , x2 , x3 ,., xN 隨機(jī)變量“本征值”f : p1 , p2 , p3 ,., pN “本征值”概率密度歸一化，幾率守恒分布函數(shù)0 £ F ( x) = å pN£1åp = 1iixi £xi=1則取0,1均勻分布隨機(jī)數(shù)x，按分布函數(shù)不等式抽樣x(j -1 ) &

17、#163; x <F()jh = x jFx28 離散變量舉例:3MeV光子入射屏蔽鉛板的全吸收反應(yīng)過程反應(yīng)類型X 反應(yīng)截面反應(yīng)幾率:：:光電效應(yīng)s1e1= s1/sCompton散射s2e2= s2/s電子對(duì)產(chǎn)生s3e3= s3/s其中，s1+s2+s3=s，e1 + e2 + e3 = 100%x1NNx1 < e1?x12?電子對(duì)產(chǎn)生YYCompton光電效應(yīng)291. 反函數(shù)直接抽樣法：類比于離散型隨機(jī)抽樣，求分布函數(shù)x()f ( x)dx £ 1ò0 £ F=x-¥則令x = F (h)h £ x的概率為：h = F -1

18、(x )該子樣中= p= pF (x )F(x ) £x £ F (x)=0 × dx +1× dx = F (x)0òòh-£ x1px-¥030B. 連續(xù)型分布例對(duì)指數(shù)分布的直接抽樣ìle -lx , x > 0, l > 0( )x = íf0, 其它.î積分得到分布函數(shù)F (x) =f (t )xxòòl-l-ldt =dt = 1 -txee-¥0令x = F (h) = 1 - e -lh則指數(shù)分布的隨機(jī)變量抽樣為h = - 1 l

19、n (1- x )- 1 ln xll312. 變換抽樣法：隨機(jī)變量x按f(x)抽樣已知，其函數(shù)y=y(x)的概率密度有f (x)dx = g( y)dydxg( y) = f (x-1( y) ×dy二維情況：聯(lián)合分布密度函數(shù)g(u,v)的隨機(jī)變量(h,d)抽樣方法已知，則其函數(shù)x = g1u( ,)vy= g 2 ( , vu)¶u¶u抽樣可依據(jù)¶x¶v¶y¶vf (x, y) = g(u(x, y), v(x, y) ׶x¶y32例Gaussian正態(tài)分布(x - m )2ì&

20、#252;ýþ 11f ( )x =s expí-2s2p2î可由標(biāo)準(zhǔn)正態(tài)分布代替f (x) =exp- x 2 / 21h = sd + m2p解： 1）相互獨(dú)立均勻抽樣 , ,2) 構(gòu)造二維隨機(jī)變量. .,=×1üüu = exp ì- 1+ y2x2íýï-2 ln u cos(2p v)-2 ln u sin (2p v)ìïx =î12p2þýïïþíïî y =ta

21、n-1vy / x33則，x,y 聯(lián)合分布密度函數(shù)為f ( x, y ) = g (u (代入Jacobi行列式，有, y ), v ( x, y )× Jìü12p12()f ( x, y ) = f ( x)× f ( y )exp-x +22yíýîþ其中，exp- x 2 / 2exp- y 2 / 211f (x) =f (y) =2p2p34改進(jìn)的Maraglia方法:u,v0,1(1)產(chǎn)生0,1區(qū)間上的獨(dú)立均勻分布隨機(jī)數(shù)u,vw = (2u - 1)2 + (2v - 1)2(2)計(jì)算w=(2u-1)

22、2+(2v-1)2如果 w>1，回到步驟(1)；否則，執(zhí)行(4)(3)w >1?z = - 2 ln(w)/ w1 / 2(4)計(jì)算Y舍選效率79%(5)得到二維相互獨(dú)立標(biāo)準(zhǔn)正態(tài)分布抽樣x = uz, y = vzz=-2lnw/w1/2X=uz35N3. 舍選抽樣法：（是第二類舍選法的特殊情況。如抽樣概率密度函數(shù)g(x,y)為均勻分布，則為第一類舍選法)( x) = 1L =maxl定義域xÎ a ,b x1,x2d =a+(b-a)x1x 2 £ lf (d )Yhx = dN36第一類舍選例對(duì)隨機(jī)變量抽樣1l0££.f( x) 2 &#

23、236;x,0x1,= maxf(=)x=2L=íx Î a , b,其它î( x) =xx =2，x直接抽樣反函數(shù)法F解：一舍改進(jìn)x1,xh= max(x ,x )2x12廣 推ìnxn-1 , x Î0,1, n =21, .f (x) = íNx£ x0,其它.î21Y( 1n,.).,mx =ax,2h= xx137（是第三類舍選法的特殊情況)判據(jù)函數(shù)h ( x )òf ( x) L抽樣效率g ( x, y)dy-¥抽樣密度函數(shù)(y,) = g1 ( ) x2( )g，且yy為均勻分布如

24、x,y相互g獨(dú)立x，( y ) = ì1,y Î0,1在y的定義域上0 £ h (x ) £ 1gí20, 其它î則有h ( x )f ( x ) = Lg1 ( x ) ò g 2 ( y)dy = L × h( x) × g1 ( x)-¥第二類舍選法的形式效率判據(jù)抽樣38第二類舍選則對(duì)任意密度函數(shù)構(gòu)造 f ( x)h(x) × f ( x )( x ) = L× g ( x ) º L × h ( x ) × g ( x ) ,fLg (

25、x)( x )L效率判據(jù)抽樣x, 按g(x)抽取hgY均勻分布g2(y)=1 à 0h(x)1，即f (x) > 1Lmaxh(x)x定義域Ny=xh(hg)?密度概率函數(shù)g(x)，易于抽樣Y抽樣效率1/Lhx=hg39例采用第二類舍選抽樣標(biāo)準(zhǔn)正態(tài)分布ìüx 21f ( x ) =exp í-(-¥ < x < +¥)ý2p2îþ解：不存在反函數(shù)，湊第二類舍選函數(shù)形式hg= - ln x1g(x) º e-x指數(shù)直接抽樣：h(x) º exp-( x -1)2N2判據(jù)

26、：(h-1)2 £ -2 ln xg22epL º效率：Yhx = hg(0 < x < +¥)偶函數(shù)：40不易直接抽取隨機(jī)變量x分布f(x), 而與隨機(jī)變量y聯(lián)合分布密度g(x,y)存在+¥ò g ( x, y)dy-¥h ( x )ò-¥f ( x)Lg ( x, y)dy聯(lián)合分布第三類舍選法的形式其中 函數(shù)h(x)取值在y定義域上 常數(shù)L的定義為ù-1é+¥ h( x)L = ê ò ò g(x, y)dxdyú1ë-

27、¥ -¥û41數(shù)學(xué)上構(gòu)造第三類舍選第三類舍選抽樣步驟： 由聯(lián)合分布密度函數(shù)g(x,y)抽取隨機(jī)向量 (h x ,h y )按g(x,y)抽取(hx,hy)£ h(h x ) 判別h y是否成立; 若不成立,返回 取分布密度函數(shù)f(x)的抽樣值h = h xNhyh(hx)?抽樣效率1/LYh=hx42Bayes公式則有phx £ x,hy £ h (hx ) phy £ h (hx )£ x hy £ h hx=pxh (t1 )ò-¥xò dt1-¥g (t1 ,

28、 t2 )dt2éùh (t1 )ò-xò ê=ê Lg (t1 , t2 )dt2 ú dt1+¥òh (t1 )òúëûdtg (t , t )dt1122-¥-¥p h £ x43子樣hx概率等于在hyh(hx)條件下，hxx出現(xiàn)的概率例：各向同性方位角余弦的抽樣物理含義f 0,2p各向同性 à h=cos(2px)解：密度函數(shù)ì 11< 1xf (x) = ïpí1- x2ï

29、;î0其它xF (x) = ò-¥xcos-1 xp11需要計(jì)算cos函數(shù)f (t)dt = p òdt =1-1- t 2-1直接抽樣- cosh-1x = 1h = cos(px )p44解獨(dú)立抽取均勻分布x1 與x2，定義1 y(1+ x)2x- x22x=x = 12 1x+ x2222211x2=y(1- x)2y = x+ x21(x1,x2)到(x,y)變換抽樣g(x, y) = 11(x ,x) ×= 1×JJ1241- x2則構(gòu)造(x,y)聯(lián)合密度函數(shù)ì1x < 1, 0 < y < 2,

30、(x,) = ï 41- x2íïî0,其它45類比構(gòu)造h(x)與Lh( x)11L ò g(x, y)dyf (x) =1- x2-¥h( x)h( x)11= L ò0ò0d= Ld1- x4 1- x224則有L = 4 ,h(x) = 1p46抽樣：改進(jìn)：x 1 , x2x 1 , x 2x- x22x = 12 x = x1y = 2x 2 - 1x+ x2222N1Ny =+212x2 + y2 £1Yy £ h(x) = 1Yh = x效率1/L=p/4- y22h = x+ y2

31、x22x1x 2h¢ = 2xyh¢ =21+22x2+ y2474. 復(fù)合抽樣法：復(fù)合分布密度函數(shù)+ ¥òf ( x ) =g ( x | y ) h ( y ) d y-¥ 根據(jù)密度函數(shù)h(y)抽樣 à yh 根據(jù)條件密度函數(shù)g (x | yh ) 抽樣 àxgf= xg ( x| y)h48復(fù)合概率密度函數(shù)f ( x ) = ånhn ( x ) × pn其中+¥ò hn (x) = 10 < pn < 1,pn = 1;n-¥抽樣步驟： 取0,1區(qū)間上均勻

32、分布隨機(jī)數(shù)x , 解下不等式求nn-1nå pi£ å pi<i=1i=1=找到對(duì)應(yīng)的hn (x)，對(duì)其抽樣得到hn49加分布抽樣例：球殼均勻分布的抽樣23 rR0 £ r £ R1f ( r ) =,R-330R1 積分求分布函數(shù)，直接抽樣h = R ()1 / 3xR- R+333100 變量代換r = (R - R )x + R ,l =+ R R+ R 22R1001100密度函數(shù)(R - R3R (R - R ) 23R 2f (x) =+ 0102x + 0 ×1 103x 2lll5051復(fù)合概率密度函數(shù)f (x)

33、 = A1g1 (x) - A2 g2 (x),A1, A2 > 0取g 2 ( x ) ,m =m ingx ,x >0概率密度函數(shù)21g ( x )xÎ a ,b 1則ég ( x) ù0 < f ( x) = g1( x) ê A1- A2 2ú £ g1 ( x)( A1 - A2 m)g1 ( x) ûë即f (x)A1A2g2 (x) £ 10 < h (x)=-1( A - A m)g (x)A - A mA - A mg (x)1211212152減分布抽樣亦即第二

34、類舍選f (x) = ( A1- A2m) × h1(x) × g1(x)效率“L”>1?(0,1，判據(jù)概率密度函數(shù)，抽樣另一種第二類舍選- A2 m × h ( x) × gA1f ( x) =( x)22m其中A1m g1 (x) -A2m0 < h (x) =£ 12A - A m g (x)A - A m1221253f (x) = å Hn (x)gn (x)n考慮n=2情況f ( x) = H1 ( x) × g1 ( x) + H 2 ( x) × g 2 ( x)(2.3.69)令bb&

35、#242;ò+ p2 = 1p1 =p2 =p1H1(x)g1(x)dx,H2 (x)g2 (x)dx,aa則加分布H (x)H2 (x) gp g¢(x) + p g¢ (x)f (x) = pg (x) + p(x) =111221122pp12問題：需要計(jì)算p1與p2積分54乘加分布抽樣另一種方式ì)ü(x)M1H (x)M 2H2 (x) gf (x) = (M+ Mg (x) +1í Mýþ1212+ MM + MMMî121122Mi = maxHi (x)其中xÎa,b令MM=p2

36、 =p1 1, 2M1 + M 2M1 + M 2取H (x)L = L= M+ M0 < h (x) =£ 1 i1212iMi則有f (x) = p1 L1 × h1 (x) × g1 (x) + p2 L2 × h2 (x) × g2 (x)加分布+第二類舍選55f (x) = H1 (x) × g1 (x) - H2 (x) × g2 (x),x Îa,b令H2 (x)g2 (x)m =M = max H (x)min11H (x)g (x)xÎ a,bxÎa,b11有é

37、H (x)g (x) ù0 < f (x) = H1 (x)g1 (x)ê1 - 22ú £ H1 (x)g1 (x)(1 - m) £ M1 (1 - m)g1 (x)H1 (x)g1 (x) ûë再令é Hù1H ( x) g ( x)0 < h ( x) =( x) -£ 122M (1 - m) êú11g ( x)ëû11給出f (x) = M1 (1- m) × h1 (x) × g1 (x)第二類舍選56乘減分

38、布抽樣4. 特殊抽樣方法：- 問題：根據(jù)分布函數(shù)F ( x)直接抽樣最為便捷，但其反函數(shù) F -1 ( y) 難以解析求出- 解答：分布函數(shù) F ( x)F -1( y) » Q( y)性質(zhì) à 近似替代lim F -1 ( y) = -¥lim F -1 ( y) » +¥y Î 0,1y ® 0y ®1最小二乘擬合F -1 ( y) » Q( y) = a + by + cy 2 + a (1 - y)2 ln y + by 2 ln(1 - y)yà0,yà1,57反函數(shù)近似例：

39、標(biāo)準(zhǔn)正態(tài)分布N(0,1)- x22- x2211F -1 ( y )Q( y )y = òdx » òeedx-¥-¥22= 1,2 . 199 )將0,1區(qū)間200等份 (k- x221Q ( y )yk= òedx,2p-¥數(shù)值計(jì)算出對(duì)應(yīng)的值，然后利用逐步回歸法計(jì)算出：a = -0.8268,b = 1.6736,c = 0,a =b =.3315,.331558期望快速正態(tài)分布隨機(jī)變量抽樣；利用均勻分布偽隨機(jī)數(shù)x Î0,1極限近似法21E (x ) =,D(x ) =212x1+x 2+ . + xn則定義

40、Rn=，則+¥+¥n2nE ( Rn ) = ò nx × f ( x)dx = ò nx ×1dxD ( Rn ) =12=,2-¥-¥據(jù)中心極限定理，有D(Rn )nnnnN (Rn ;,) : E(Rn ) = E(Rn ) =,D(Rn ) =2 122n1259引入服從標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量- nR1nx2-ò2e-t / 2dt = N (0,1)lim p(d£ x) =dN (d ; 0,1)=nnn2pn12¥®n¥通常取n=12，則d12= R12

41、 - 6< 6注意：以上n=12處理僅能描述x定義域范圍內(nèi)的正態(tài)分布抽樣605. 多維隨機(jī)向量的抽樣方法:f (x1 , x2 ,× × ×, xn )，定義域?yàn)槎嗑S相互關(guān)聯(lián)隨機(jī)變量的聯(lián)合分布密函數(shù)a1 £ x1 £ b1 , a2 £ x2如存在有限上限b2 ,× × ×, an £ xn £ bn L = sup f (x1 , x2 ,× × ×, xn ) < +¥則可采取第一類舍選< 1(b)- a )x + a ,

42、(b- a )x+ a ,×××, (b)xx- a+ an+111112222nnnnL抽樣效率1E =nLÕ(bi - ai )i=161舍選法三維隨機(jī)向量的聯(lián)合分布密度函數(shù)可寫為條件密度法f (x1 , x2 , x3 ) = f1 (x1 ) × f 2 (x2 | x1 ) × f3 (x3| x1 , x2 )其中Bayes公式，有+¥ +¥üf1 (x1 ) = ò òf (x1 , x2 , x3 )dx2 dx3ïï-¥ -¥+

43、¥f2 (x2 | x1 ) = ò-¥ïf (x1 , x2 , x3 )dx3 / fx )ýï(x | x )ï(|)f)ïïþ2162抽樣步驟h1 = x1f1 (x1 )(1) 由為分布密度函數(shù)產(chǎn)生抽樣值在h1 =x1條件下，由分布密度函數(shù) f2 (x2 |h1) 產(chǎn)生抽樣值 h2= x2(2)h3 = x3在 h1 =x1,h2 =x2條件下，由分布密度函數(shù) f3 (x3 |1 ,h2 )產(chǎn)生抽樣值(3)可進(jìn)一步推廣至n 維隨機(jī)向量f (x1 , x2 ,× ×

44、×, xn ) = f1 (x1 ) × f2 (x2 | x1 ) × f3 (x3 | x1 , x2 ) × × × fn (xn | x1 , x2 ,× × ×xn-1 )( x , x , x ,×××, x)T得到隨機(jī)變量的一組抽樣向量123nh = (h ,h ,h , ×××,h)T123n63例 中子入射角服從聯(lián)合分布密度函數(shù)ìf (j ,q ) = 1 (1 +3 sinj sinq )sinj sin2 q ,&

45、#239;a³ 0,0 £íïî/ 2 >£/ 2,= (3 + 23)/12.x = cosj, y = cosq求余弦抽樣解：變量代換(j,q)à(x,聯(lián)合密度函數(shù)為121- y2 é1+3(1- x2 )(1- y2 ) ù = f (x) f ( y | x)f (x, y) =ëû12(3 + 23)則有æ pö+¥1223òf (x) =f (x, y)dy =+1- x2ç÷1(3 + 23)p43

46、2;ø-¥1 - y 2 (1 +(1 - x 2 )(1 - y 2 ) )f ( y | x) =132+ 231 - x 24364先對(duì) f1 (x) 抽樣，加分布x=hìïí× 4f (x) = p + p1 - x 2112p=3, 3 + 2 3= 2 33 + 2 3ïppïî12Nox < 3 /(3 + 2 3)? 1Yesx2 + x< 1?223NoYes65h = x2產(chǎn)生隨機(jī)數(shù) x2 ,x3產(chǎn)生隨機(jī)數(shù)x1f2 ( y | x = h)加分布抽樣y=d再按p4pf (

47、y | x = h) =4×1- y22p231-h 2+No423<p / 4?x4p2 331-h2+41- x22× 3233+(1- y )2p + 231-h 2x 2 + x < 1?4356NoNox5 + x < 1?26YesYes66d = x5產(chǎn)生隨機(jī)數(shù) x5,x6產(chǎn)生隨機(jī)數(shù) x5,x6產(chǎn)生隨機(jī)數(shù)x4相互獨(dú)立111exp - y / 2×exp - y / 2×××exp - y / 2f ( y , y ,×××, y ) =22212n12n2p2p2p可進(jìn)行n

48、次一維變量正態(tài)分布抽樣 hy ，獲得抽樣樣本向量i= (hy1 ,hy2 ,hy3 ,×××,hynhryT67高維正態(tài)分布隨機(jī)向量抽樣相互關(guān)聯(lián)ìrr ü- 112- nrrf (x1 , x2 ,× × ×xn ) = (2p )× M(c - m ) M(c - m )ýT× expí-12îþ其中cr = ( x , x ,×××, x)mr = E hr = (m , m ,×××, m)

49、T12n12n正定對(duì)稱 n 階協(xié)方差矩陣M的矩陣元= E(-i )()= sjs ij-ijji存在非奇異的下三角矩陣AM = AAT則關(guān)聯(lián)抽樣可由無關(guān)抽樣表示標(biāo)準(zhǔn)正態(tài)分布抽樣hx = m + Ahry68例二維正態(tài)分布的抽樣r()Th = h,h解設(shè)服從二維正態(tài)分布12æs 11s 12ör()Tm =m , mM = çs÷,s12è2122 ø其協(xié)方差矩陣A矩陣為æös 110ç÷A = ç÷ss s- s2ç12 112212 ÷s 11s 11èø抽樣結(jié)果ìhx= m1s 11h y+ï111í=+ (2 )1 /