蒙特卡羅方法1課件

計(jì)算物理3/lesson/ComputationalPhysics蒙特卡羅方法計(jì)算物理3/less1蒙特卡羅方法蒲豐投針收斂性、誤差和優(yōu)缺點(diǎn)任意分布的隨機(jī)數(shù)粒子輸運(yùn)問(wèn)題隨機(jī)過(guò)程模擬梅氏抽樣√蒙特卡羅方法蒲豐投針√2蒲豐投針(1/5)蒙特卡羅方法又稱隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法以概率統(tǒng)計(jì)理論為基礎(chǔ)的能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程解決一般數(shù)值方法難以解決的問(wèn)題隨著電子計(jì)算機(jī)的發(fā)展而發(fā)展首先在核武器的試驗(yàn)與研制中得到了應(yīng)用蒲豐投針?lè)▏?guó)數(shù)學(xué)家蒲豐的1777年出版的著作：“在平面上畫有一組間距為d的平行線，將一根長(zhǎng)度為

l

(l<d)

的針任意擲在這個(gè)平面上，求此針與平行線中任一條相交的概率?！薄唐沿S投針(1/5)蒙特卡羅方法蒲豐投針√3蒲豐投針(2/5)步驟在桌面上畫出間距為

2d的平行線準(zhǔn)備長(zhǎng)度為

2l

(l<d)

的針向桌面隨機(jī)投針如果針與平行線相交，則計(jì)數(shù)器n加

1計(jì)算：計(jì)數(shù)器n與總投針數(shù)N的比例(視作相交概率P)概率分析

P

=

?各條平行線地位等同，僅考慮某條平行線附近的情況平行線方向的x坐標(biāo)對(duì)概率沒(méi)影響針的中點(diǎn)的

y坐標(biāo)在線之間等概率落入(均勻分布在

[0,d])，僅當(dāng)

yl才可能針-線相交針-線的夾角

q

均勻分布在

[0,

p]，q

與

y獨(dú)立xy√蒲豐投針(2/5)步驟概率分析P=?xy√4蒲豐投針(3/5)xy概率

P

=

2l/

(pd)，可計(jì)算圓周率實(shí)驗(yàn)者時(shí)間針長(zhǎng)總投數(shù)相交數(shù)p

值Wolf18500.85,0002,5323.1596Smith18550.63,2041,2183.1554DeMorgan,C18601.06003823.137Fox18840.751,0304893.1595Lezzerini19010.833,4081,8083.1415929Reina19250.54192,5208593.1795√蒲豐投針(3/5)xy概率P=2l/(pd)，5蒲豐投針(4/5)關(guān)于蒙特卡羅方法的分析和總結(jié)基本思想確定所求問(wèn)題的解是某事件的概率(或某隨機(jī)變量的數(shù)學(xué)期望、或與概率/數(shù)學(xué)期望有關(guān)的量，如

p

)通過(guò)試驗(yàn)方法，得出事件發(fā)生的頻率(或該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值，如

P)，求解數(shù)學(xué)期望與概率：當(dāng)隨機(jī)變量的取值僅為1或0時(shí)，它的數(shù)學(xué)期望就是事件的概率；反之亦然數(shù)學(xué)期望與算術(shù)平均值用隨機(jī)試驗(yàn)的方法計(jì)算積分，將積分看作服從分布密度函數(shù)f(r)的隨機(jī)變量

g(r)的數(shù)學(xué)期望通過(guò)試驗(yàn)，得到N個(gè)觀察值

r1,

r2,

…,rN

(從

f(r)

中抽取N個(gè)子樣

r1,

r2,

…,rN

)，求

g

(r)

的算術(shù)平均值√蒲豐投針(4/5)關(guān)于蒙特卡羅方法的分析和總結(jié)通過(guò)試驗(yàn)，得到6蒲豐投針(5/5)試驗(yàn)方法和次數(shù)試驗(yàn)方法不一定可行精確的近似解需要巨量的試驗(yàn)次數(shù)，但人工方法作大量的試驗(yàn)相當(dāng)困難，甚至是不可能的用計(jì)算機(jī)模擬隨機(jī)試驗(yàn)過(guò)程，完成巨量的次數(shù)抽象x的分布密度函數(shù)：q

的分布密度函數(shù)：產(chǎn)生任意的(x,q)

=

由

f1(x)抽樣

x

+

由

f2(q)抽樣

q對(duì)應(yīng)的隨機(jī)變量：數(shù)學(xué)期望與算術(shù)平均值新的問(wèn)題：誤差？收斂？√蒲豐投針(5/5)試驗(yàn)方法和次數(shù)q的分布密度函數(shù)：產(chǎn)生任意7收斂性、誤差和優(yōu)缺點(diǎn)(1/4)收斂性求解：以隨機(jī)變量X

的簡(jiǎn)單子樣X(jué)1,

X2,…,

XN

的算術(shù)平均值，作為求解的近似值近似值的收斂性大數(shù)定理：當(dāng)試驗(yàn)次數(shù)足夠多時(shí)，事件發(fā)生的頻率無(wú)窮接近于該事件發(fā)生的概率如果

X1,…,

XN

獨(dú)立分布，且期望值有限(

E(X)<

)，那么

當(dāng)隨機(jī)變量X

的簡(jiǎn)單子樣數(shù)N充分大時(shí)，其算術(shù)平均值以概率

1

收斂于期望值

E(X)√收斂性、誤差和優(yōu)缺點(diǎn)(1/4)收斂性近似值的收斂性當(dāng)隨機(jī)8收斂性、誤差和優(yōu)缺點(diǎn)(2/4)誤差概率論的中心極限定理：如果隨機(jī)變量序列X1,

X2,…,

XN獨(dú)立分布，且具有有限非零的方差s

2，即其中的

f(x)是X的分布密度函數(shù)，則當(dāng)N足夠大時(shí)(蒙特卡羅方法)不等式的概率約

1-a誤差定義為，收斂速度為a0.5000.0500.003la0.6741.9602.968√收斂性、誤差和優(yōu)缺點(diǎn)(2/4)誤差其中的f(x)是9收斂性、誤差和優(yōu)缺點(diǎn)(3/4)蒙特卡羅方法的誤差為概率誤差均方差

s

是未知的，估計(jì)值為減少誤差的技巧(在確定的置信度

a前提下)誤差

e

與試驗(yàn)次數(shù)的開(kāi)根號(hào)

N1/2

成反比：精度一個(gè)數(shù)量級(jí)，次數(shù)N兩個(gè)數(shù)量級(jí)——巨大的代價(jià)誤差

e

與均方差

s

成正比：精度一個(gè)數(shù)量級(jí)，均方差

s

一個(gè)數(shù)量級(jí)——可接受的代價(jià)效率降低方差增加觀察子樣的時(shí)間固定時(shí)間內(nèi)樣本數(shù)減少代價(jià)增加蒙特卡羅方法中的效率：由均方差

s

和觀察一個(gè)子樣的費(fèi)用(計(jì)算機(jī)時(shí))c衡量

=

s2

c√收斂性、誤差和優(yōu)缺點(diǎn)(3/4)蒙特卡羅方法的誤差為概率誤差減10收斂性、誤差和優(yōu)缺點(diǎn)(4/4)優(yōu)點(diǎn)較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程受幾何條件限制小收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力誤差容易確定程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)缺點(diǎn)收斂速度慢誤差具有概率性在某些問(wèn)題中，計(jì)算結(jié)果與系統(tǒng)大小有關(guān)主要應(yīng)用范圍粒子輸運(yùn)，統(tǒng)計(jì)物理，典型數(shù)學(xué)，真空技術(shù)，激光技術(shù)，醫(yī)學(xué)，生物，探礦等√收斂性、誤差和優(yōu)缺點(diǎn)(4/4)優(yōu)點(diǎn)√11任意分布的隨機(jī)數(shù)(1/11)隨機(jī)數(shù)(蒙特卡羅方法的關(guān)鍵)基本的：均勻地在

(0,1)

內(nèi)分布任意的：非均勻地在

(a,b)

內(nèi)分布方法：直接的(離散、連續(xù))、舍選的(簡(jiǎn)單、乘、乘加)直接抽樣方法離散型：產(chǎn)生隨機(jī)量

x

的抽樣值

xi

，概率為Pi

(i=

1,

2,…)方法：先計(jì)算x的分布函數(shù)

對(duì)于產(chǎn)生的R，如果滿足Fi-1<

R

Fi，則令抽樣x=

xi

證明：√任意分布的隨機(jī)數(shù)(1/11)隨機(jī)數(shù)(蒙特卡羅方法的關(guān)鍵)方法12任意分布的隨機(jī)數(shù)(2/11)例：二項(xiàng)分布

直接抽樣方法：

產(chǎn)生R，如果滿足，則令x

=

n例：泊松分布

直接抽樣方法：

產(chǎn)生R，如果滿足，則令x=n√任意分布的隨機(jī)數(shù)(2/11)例：二項(xiàng)分布直接抽樣13任意分布的隨機(jī)數(shù)(3/11)連續(xù)型：產(chǎn)生隨機(jī)量

x

的抽樣值，概率密度函數(shù)為f(x)方法：先計(jì)算x的分布函數(shù)

對(duì)于隨機(jī)數(shù)R'，解方程

R'

=

F(x')，得到

x'

=

F-1(R')，則令x=

x'

證明：√任意分布的隨機(jī)數(shù)(3/11)連續(xù)型：產(chǎn)生隨機(jī)量x的抽樣值14任意分布的隨機(jī)數(shù)(4/11)例：均勻介質(zhì)中，粒子運(yùn)動(dòng)的自由程S是隨機(jī)量，其概率密度函數(shù)為

f(S)

=

Se

-SS；其中的S

=

ns是宏觀總截面，s是原子截面，n為介質(zhì)中的原子數(shù)密度分布函數(shù)：介質(zhì)中粒子的碰撞過(guò)程

隨機(jī)量

R中的抽樣過(guò)程例：散射方位角余弦分布√任意分布的隨機(jī)數(shù)(4/11)例：均勻介質(zhì)中，粒子運(yùn)動(dòng)的自由程15舍任意分布的隨機(jī)數(shù)(5/11)舍選抽樣方法直接法的特點(diǎn)：簡(jiǎn)單，但分布函數(shù)的反函數(shù)不一定有解析解簡(jiǎn)單分布定理：如果Z是

[a,b]

上均勻分布的隨機(jī)數(shù)，那么利用條件f(Z)/MR

(M是f(Z)

的上界)選出的Z將是

[a,b]

上概率密度函數(shù)為f(Z)

的隨機(jī)數(shù)

證明：如右圖所示dZ(Z',R)RZf(Z)

/

Mab10Z為橫坐標(biāo)，R為縱坐標(biāo)，實(shí)曲線為函數(shù)f(Z)/MR在(0,1)內(nèi)、Z在[a,b]內(nèi)均勻分布隨機(jī)點(diǎn)(Z',R)在虛框內(nèi)均勻分布隨機(jī)點(diǎn)落入窄條的概率=兩面積之比選已知概率密度函數(shù)抽樣方法√舍任意分布的隨機(jī)數(shù)(5/11)舍選抽樣方法證明：如右圖所16任意分布的隨機(jī)數(shù)(6/11)簡(jiǎn)單分布抽樣方法的流程產(chǎn)生隨機(jī)數(shù)x和R計(jì)算Z

=

a+(b-a)xf(Z)

/

MR?Z'

=

ZYN效率(有一部分隨機(jī)數(shù)被舍棄)舍dZ(Z',R)RZf(Z)

/

Mab10選例：受限的散射方位角余弦分布

效率√任意分布的隨機(jī)數(shù)(6/11)簡(jiǎn)單分布抽樣方法的流程產(chǎn)生隨機(jī)數(shù)17任意分布的隨機(jī)數(shù)(7/11)乘分布：f(x)

有銳鋒時(shí)效率很低，需要改進(jìn)的簡(jiǎn)單分布方法：將函數(shù)寫成

f(Z)

=

h(Z)

f1(Z)，由容易抽樣的

f1(Z)抽樣出Z，代入h(Z)，如果滿足條件

h(Z)/MR

(M是

h(Z)

的上界)，那么得到概率密度函數(shù)為f(Z)

的抽樣值舍

證明：如右圖所示dF1(F1',R)RF1(Z)h(Z)

/

M110f1(Z)的分布函數(shù)F1(Z)

為橫坐標(biāo)，R為縱坐標(biāo)，實(shí)曲線為函數(shù)

h(Z)/MR和F1

在(0,1)內(nèi)均勻分布隨機(jī)點(diǎn)(F1',R)在虛框內(nèi)均勻分布隨機(jī)點(diǎn)落入窄條的概率=兩面積之比選效率√任意分布的隨機(jī)數(shù)(7/11)乘分布：f(x)有銳鋒時(shí)效率很18任意分布的隨機(jī)數(shù)(8/11)乘分布抽樣方法的流程產(chǎn)生隨機(jī)數(shù)F1和R利用F1，由f1(Z)抽樣Zh(Z)

/

MR?Z'

=

ZYN例：半正態(tài)分布，概率密度函數(shù)√任意分布的隨機(jī)數(shù)(8/11)乘分布抽樣方法的流程產(chǎn)生隨機(jī)數(shù)利19任意分布的隨機(jī)數(shù)(9/11)乘加分布方法：如果隨機(jī)量x具有以下乘加形式

其中f1(x)和f2(x)滿足密度函數(shù)的要求，即

并且h1(x)和h2(x)滿足

那么從f(x)的以下變形，可以得到乘加分布

其中

和

為加權(quán)因子

證明：參照直接法和舍選法(略)√任意分布的隨機(jī)數(shù)(9/11)乘加分布其中f1(x)和20任意分布的隨機(jī)數(shù)(10/11)乘加分布抽樣方法的流程產(chǎn)生隨機(jī)數(shù)g,

F,R利用F由f1(x)抽樣x1gg1

?YNh1(x1)/M1

R?利用F由f2(x)抽樣x2h2(x2)/M2

R?YYNNx

=

x1x

=

x2效率√任意分布的隨機(jī)數(shù)(10/11)乘加分布抽樣方法的流程產(chǎn)生隨機(jī)21任意分布的隨機(jī)數(shù)(11/11)例：散射光子能量抽樣。能量為

E0

的入射光子，經(jīng)原子散射后的能量E按某種概率分布。如果令

x

=

E0/E，那么概率密度函數(shù)

f(x)

為(其中K=K(E0)為歸一因子)√任意分布的隨機(jī)數(shù)(11/11)例：散射光子能量抽樣。能量為22粒子輸運(yùn)問(wèn)題(1/7)蒙特卡羅模擬粒子輸運(yùn)是隨機(jī)過(guò)程，運(yùn)動(dòng)規(guī)律是大量粒子運(yùn)動(dòng)的統(tǒng)計(jì)蒙特卡羅模擬：模擬一定數(shù)量粒子在介質(zhì)中的運(yùn)動(dòng)，再現(xiàn)粒子運(yùn)動(dòng)的統(tǒng)計(jì)規(guī)律例(平板介質(zhì)模型)：由單一物質(zhì)組成的均勻介質(zhì)，厚度為H，能量為

E00

的平行光子束垂直射入板內(nèi)，求光子對(duì)板的投射率H抽樣：自由程、作用類型、散射能量、散射方向自由程抽樣(PPT15)：粒子運(yùn)動(dòng)的自由程S是隨機(jī)量，其概率密度函數(shù)為

f(S)

=

Se

-SS；其中的S

=

ns是宏觀總截面，s是原子截面，n為介質(zhì)中的原子數(shù)密度分布函數(shù)：F(S')

=

1

-

e-SS'粒子的碰撞過(guò)程

隨機(jī)量

R中的抽樣過(guò)程：S'

=

-lnR'

/

S√粒子輸運(yùn)問(wèn)題(1/7)蒙特卡羅模擬例(平板介質(zhì)模型)：由單一23粒子輸運(yùn)問(wèn)題(2/7)作用類型抽樣光子-介質(zhì)：散射(康普頓散射)和吸收(光電效應(yīng))

散射截面

吸收截面

總截面：St

=

Ss

+

Sa光子隨機(jī)地被散射或吸收的過(guò)程隨機(jī)量

R中的抽樣過(guò)程：如果

Ss

/

St

R，則散射，否則為吸收√粒子輸運(yùn)問(wèn)題(2/7)作用類型抽樣吸收截面總截面：S24粒子輸運(yùn)問(wèn)題(3/7)散射能量抽樣(PPT22)：能量為

E0

的入射光子，經(jīng)原子散射后的能量E按某種概率分布。如果令

x

=

E0/E，那么概率密度函數(shù)

f(x)

為(其中K=K(E0)為歸一因子)√粒子輸運(yùn)問(wèn)題(3/7)散射能量抽樣(PPT22)：能量為E25粒子輸運(yùn)問(wèn)題(4/7)散射方向抽樣坐標(biāo)系：以入射方向(由

q

和

f

確定)為參考系，散射方向

由散射角

q'

和散射方位角

f'

確定xyz流程圖輸入g0,E0,E計(jì)算m

=1+1/E0-1/E計(jì)算g產(chǎn)生隨機(jī)數(shù)R計(jì)算f'

=

2pR√粒子輸運(yùn)問(wèn)題(4/7)散射方向抽樣坐標(biāo)系：以入射方向(26粒子輸運(yùn)問(wèn)題(5/7)直接模擬方法模擬粒子在介質(zhì)中的真實(shí)物理過(guò)程粒子在介質(zhì)中的狀態(tài)：空間位置，能量和運(yùn)動(dòng)方向碰撞點(diǎn)的狀態(tài)參數(shù)

(Zm,

Em,

gm)

表示從源出發(fā)的粒子在介質(zhì)中經(jīng)過(guò)m次碰撞后的狀態(tài)ZHOZ0E0q0Z1E1q1Z2E2q2粒子的運(yùn)動(dòng)過(guò)程

=

碰撞點(diǎn)的狀態(tài)序列

模擬粒子的運(yùn)動(dòng)過(guò)程

=確定狀態(tài)序列的問(wèn)題√粒子輸運(yùn)問(wèn)題(5/7)直接模擬方法碰撞點(diǎn)的狀態(tài)參數(shù)(Zm,27粒子輸運(yùn)問(wèn)題(6/7)記錄和計(jì)算的內(nèi)容穿透率和誤差估計(jì)(設(shè)N為入射粒子數(shù)，N1

為透射數(shù))透射粒子的能量和方向分布離散化能量：Emin

=

E0

<E1

<

<Ei

<

<

EI離散化角度：0

=

q0

<

q1

<

<

qj

<

<

qJ

=

p

/

2能量和方向分布粒子的軌跡粒子運(yùn)動(dòng)關(guān)于Z軸對(duì)稱，只需要記錄(Zim,

qim)流程圖√粒子輸運(yùn)問(wèn)題(6/7)記錄和計(jì)算的內(nèi)容透射粒子的能量和方向分28粒子輸運(yùn)問(wèn)題(7/7)√輸入N,

H,E00,

Emin,g00光子i的初態(tài)m=0,Z=0,

E=E00,g=g00利用E計(jì)算截面Ss,St利用Ss抽樣自由程S計(jì)算碰撞點(diǎn)位置ZZ+Sg利用Ss,St抽樣作用類型計(jì)數(shù)器N1N1+1利用E0

E抽樣散射能量E利用g0g,

E0

,

E抽樣散射方向gmm+1ii+1輸出P=N1/N,

s,

DPZ

0

?類型=吸收

?E

Emin

?YNNYNYm>M?NYi>N?NYZ

H

?NY粒子輸運(yùn)問(wèn)題(7/7)√輸入N,H,光子i的初態(tài)m=029隨機(jī)過(guò)程模擬(1/3)蒙特卡羅模擬隨機(jī)過(guò)程的兩種情況已知隨機(jī)過(guò)程的概率分布函數(shù)隨機(jī)過(guò)程的統(tǒng)計(jì)特征建立與隨機(jī)過(guò)程的模型，形成隨機(jī)量，并使其數(shù)字特征(概率、平均值、方差等)是問(wèn)題的解由已知的概率分布函數(shù)進(jìn)行大量的抽樣統(tǒng)計(jì)處理抽樣結(jié)果，給出問(wèn)題的解及其誤差已知隨機(jī)現(xiàn)象的觀測(cè)數(shù)據(jù)概率分布函數(shù)分析現(xiàn)象的特征，假設(shè)隨機(jī)量服從某種分布，建立概率模型由觀測(cè)數(shù)據(jù)推斷分布中的參數(shù)按推斷的分布進(jìn)行大量的抽樣比較抽樣值和觀測(cè)值，根據(jù)誤差，判斷假設(shè)的準(zhǔn)確性√隨機(jī)過(guò)程模擬(1/3)蒙特卡羅模擬隨機(jī)過(guò)程的兩種情況√30隨機(jī)過(guò)程模擬(2/3)例：a

粒子衰變的蒙特卡羅模擬隨機(jī)現(xiàn)象的觀測(cè)數(shù)據(jù)概率分布函數(shù)觀測(cè)數(shù)據(jù)：每隔

Dt<<半衰期

測(cè)量一次放射的

a

粒子數(shù)，共測(cè)量

N=2608

次，測(cè)得

k

個(gè)粒子的次數(shù)為

nkknkknk05762731203