1、计算物理计算物理http:/125.217.162.13/lesson/ComputationalPhysicshttp:/125.217.162.13/lesson/ComputationalPhysics蒙特卡罗方法蒙特卡罗方法蒙特卡罗方法蒙特卡罗方法n蒲丰投针蒲丰投针n收敛性、误差和优缺点收敛性、误差和优缺点n任意分布的随机数任意分布的随机数n粒子输运问题粒子输运问题n随机过程模拟随机过程模拟n梅氏抽样梅氏抽样蒲丰投针蒲丰投针(1/(1/5 5) )n蒙特卡罗方法蒙特卡罗方法n又称随机抽样技巧或统计试验方法又称随机抽样技巧或统计试验方法n以概率统计理论为基础的以概率统计理论为基础的n能够
2、比较逼真地描述事物的特点及物理实验过程能够比较逼真地描述事物的特点及物理实验过程n解决一般数值方法难以解决的问题解决一般数值方法难以解决的问题n随着电子计算机的发展而发展随着电子计算机的发展而发展n首先在核武器的试验与研制中得到了应用首先在核武器的试验与研制中得到了应用n蒲丰投针蒲丰投针n法国数学家蒲丰的法国数学家蒲丰的17771777年出版的著作：年出版的著作：“在平面上画有一组间距为在平面上画有一组间距为 d 的平行线，的平行线，将一根长度为将一根长度为 l ( (l d) ) 的针任意掷在这个的针任意掷在这个平面上，求此针与平行线中任一条相交平面上，求此针与平行线中任一条相交的概率。的概
3、率。”蒲丰投针蒲丰投针( (2 2/ /5 5) )n步骤步骤n在桌面上画出在桌面上画出间距为间距为 2 2d 的平行线的平行线n准备长度为准备长度为 2 2l ( (l d) ) 的针的针n向向桌面桌面随机随机投投针针n如果针与平行线相交，则计数器如果针与平行线相交，则计数器 n 加加 1 1n计算：计数器计算：计数器 n 与总与总投投针数针数 N 的比例的比例( (视作相交概率视作相交概率 P ) )n概率分析概率分析 P = = ? ?n各条平行线地位等同，仅考虑某条各条平行线地位等同，仅考虑某条平行线附近的情况平行线附近的情况n平行线方向的平行线方向的 x 坐标坐标对对概率概率没影响没
4、影响n针的中点的针的中点的 y 坐标在线之间等概率坐标在线之间等概率落入落入( (均匀分布在均匀分布在 0,0, d ) )，仅当，仅当 y l 才可能才可能针针- -线线相交相交n针针- -线的夹角线的夹角 q q 均匀分布在均匀分布在 0,0, p p ，q q 与与 y 独立独立xy蒲丰投针蒲丰投针(3/(3/5 5) )xyn概率概率 P = = 2 2l / / ( (p p d) ) ，可计算圆周率，可计算圆周率实验者实验者时间时间针长针长总投数总投数相交数相交数p p 值值WolfWolf185018500.80.85,0005,0002,5322,532 3.15963.159
5、6SmithSmith185518550.60.63,2043,2041,2181,218 3.15543.1554De Morgan,CDe Morgan,C186018601.01.0600600382382 3.1373.137FoxFox188418840.750.751,0301,030489489 3.15953.1595LezzeriniLezzerini190119010.830.833,4083,4081,8081,808 3.14159293.1415929ReinaReina192519250.54190.54192,5202,520859859 3.17953.1795
6、蒲丰投针蒲丰投针(4/(4/5 5) )n关于关于蒙特卡罗方法的蒙特卡罗方法的分析和总结分析和总结n基本思想基本思想n确定确定所求问题的解是某事件的概率所求问题的解是某事件的概率( (或某随机变量或某随机变量的数学期望的数学期望、或与概率或与概率/ /数学期望有关的量数学期望有关的量，如，如 p p ) )n通过试验方法，得出事件发生的频率通过试验方法，得出事件发生的频率( (或该随机变或该随机变量若干个具体观察值的算术平均值量若干个具体观察值的算术平均值，如，如 P ) )，求解，求解n数学期望与概率：当随机变量的取值仅数学期望与概率：当随机变量的取值仅为为 1 1 或或 0 0 时时，它的
7、数学期望就是事件的概率；反之亦然它的数学期望就是事件的概率；反之亦然n数学期望与算术平均值数学期望与算术平均值n用随机试验的方法计算积分，将积分看作服从分布用随机试验的方法计算积分，将积分看作服从分布密度函数密度函数 f ( (r) ) 的随机变量的随机变量 g ( (r) ) 的的数学期望数学期望 =rrfrggd)()(n通过试验，得到通过试验，得到 N 个观察值个观察值 r1 1, , r2 2, , , rN ( (从从 f ( (r) ) 中中抽取抽取 N 个个子样子样 r1 1, , r2 2, , , rN ) )，求，求 g ( (r) ) 的算术平均值的算术平均值NrggNi
8、iN/ )(1=蒲丰投针蒲丰投针(5/(5/5 5) )n试验方法和次数试验方法和次数n试验方法不一定可行试验方法不一定可行n精确的近似解需要巨量的试验次数，但人工方法作精确的近似解需要巨量的试验次数，但人工方法作大量的试验相当困难，甚至是不可能的大量的试验相当困难，甚至是不可能的n用计算机模拟随机试验过程，完成巨量的次数用计算机模拟随机试验过程，完成巨量的次数n抽象抽象nx 的分布密度函数：的分布密度函数：其它当 , 0 ;0 ,/1)(1dxdxf=nq q 的分布密度函数：的分布密度函数：其它当 , 0 ;0 ,/1)(2pqpq=fn产生任意的产生任意的 ( (x, ,q q ) )
9、= = 由由 f1 1( (x) )抽样抽样 x + + 由由 f2 2( (q q) )抽样抽样 q q=其它, 0sin, 1),(qqlxxsn对应的对应的随机变量：随机变量：n数学期望与算术平均值数学期望与算术平均值=NiNiiNdslPxsNsdlxfxfxsP1212),(1 ,2dd)()(),(pqpqqqn新的问题：误差？收敛？新的问题：误差？收敛？收敛性、误差和收敛性、误差和优缺优缺点点(1/(1/4 4) )n收敛性收敛性n求解：以求解：以随机变量随机变量 X 的简单的简单子样子样 X1 1, , X2 2, , , , XN N 的算术平均的算术平均值，作为求解的近似值
10、值，作为求解的近似值=NiiNXNX11n近似值的收敛性近似值的收敛性n大数定理：当试验次数足够多时，事件发生的频率无大数定理：当试验次数足够多时，事件发生的频率无穷接近于该事件发生的概率穷接近于该事件发生的概率n如果如果 X1 1, , , , XN N 独立分布，且期望值有限独立分布，且期望值有限( ( E( (X) ) )，那么那么()1)(lim=XEXPNN 当当随机变量随机变量 X 的简单子样数的简单子样数 N 充分大时，其算术平均充分大时，其算术平均值以概率值以概率 1 1 收敛于期望值收敛于期望值 E( (X) )收敛性、误差和收敛性、误差和优缺优缺点点(2/(2/4 4) )
11、n误差误差n概率论的中心极限定理：如果随机变量序列概率论的中心极限定理：如果随机变量序列 X1 1, , X2 2, , , , XN N 独立分布，且具有有限非零的方差独立分布，且具有有限非零的方差 s s 2 2，即，即=xxfXExd)()(022s 其中的其中的 f( (x) ) 是是 X 的分的分布密度函数，则布密度函数，则=xxtNNtexXEXNP 2/d21| )(|lim2psn当当 N 足够大时足够大时( (蒙特卡罗方法蒙特卡罗方法) )n不等式不等式 的概率约的概率约 11a aNXEXN/| )(|san误差定义为误差定义为 ，收敛速度为，收敛速度为N/sa=)( NO
12、a a0.5000.5000.0500.0500.0030.003 a a0.6740.6741.9601.9602.9682.968收敛性、误差和收敛性、误差和优缺优缺点点(3/(3/4 4) )n蒙特卡罗方法的误差为概率误差蒙特卡罗方法的误差为概率误差n均方差均方差 s s 是未知的，估计值为是未知的，估计值为222112)1(1XXXNXNNiiNii=sn减少误差的技巧减少误差的技巧( (在确定的置信度在确定的置信度 a a 前提前提下下) )N/sa=n误差误差 与试验次数的开根号与试验次数的开根号 N1/21/2 成反比：精度成反比：精度 一个数一个数量级，次数量级，次数 N 两个
13、数量级两个数量级巨大的代价巨大的代价n误差误差 与均方差与均方差 s s 成正比：精度成正比：精度 一个数量级，均方差一个数量级，均方差 s s 一个数量级一个数量级可接受的代价可接受的代价n效率效率n降低方差降低方差增加观察子样的时间增加观察子样的时间固定时间内样本数固定时间内样本数减少减少代价增加代价增加n蒙特卡罗方法中的效率：由均方差蒙特卡罗方法中的效率：由均方差 s s 和观察一个子样和观察一个子样的费用的费用( (计算机时计算机时) ) c 衡量衡量 = = s s 2 2 c收敛性、误差和收敛性、误差和优缺优缺点点(4/(4/4 4) )n优点优点n较逼真地描述具有随机性质的事物的
14、特点及物理实验过程较逼真地描述具有随机性质的事物的特点及物理实验过程n受几何条件限制小受几何条件限制小n收敛速度与问题的维数无关收敛速度与问题的维数无关n具有同时计算多个方案与多个未知量的能力具有同时计算多个方案与多个未知量的能力n误差容易确定误差容易确定n程序结构简单，易于实现程序结构简单，易于实现n缺点缺点n收敛速度慢收敛速度慢n误差具有概率性误差具有概率性n在某些问题中，计算结果与系统大小有关在某些问题中，计算结果与系统大小有关n主要应用范围主要应用范围n粒子输运，统计物理，典型数学，真空技术，激光技术，粒子输运，统计物理，典型数学，真空技术，激光技术，医学，生物，探矿等医学，生物，探矿
15、等任意分布的随机数任意分布的随机数(1/(1/1111) )n随机数随机数( (蒙特卡罗方法的关键蒙特卡罗方法的关键) )n基本的：基本的：均匀均匀地在地在 (0, 1)(0, 1) 内分布内分布n任意的：任意的：非均匀非均匀地在地在 ( (a, , b) ) 内分布内分布n方法：直接的方法：直接的( (离散、连续离散、连续) )、舍选的、舍选的( (简单、乘、乘加简单、乘、乘加) )n直接抽样方法直接抽样方法n离散型：产生随机量离散型：产生随机量 x 的抽样值的抽样值 xi ，概率为，概率为 Pi ( (i = = 1, 1, 2, 2, ) )n方法：先计算方法：先计算 x 的分布函数的分
16、布函数 =ijjiPF1 对于产生的对于产生的 R，如果满足，如果满足 Fi11 R Fi，则令抽样，则令抽样 x = = xi 证明：证明：任意分布的随机数任意分布的随机数(2/(2/1111) )n例：二项分布例：二项分布nNnnNnppCP=)1 ( 直接抽样方法直接抽样方法： 产生产生 R，如果满足，如果满足 ，则令，则令 x = = n=niiniiPRP010n例：泊松分布例：泊松分布! nePnn= 直接抽样方法直接抽样方法： 产生产生 R，如果满足，如果满足 ，则令，则令x= =n=niiniiieRi010! !任意分布的随机数任意分布的随机数(3/(3/1111) )n连续
17、型：产生随机量连续型：产生随机量 x 的抽样值，概率密度函数为的抽样值，概率密度函数为 f ( (x) )n方法：先计算方法：先计算 x 的分布函数的分布函数=xxxfxF d)()( 对于随机数对于随机数 R，解方程，解方程 R = = F( (x) )，得到，得到 x = = F11( (R) )，则令则令 x = = x 证明：证明：任意分布的随机数任意分布的随机数(4/(4/1111) )n例：均匀介质中，粒子运动的自由程例：均匀介质中，粒子运动的自由程 S 是随机量，其是随机量，其概率密度函数为概率密度函数为 f ( (S) ) = = S Se SSS；其中的；其中的 S S =
18、= ns s 是宏观总是宏观总截面，截面，s s 是原子截面，是原子截面，n 为介质中的原子数密度为介质中的原子数密度n分布函数：分布函数：n介质中粒子的碰撞过程介质中粒子的碰撞过程 随机量随机量 R 中的抽样过程中的抽样过程n例：散射方位角余弦分布例：散射方位角余弦分布)cos(11 ,111)(2R xxxxf=pp舍舍任意分布的随机数任意分布的随机数(5/(5/1111) )n舍选抽样方法舍选抽样方法直接法的特点：简单，但分布函数的反函数不一定有解析解直接法的特点：简单，但分布函数的反函数不一定有解析解n简单分布简单分布n定理：如果定理：如果 Z 是是 a, , b 上均匀分布的随机数，
19、那么利用上均匀分布的随机数，那么利用条件条件 f( (Z)/ )/M R ( (M 是是 f( (Z) ) 的上界的上界) )选出的选出的 Z 将是将是 a, , b 上概率密度函数为上概率密度函数为 f( (Z) ) 的随机数的随机数 证明：如右图所示证明：如右图所示dZ( (Z, , R) )RZf ( (Z) ) / / Mab1 10 0nZ 为横坐标，为横坐标，R 为纵坐标，实曲线为为纵坐标，实曲线为函数函数 f( (Z)/ )/MnR 在在 ( (0, , 1) ) 内、内、Z 在在 a, , b 内均匀分布内均匀分布 随机点随机点 ( (Z, , R) ) 在虚框内在虚框内均匀分
20、布均匀分布n随机点落入窄条的概率随机点落入窄条的概率 = = 两面积之比两面积之比选选n已知概率密度函数已知概率密度函数抽样方法抽样方法任意分布的随机数任意分布的随机数(6/(6/1111) )n简单分布抽样方法的流程简单分布抽样方法的流程产生随机数产生随机数x 和和 R计算计算Z = = a + (+ (b a) ) x f ( (Z) ) / / M R ? ?Z = = ZYNn效率效率( (有一部分随机数被舍弃有一部分随机数被舍弃) )舍舍dZ( (Z, , R) )RZf ( (Z) ) / / Mab1 10 0选选n例：受限的散射方位角余弦分布例：受限的散射方位角余弦分布2121
21、 ,113)(2=xxxfp 效率效率任意分布的随机数任意分布的随机数( (7 7/ /1111) )n乘分布：乘分布：f( (x) ) 有锐锋时效率很低，需要改进的简单分布有锐锋时效率很低，需要改进的简单分布n方法：将函数写成方法：将函数写成 f( (Z) ) = = h( (Z) ) f1 1( (Z) )，由容易抽样的，由容易抽样的 f1 1( (Z) ) 抽样出抽样出 Z，代入，代入 h( (Z) )，如果满足条件如果满足条件 h( (Z) )/ /M R ( (M 是是 h( (Z) ) 的上界的上界) )，那么得到概率密度函数为，那么得到概率密度函数为 f( (Z) ) 的抽样值的
22、抽样值舍舍 证明：如右图所示证明：如右图所示dF1 1( (F1 1, , R) )RF1 1( (Z) )h ( (Z) ) / / M1 11 10 0nf1 1( (Z) ) 的分布函数的分布函数 F1 1( (Z) ) 为横坐标，为横坐标，R 为纵坐标，实曲线为函数为纵坐标，实曲线为函数 h( (Z) )/ /MnR 和和 F1 1 在在 ( (0, , 1) ) 内均匀分布内均匀分布 随机随机点点 ( (F1 1, , R) ) 在虚框内在虚框内均匀分布均匀分布n随机点落入窄条的概率随机点落入窄条的概率 = = 两面积之比两面积之比选选n效率效率任意分布的随机数任意分布的随机数(8/
23、(8/1111) )n乘分布抽样方法的流程乘分布抽样方法的流程产生随机数产生随机数F1 和和 R利用利用 F1，由由 f1 1( (Z) ) 抽样抽样 Z h ( (Z) ) / / M R ? ?Z = = ZYNn例：半正态分布，概率密度函数例：半正态分布，概率密度函数0 ,2)(221=xexfxp任意分布的随机数任意分布的随机数(9/(9/1111) )n乘加分布乘加分布n方法：如果随机量方法：如果随机量 x 具有以下乘加形式具有以下乘加形式bxaxfxhxfxhxf+= ),()()()()(2211 其中其中 f1 1( (x) ) 和和 f2 2( (x) ) 满足密度函数的要求
24、，即满足密度函数的要求，即=babaxxfxxfxfxf 2 121d)(d)(1 , 0)(),( 并且并且 h1 1( (x) ) 和和 h2 2( (x) ) 满足满足2211)( ,)(MxhMxh 那么从那么从 f( (x) ) 的以下变形，可以得到的以下变形，可以得到乘加分布乘加分布)()()()()()(2222111121xfMxhgxfMxhgMMxf+= 其中其中 和和 为加权因子为加权因子2111MMMg+=2122MMMg+= 证明：参照直接法和舍选法证明：参照直接法和舍选法( (略略) )任意分布的随机数任意分布的随机数(10/(10/1111) )n乘加分布抽样方法
25、的流程乘加分布抽样方法的流程产生随机数产生随机数g, , F, , R利用利用 F由由 f1 1( (x) ) 抽样抽样 x1 1 g g1 1 ? ?YN h1 1( (x1 1)/ )/M1 1 R ? ?利用利用 F由由 f2 2( (x) ) 抽样抽样 x2 2 h2 2( (x2 2)/ )/M2 2 R ? ?YYNNx = = x1 1x = = x2 2n效率效率211MM +=任意分布的随机数任意分布的随机数(11/(11/1111) )n例：散射光子能量抽样。能量为例：散射光子能量抽样。能量为 E0 0 的入射光子，经原的入射光子，经原子散射后的能量子散射后的能量 E 按某
26、种概率分布。如果令按某种概率分布。如果令 x = = E0 0/ /E，那么概率密度函数那么概率密度函数 f( (x) ) 为为( (其中其中 K = = K( (E0 0) ) 为归一因子为归一因子) )121 ,) 1(1 1)1(1)(0322200+=ExxxxExEKxf粒子输运问题粒子输运问题(1/(1/7 7) )n蒙特卡罗模拟蒙特卡罗模拟n粒子输运是随机过程，运动规律是大量粒子运动的统计粒子输运是随机过程，运动规律是大量粒子运动的统计n蒙特卡罗模拟：模拟一定数量粒子在介质中的运动，再现蒙特卡罗模拟：模拟一定数量粒子在介质中的运动，再现粒子运动的统计规律粒子运动的统计规律n例例(
27、 (平板介质模型平板介质模型) )：由单一物质组成的均：由单一物质组成的均匀介质，厚度为匀介质，厚度为 H，能量为，能量为 E0000 的平行光的平行光子束垂直射入板内，求光子对板的投射率子束垂直射入板内，求光子对板的投射率Hn抽样：自由程、作用类型、散射能量、散射方向抽样：自由程、作用类型、散射能量、散射方向n自由程抽样自由程抽样(PPT15)(PPT15)：粒子运动的自由程粒子运动的自由程 S 是随机量，是随机量，其概率密度函数为其概率密度函数为 f ( (S) ) = = S Se SSS；其中的；其中的 S S = = ns s 是宏观是宏观总截面，总截面，s s 是原子截面，是原子截
28、面，n 为介质中的原子数密度为介质中的原子数密度n分布函数：分布函数：F( (S) ) = = 1 1 eSSSn粒子的碰撞过程粒子的碰撞过程 随机量随机量 R 中中的的抽样过程：抽样过程：S = = lnR / / S S粒子输运问题粒子输运问题(2/(2/7 7) )n作用类型抽样作用类型抽样n光子光子- -介质：散射介质：散射( (康普顿散射康普顿散射) )和吸收和吸收( (光电效应光电效应) ) 散射截面散射截面+=S1,)2/12(ln1953. 01),5/2621 (5207. 0)(12EEEEEEEs 吸收截面吸收截面=S1,107916. 01,102985. 0)(142
29、/73EEEEEa 总截面：总截面：S St = = S Ss + + S San光子随机地被散射或吸收的过程光子随机地被散射或吸收的过程 随机量随机量 R 中的中的抽抽样过程：如果样过程：如果 S Ss / / S St R，则散射，否则为吸收，则散射，否则为吸收粒子输运问题粒子输运问题(3/(3/7 7) )n散射能量抽样散射能量抽样(PPT22)(PPT22)：能量为：能量为 E0 0 的入射光子，经原的入射光子，经原子散射后的能量子散射后的能量 E 按某种概率分布。如果令按某种概率分布。如果令 x = = E0 0/ /E，那么概率密度函数那么概率密度函数 f( (x) ) 为为( (
30、其中其中 K = = K( (E0 0) ) 为归一因子为归一因子) )121 ,) 1(1 1)1(1)(0322200+=ExxxxExEKxf粒子输运问题粒子输运问题(4/(4/7 7) )n散射方向抽样散射方向抽样n坐标系：以入射方向坐标系：以入射方向 ( (由由 q q 和和 f f 确定确定) )为参考系，散射方向为参考系，散射方向 由由散射角散射角 q q 和和散射方位角散射方位角 f f 确定确定0 xyz00q0fqfqfn流程图流程图输入输入g g0 0, E0 0, E计算计算m m = =1+1/+1/E0 01/1/E计算计算 g g产生随机数产生随机数 R计算计算
31、f f = = 2 2p pR粒子输运问题粒子输运问题(5/(5/7 7) )n直接模拟方法直接模拟方法n模拟粒子在介质中的真实物理过程模拟粒子在介质中的真实物理过程n粒子在介质中的状态：空间位置，能量和运动方向粒子在介质中的状态：空间位置，能量和运动方向n碰撞点的状态参数碰撞点的状态参数 ( (Zm, , Em, , g gm) ) 表示从源出发的粒子在表示从源出发的粒子在介质中经过介质中经过 m 次碰撞后的状态次碰撞后的状态ZHOZ0E0q0Z1E1q1Z2E2q2n粒子的运动过程粒子的运动过程 = = 碰撞点的状态序列碰撞点的状态序列是末态是初态MmmEEEEZZZZMMMMMM= ,0
32、 ,110110110gggg 模拟粒子的运动过程模拟粒子的运动过程 = = 确定确定状态序列的问题状态序列的问题粒子输运问题粒子输运问题(6/(6/7 7) )n记录和计算的内容记录和计算的内容n穿透率和误差估计穿透率和误差估计( (设设 N 为入射粒子数，为入射粒子数，N1 1 为透射数为透射数) )n透射粒子的能量和方向分布透射粒子的能量和方向分布n离散化能量：离散化能量：Emin = = E0 0 E1 1 Ei EIn离散化角度：离散化角度：0 0 = = q q0 0 q q1 1 q qj q qJ = = p p / / 2 2n能量和方向分布能量和方向分布n粒子的轨迹粒子的轨
33、迹n粒子运动关于粒子运动关于 Z 轴对称，只需要记录轴对称，只需要记录 ( (Zim, , q qim) )n流程图流程图粒子输运问题粒子输运问题( (7 7/ /7 7) )输入输入N, , H, ,E0000, , Emin, , g g0000光子光子 i 的的初态初态m=0,=0,Z=0,=0, E= =E0000, , g g= =g g0000利用利用 E 计算截面计算截面S Ss , S, St利用利用 S Ss 抽样抽样自由程自由程 S计算碰撞点位置计算碰撞点位置ZZ+ +S g g利用利用 S Ss , S, St 抽样抽样作用类型作用类型计数器计数器N1 1N1 1+1+1
34、利用利用 E0 0 E 抽样抽样散射能量散射能量 E利用利用 g g0 0g g, , E0 0 , , E抽样散射方向抽样散射方向 g gmm+1+1ii+1+1输出输出P= =N1 1/ /N, , s s, , D DPZ 0 0 ? ? 类型类型= =吸收吸收 ? ? E Emin ? ?YNNYNY m M ? ?NY i N ? ?NYZ H H ? ?NY随机过程模拟随机过程模拟(1/(1/3 3) )n蒙特卡罗模拟蒙特卡罗模拟随机过程随机过程的两种情况的两种情况n已知随机过程的概率分布函数已知随机过程的概率分布函数随机过程的统计特征随机过程的统计特征n建立与建立与随机过程的模型
35、，形成随机量，并使其数字特随机过程的模型，形成随机量，并使其数字特征征( (概率、平均值、方差等概率、平均值、方差等) )是问题的解是问题的解n由已知的概率分布函数进行大量的抽样由已知的概率分布函数进行大量的抽样n统计处理抽样结果，给出问题的解及其误差统计处理抽样结果，给出问题的解及其误差n已知随机现象的观测数据已知随机现象的观测数据概率分布函数概率分布函数n分析现象的特征，假设随机量服从某种分布，建立概分析现象的特征，假设随机量服从某种分布，建立概率模型率模型n由观测数据推断分布中的参数由观测数据推断分布中的参数n按推断的分布进行大量的抽样按推断的分布进行大量的抽样n比较抽样值和观测值，根据
36、误差，判断假设的准确性比较抽样值和观测值，根据误差，判断假设的准确性随机过程模拟随机过程模拟(2/(2/3 3) )n例：例：a a 粒子衰变的蒙特卡罗模拟粒子衰变的蒙特卡罗模拟n随机现象的观测数据随机现象的观测数据概率分布函数概率分布函数n观测数据：每隔观测数据：每隔 D Dt半衰期半衰期 测量一次测量一次放射的放射的 a a 粒子数，共测量粒子数，共测量 N=2608=2608 次，次，测得测得 k 个粒子的次数为个粒子的次数为 nkknkknk 0 57 6273 1203 7139 2383 8 45 3525 9 27 453210 16 5408n分析现象，建立模型分析现象，建立模
37、型nD Dt 内的衰变数内的衰变数 k 是随机事件的发生次数，是随机量是随机事件的发生次数，是随机量nk 在有限的平均值在有限的平均值 a 上下波动上下波动87. 3)1610571 (26081100+=kkNNnkkank 可看作大量可看作大量( ( N 足够大足够大) )独立独立( (每个原子核的衰变每个原子核的衰变不受其它原子核的影响不受其它原子核的影响) )试验的结果试验的结果n事件有相同的很小的概率事件有相同的很小的概率( (原子核衰变的机会相等原子核衰变的机会相等) )nk 的上述特征表明：的上述特征表明： k 近似服从泊松分布近似服从泊松分布10, 0 ,87. 3 , !/=
38、kakaePkak随机过程模拟随机过程模拟(3/(3/3 3) )n蒙特卡罗模拟蒙特卡罗模拟n产生产生 26082608 个泊松分布的随机数，统计其中数值等于个泊松分布的随机数，统计其中数值等于 k 的个数并与观测值的个数并与观测值 nk 比较比较梅氏抽样梅氏抽样(1/(1/1 1) )n特别适合于多维随机量的系统特别适合于多维随机量的系统n设设 x 是多维随机量，是多维随机量，f( (x) ) 是概率密度函数，流程是概率密度函数，流程梅氏梅氏1 1步步输入输入x, , d dx产生多维随机数产生多维随机数Rx= =x+ +d dx( (R0.5)0.5)产生随机数产生随机数r f( (x)/
39、 )/f( (x) ) r ? ?x= =x返回返回 xYN输入输入Nt, , Ng, , Nf, , d dx初始化初始化 x由由 x 走梅氏走梅氏 Nt 步步间隔间隔 Nf 抽样：抽样：由由 x 走梅氏走梅氏 Nf 步步由由 x 计算计算 u 和和 u2 2AA+ +u, , BB+ +u2 2 m Ng ? ?YNmm+1+1计算平均值计算平均值,方差方差, 误差误差主程序主程序结束结束作业作业n用蒲丰投针在计算机上计算用蒲丰投针在计算机上计算 p p 值，取值，取 d = = 4, 4, l = = 3 3n分子热运动的速率概率密度函数是分子热运动的速率概率密度函数是=kTmvxxexxxfx20 ,d2d)(2p 取取 f1 1( (x) ) = = 2 2e22x/3 /3/3 /3，找出，找出 h( (x) ) 和和 M ，用直接法对，用直接法对 f1 1( (x) ) 抽样并写出程序抽样并写出程序
