7.3 蒙特卡罗模拟蒙特卡罗 (Monte－ Carlo)模拟，又称蒙特卡罗方法、统计试验法等,
M－ C模拟是静态模拟，描述特定时间点上的系统行为。
模拟过程中不出现时间参数。
基本思想,把随机事件
（变量）的概率特征与数学分析的解联系起来,
概率特征,随机事件的概率和随机变量的数学期望等。
用试验方法确定一,蒙特卡罗法计算定积分例 7.3.1 用 M－ C 模拟求圆周率 π的估计值 。
设二维随机变量
(X,Y)在正方形内服从均匀分布,
(X,Y)落在圆内的概率为,
1
10P｛ X2+Y2≤1}=
4
计算机上做 n次掷点试验：
产生 n 对二维随机点 (xi,yi),i＝ 1,2,…,n,
若有 k 个点落在 l／ 4圆内其中,xi 和 yi 是 RND 随机数对,
相当于第 i
个随机点落在 1/4圆内,
检查每对数是否满足,
x2i+y2i≤1
随机事件“点落入 1/4圆内”的频率为 k/n 。
根据概率论中的大数定律，事件发生的频率依概率收敛于事件发生的概率 p,即有
1}{lim
pP nk
n
得圆周率 π的估计值为
=4k/n
且当 试验次数足够大 时，其精度也随之提高 。
分析,实际上概率值为
4
110 2 dxx
恰为 1/4圆的面积频率法,利用随机变量落进指定区域内的频率来计算定积分。
平均值法,利用随机变量的平均值（数学期望 )
来计算定积分。
ba dxxfI )(
平均值法 的算法如下：
(1)产生 RND 随机数,r1,r2,…,rn；
（ 2） 令 ui=a＋ (b－ a)ri,i=1,2,…,n；
（ 3） 计算 作为 I的估计值 。
n
i
iufn
ab
1
)(
原理分析：
设随机变量 ζ1,ζ2,…,ζn相互独立,且
ζi～ U(0,1)。
｛ f（ ξi）｝,i=1,2,…,n 相互独立同分布。
ab
Idxxf
ab
b
a?
)(1
E[f(ξi)]=
由大数定律知
n
i
i
n ab
I
f
n
1
)(1lim?
a.e.
即以概率为 1 成立。
当 n 足够大时，得近似公式：
n
i
ifnab
b
a dxxf
1
)(1)()(?
I=
注,本质上平均值法是用样本平均值作为总体教学期望的估计。
二,蒙特卡罗模拟试验次数的确定
M－ C 模拟是一种试验近似方法,试验次数如何确定？？
希望：模拟次数较少、
模拟精度较高频率法的讨论频率法是用事件 A出现的频率作为概率 p 的估计,
n
kp n
问题,试验次数 n 多大时，对给定的置信度
1－ α（ 0<α<1）,估计精度达到 ε。
即问：取多大的 n 使
1? p
n
k
PppP n
成立？
答案：
2
2
)1(
z
pp?
n=
其中,zα是正态分布的临界值。
证明平均值法 在给定 α和 ε下所需的试验次数的估计式为
0
1
2
)(
10
12
2
/
22
n
i
xix
n
S
Szn
试验次数估计式的分析
0
1
2
)(
10
12
2
/
22
n
i
xix
n
S
Szn
为估计概率 p做模拟,却又需要用 p去估计模拟次数 n。
如何计算 S2？
解决方法,先做 n0 次模拟（称为学习样本 )，根据学习样本
（ 1） 先求出 p的估计，再估计模拟次数 n:
2
2
)1(
z
pp?
n=1.
2
2
)?1(?
zpp?
n =
（ 2） 计算出的样本方差 S2,用来估计 n。
2,M － C模拟的估计精度 ε与试验次数 n的平方根成反比，若精度 ε提高 10倍，则试验次数 n
要增大 100倍 。
P107表 6.4中列出了置信度为 0.95 时，在不同精度 ε 及概率 p条件下频率法所需试验次数。
对该表进行分析，能得到什么结论？
1,精度提高，试验次数大幅提高；
2,事件发生概率越接近 0.5，试验次数越高；
核反应堆屏蔽层是用一定厚度的铅包围反应堆，用以阻挡或减弱反应堆发出的各种射线。在各种射线中，中子对人体伤害极大，因此，在屏蔽层的设计中，了解中子穿透屏蔽层的概率，对反应堆的安全运行至关重要 。 1.问题背景
D
返回穿透吸收三种状态例 7.3.2 核反应堆屏蔽层设计问题
2,简化假设：
*1 假定屏蔽层平行板厚度为 D=3d,其中 d
为两次碰撞之间中子的平均游动距离；
*2 假设在第 10次碰撞以后，中子速度下降到为某一很小数值而终止运动（被引收）。
因为每次碰撞后，中子因损失一部分能量而速度下降。
*3 假定中子在屏蔽层内相继两次碰撞之间游动的距离服从指数分布 ；
*4 中子经碰撞后的弹射角 θ～ U(0,2π).
阐述中子的运动，
为模拟做理论准备
3,中子运动的数学描述引进变量,
弹射角 θi — 第 i 次碰撞后中子的运动方向与 x
轴正向的夹角。
xi — 第 i 次碰撞后中子所处位置与屏蔽层内壁的距离 。
D
0 x
xi
θiR i— 中子在第 i 次碰撞前后的游动距离。
三个变量均为随机变量
（ 1） xi<0,中子返回反应堆 ；
（ 2） xi＞ D,中子穿透屏蔽层；
经过第 i 次碰撞，中子在屏蔽层内的位置是
xi=xi－ 1+Ricosθi,i=1,2,…,10,
中子在屏蔽层里随机游动，第 i 次碰撞以后，
按照它的位置坐标 xi,可能有以下三种情况发生：
4,模拟过程
(1) 产生 RND随机数对 (ri,ui )；
（ 3） 0＜ xi＜ D,若 i<10,中子在屏蔽层内继续运动，
若 i=10,中子被屏蔽层吸收。
中子三状态判别准则
(2) 将 (ri,ui )代入公式计算
,2
,ln
ii
ii
u
rdR
(i=1,2,3,…,10)
(3) 将 （ Ri,θi) 代入公式
xi = xi－ 1+Ricosθi,i=1,2,…,10
(4) 判断中子是否穿透屏蔽层。
5,模拟结果分析要求穿透屏蔽层的概率数量级为 10－ 6～ 10－ 10，
按假设条件得到一次 模拟结果 如下,
第 i次中子的移动距离和弹射角。
计算出第 i 次碰撞中子与内壁的距离 xi。
中子数 (个 ) 穿透 (%) 吸收 (%) 返回 (%)
100
1000
3000
5000
30.0
26.0
26.5
26.3
28.0
23.4
21.8
22.0
42.0
50.6
51.7
51.7
中子穿透屏蔽层的百分比超过了 1/4,模拟结果表明屏蔽层厚度 D=3d不合适 。
多厚的屏蔽层才能使穿透的概率
W＜ 10－ 6？
问 题,
如何解决这个问题？ 思路？
1,计算机收索法增大屏蔽层的厚度，如 D=6d,12d,24d、
36d,…,交由计算机进行模拟，并搜索到所求解。
2,分析法
1 2 m
……
D D D
设计屏蔽层的厚度,x=mD
将屏蔽层视为 m
层厚度均为 D的平行板 。
由于碰撞的能量损失，中子穿过屏蔽层的平均速度会逐层下降。
设 WD 是中子穿过厚度为 D屏蔽层的概率，则穿过整个屏蔽层的概率 W为
mDW
W＜
利用模拟结果,当 D=3d,WD≈0.25,令
(WD)m＜ 10－ 6,或 (m)4＞ 10 6
9 6 7.9
6 0 2.0
6
4l o g
10l o g 6
m
取屏蔽层的厚度 x=10D=30d，
可使穿透屏蔽层的概率 w＜ 10－ 6
注,模拟 5000个中子 的运动,用穿透屏蔽的频率估计穿透概率,由表 8.2可知精度大约只有 1％,
模拟精度太低,应适当增大模拟次数 。
总结,模拟的意义？
1.模拟方法本质上是试验性的，模拟系统是现实系统的仿真。
例中，每模拟一次相当于对一个中子的运动做一次,试验,或,观察,。
2,是对思维结果的一种验证。
3,模拟本质上是一种求解问题的试验方法，需要进行较多次数的重复模拟，并且对试验结果还需进行统计分析。
4,上千万次的模拟计算工作可以借助计算机完成，
而且运算速度是非常快。
M－ C模拟是静态模拟，描述特定时间点上的系统行为。
模拟过程中不出现时间参数。
基本思想,把随机事件
（变量）的概率特征与数学分析的解联系起来,
概率特征,随机事件的概率和随机变量的数学期望等。
用试验方法确定一,蒙特卡罗法计算定积分例 7.3.1 用 M－ C 模拟求圆周率 π的估计值 。
设二维随机变量
(X,Y)在正方形内服从均匀分布,
(X,Y)落在圆内的概率为,
1
10P｛ X2+Y2≤1}=
4
计算机上做 n次掷点试验：
产生 n 对二维随机点 (xi,yi),i＝ 1,2,…,n,
若有 k 个点落在 l／ 4圆内其中,xi 和 yi 是 RND 随机数对,
相当于第 i
个随机点落在 1/4圆内,
检查每对数是否满足,
x2i+y2i≤1
随机事件“点落入 1/4圆内”的频率为 k/n 。
根据概率论中的大数定律，事件发生的频率依概率收敛于事件发生的概率 p,即有
1}{lim
pP nk
n
得圆周率 π的估计值为
=4k/n
且当 试验次数足够大 时，其精度也随之提高 。
分析,实际上概率值为
4
110 2 dxx
恰为 1/4圆的面积频率法,利用随机变量落进指定区域内的频率来计算定积分。
平均值法,利用随机变量的平均值（数学期望 )
来计算定积分。
ba dxxfI )(
平均值法 的算法如下：
(1)产生 RND 随机数,r1,r2,…,rn；
（ 2） 令 ui=a＋ (b－ a)ri,i=1,2,…,n；
（ 3） 计算 作为 I的估计值 。
n
i
iufn
ab
1
)(
原理分析：
设随机变量 ζ1,ζ2,…,ζn相互独立,且
ζi～ U(0,1)。
｛ f（ ξi）｝,i=1,2,…,n 相互独立同分布。
ab
Idxxf
ab
b
a?
)(1
E[f(ξi)]=
由大数定律知
n
i
i
n ab
I
f
n
1
)(1lim?
a.e.
即以概率为 1 成立。
当 n 足够大时，得近似公式：
n
i
ifnab
b
a dxxf
1
)(1)()(?
I=
注,本质上平均值法是用样本平均值作为总体教学期望的估计。
二,蒙特卡罗模拟试验次数的确定
M－ C 模拟是一种试验近似方法,试验次数如何确定？？
希望：模拟次数较少、
模拟精度较高频率法的讨论频率法是用事件 A出现的频率作为概率 p 的估计,
n
kp n
问题,试验次数 n 多大时，对给定的置信度
1－ α（ 0<α<1）,估计精度达到 ε。
即问：取多大的 n 使
1? p
n
k
PppP n
成立？
答案：
2
2
)1(
z
pp?
n=
其中,zα是正态分布的临界值。
证明平均值法 在给定 α和 ε下所需的试验次数的估计式为
0
1
2
)(
10
12
2
/
22
n
i
xix
n
S
Szn
试验次数估计式的分析
0
1
2
)(
10
12
2
/
22
n
i
xix
n
S
Szn
为估计概率 p做模拟,却又需要用 p去估计模拟次数 n。
如何计算 S2？
解决方法,先做 n0 次模拟（称为学习样本 )，根据学习样本
（ 1） 先求出 p的估计，再估计模拟次数 n:
2
2
)1(
z
pp?
n=1.
2
2
)?1(?
zpp?
n =
（ 2） 计算出的样本方差 S2,用来估计 n。
2,M － C模拟的估计精度 ε与试验次数 n的平方根成反比，若精度 ε提高 10倍，则试验次数 n
要增大 100倍 。
P107表 6.4中列出了置信度为 0.95 时，在不同精度 ε 及概率 p条件下频率法所需试验次数。
对该表进行分析，能得到什么结论？
1,精度提高，试验次数大幅提高；
2,事件发生概率越接近 0.5，试验次数越高；
核反应堆屏蔽层是用一定厚度的铅包围反应堆，用以阻挡或减弱反应堆发出的各种射线。在各种射线中，中子对人体伤害极大，因此，在屏蔽层的设计中，了解中子穿透屏蔽层的概率，对反应堆的安全运行至关重要 。 1.问题背景
D
返回穿透吸收三种状态例 7.3.2 核反应堆屏蔽层设计问题
2,简化假设：
*1 假定屏蔽层平行板厚度为 D=3d,其中 d
为两次碰撞之间中子的平均游动距离；
*2 假设在第 10次碰撞以后，中子速度下降到为某一很小数值而终止运动（被引收）。
因为每次碰撞后，中子因损失一部分能量而速度下降。
*3 假定中子在屏蔽层内相继两次碰撞之间游动的距离服从指数分布 ；
*4 中子经碰撞后的弹射角 θ～ U(0,2π).
阐述中子的运动，
为模拟做理论准备
3,中子运动的数学描述引进变量,
弹射角 θi — 第 i 次碰撞后中子的运动方向与 x
轴正向的夹角。
xi — 第 i 次碰撞后中子所处位置与屏蔽层内壁的距离 。
D
0 x
xi
θiR i— 中子在第 i 次碰撞前后的游动距离。
三个变量均为随机变量
（ 1） xi<0,中子返回反应堆 ；
（ 2） xi＞ D,中子穿透屏蔽层；
经过第 i 次碰撞，中子在屏蔽层内的位置是
xi=xi－ 1+Ricosθi,i=1,2,…,10,
中子在屏蔽层里随机游动，第 i 次碰撞以后，
按照它的位置坐标 xi,可能有以下三种情况发生：
4,模拟过程
(1) 产生 RND随机数对 (ri,ui )；
（ 3） 0＜ xi＜ D,若 i<10,中子在屏蔽层内继续运动，
若 i=10,中子被屏蔽层吸收。
中子三状态判别准则
(2) 将 (ri,ui )代入公式计算
,2
,ln
ii
ii
u
rdR
(i=1,2,3,…,10)
(3) 将 （ Ri,θi) 代入公式
xi = xi－ 1+Ricosθi,i=1,2,…,10
(4) 判断中子是否穿透屏蔽层。
5,模拟结果分析要求穿透屏蔽层的概率数量级为 10－ 6～ 10－ 10，
按假设条件得到一次 模拟结果 如下,
第 i次中子的移动距离和弹射角。
计算出第 i 次碰撞中子与内壁的距离 xi。
中子数 (个 ) 穿透 (%) 吸收 (%) 返回 (%)
100
1000
3000
5000
30.0
26.0
26.5
26.3
28.0
23.4
21.8
22.0
42.0
50.6
51.7
51.7
中子穿透屏蔽层的百分比超过了 1/4,模拟结果表明屏蔽层厚度 D=3d不合适 。
多厚的屏蔽层才能使穿透的概率
W＜ 10－ 6？
问 题,
如何解决这个问题？ 思路？
1,计算机收索法增大屏蔽层的厚度，如 D=6d,12d,24d、
36d,…,交由计算机进行模拟，并搜索到所求解。
2,分析法
1 2 m
……
D D D
设计屏蔽层的厚度,x=mD
将屏蔽层视为 m
层厚度均为 D的平行板 。
由于碰撞的能量损失，中子穿过屏蔽层的平均速度会逐层下降。
设 WD 是中子穿过厚度为 D屏蔽层的概率，则穿过整个屏蔽层的概率 W为
mDW
W＜
利用模拟结果,当 D=3d,WD≈0.25,令
(WD)m＜ 10－ 6,或 (m)4＞ 10 6
9 6 7.9
6 0 2.0
6
4l o g
10l o g 6
m
取屏蔽层的厚度 x=10D=30d，
可使穿透屏蔽层的概率 w＜ 10－ 6
注,模拟 5000个中子 的运动,用穿透屏蔽的频率估计穿透概率,由表 8.2可知精度大约只有 1％,
模拟精度太低,应适当增大模拟次数 。
总结,模拟的意义？
1.模拟方法本质上是试验性的，模拟系统是现实系统的仿真。
例中，每模拟一次相当于对一个中子的运动做一次,试验,或,观察,。
2,是对思维结果的一种验证。
3,模拟本质上是一种求解问题的试验方法，需要进行较多次数的重复模拟，并且对试验结果还需进行统计分析。
4,上千万次的模拟计算工作可以借助计算机完成，
而且运算速度是非常快。
