第4章 蒙特卡罗仿真方法

m 事件概率的估计为： p N
其中m为N次试验中事件出现的次数。
2016/6/21
航空航天中的计算方法
Page 18
4.4 随机模拟结果的统计分析 作业： 考虑炮弹在二维铅锤平面内的运动，假设地面为一水平 面，重力加速度为9.8m/s2，大气密度为1.29kg/m3。 炮弹 质量30kg，阻力系数为0.1，迎风截面积为0.05m2。初速 度为[95,105]m/s区间的均匀分布，发射角（初速度方向 与水平面的夹角）服从均值为40deg、均方差为0.1deg的 正态分布，炮弹初始高度为零。炮弹飞行受水平方向大 小为[-5, 5]m/s 区间上均匀分布的风速影响，且假设一发 炮弹飞行过程中风速不变。为评定炮弹落点精度，请建 立炮弹运动的微分方程，并采用蒙特卡罗法统计炮弹落 点散布的均值和方差。
4.3 随机变量的抽样
X 1 Y1 X 2 J Y1 X n Y1 X 1 Y2 X 2 Y2 X n Y2 X 1 Yn X 2 Yn 0 X n Yn
则随机变量 Y1 , Y2 , , Yn 的联合概率密度函数为：
1 f g1 y1 , y2 , , yn , g21 y1 , y2 , , yn , , gn1 y1 , y2 , , yn
α tα
0.5 0.6745
0.05 1.96
0.003 3
查标准正态 分布表确定
精度提高一个数量级，试验次数N需增加2个数量级
2016/6/21 航空航天中的计算方法 Page 9
4.3 随机变量的抽样 4.3 随机变量的抽样 4.3.1 随机数的产生
随机数序列——[0,1]上服从均匀分布的随机变量的抽样 序列。 用专门的算法和程序可在计算机上生成伪随机数。 平方取中法、线性同余法、加同余法、平方同余法 确定的算法、给定初值（种子），确定的序列！！！ 需对伪随机数进行独立性、均匀性等统计性质的检验。 周期长、统计性质好的伪随机数，用于模拟试验造成的 误差不大。
2016/6/21
航空航天中的计算方法
Page 6
4.2 蒙特卡罗法的理论依据 4.2 蒙特卡罗法的理论依据 大数定理
如果随机变量X1，X2，…，XN独立同分布，且具有有限 期望值E(X)<∞，则： 1 N X N Xi P lim X N E ( X ) 1 N N i 1
2016/6/21
航空航天中的计算方法
Page 19
航空航天中的计算方法
授课教师：陈琪锋 中南大学航空航天学院
第4章 蒙特卡罗仿真方法
内容提要 4.1 4.2 4.3 4.4 蒙特卡罗法的基本概念 蒙特卡罗法的理论依据 随机变量的抽样 随机模拟结果的统计分析
参考资料： [1] 闫晓东,许志. 飞行器系统仿真实训教程, 西北工业大学出 版社，2013. Chp. 6. [2] David L. Darmofal, Computational Methods in Aerospace Engineering (Lecture Notes), MIT, 2005. Chp. 15, 16. [3] Averill M. Law, 仿真建模与分析 (Fourth Edition), 清华 大学出版社, 2009. Chp. 1.5.
2016/6/21 航空航天中的计算方法 Page 16
4.4 随机模拟结果的统计分析 4.4 随机模拟结果的统计分析 随机模拟试验的目的是获得所研究的随机现象的统计特 性。 把所要研究的随机现象用随机变量X表示，N次随机模拟 试验得到的结果（抽样）序列为 x1 , x2 , , x N 。 X的数学期望和方差分别用样本均值和样本方差来统计： 1 N x xi N i 1 1 N 2 2 S xi x N 1 i 1
4.3 随机变量的抽样 函数变换法 已知随机变量 X 1 , X 2 , , X n 的联合概率密度函数为 f x1 , x2 , , xn 。 随机变量的变换
Y1 g1 X 1 , X 2 , , X n Y2 g2 X 1 , X 2 , , X n Yn gn X 1 , X 2 , , X n


即随机变量X的简单子样的算术平均值，当子样数N充分 大时，以概论1收敛于它的期望值。 子样均值与数学期望的偏差？
2016/6/21
航空航天中的计算方法
Page 7
4.2 蒙特卡罗法的理论依据
E[ X N ]
方差： 2 E X N 2 2 1 E 2 X 1 X 2 2 X 1 X 2 N
连续随机变量X的分布函数为F(x)是值域在[0,1]上单调增 的连续函数，利用ri得到X的抽样值：
xi F 1 ( ri )
函数F-1(R)可显示表示时，才能应用。
2016/6/21 航空航天中的计算方法 Page 12
4.3 随机变量的抽样 例：指数分布随机变量的生成 参数为a的指数分布的概率密度函数： f ( x ) ae ax , x 0 累积分布函数：

j 1
j
依次产生随机数ri，若 P ( k 1) ri P ( k ) , k 1, 2, , n 则第i次抽样时随机变量X的抽样值为xk。
2016/6/21
航空航天中的计算方法
Page 11
4.3 随机变量的抽样 4.3.3 连续随机变量的抽样 直接法（反变换法）
[0,1]上服从均匀分布的随机变量R的分布函数： 0, r 0 G(r ) r , 0 r 1 1, r 1
2 2
பைடு நூலகம்变换：
x cos y sin
随机变量Ϙ和Θ的分布函数： F ( ) 1 e
, F ( ) 2
用反变换法得Ϙ和Θ的抽样： 2 ln r1 , 2 r2
x 2 ln r1 cos(2 r2 ) 得到随机变量X，Y的抽样： y 2 ln r1 sin(2 r2 )
2016/6/21 航空航天中的计算方法 Page 10
4.3 随机变量的抽样 4.3.2 离散随机变量的抽样 设随机变量X，其概率分布为：
P X xk pk , k 1, 2, , n
p
k 1
n
k
1
可按轮盘赌的方式，由产生的随机数序列进行抽样。 k 令： P (0) 0, P ( k ) p, k 1, 2, , n
1 X 1 g1 Y1 , Y2 , , Yn 1 X 2 g2 Y1 , Y2 , , Yn 1 X n gn Y1 , Y2 , , Yn
若逆变换存在且具有一阶连续偏导数，即
2016/6/21
航空航天中的计算方法
Page 14
2016/6/21 航空航天中的计算方法 Page 17
4.4 随机模拟结果的统计分析
x方差的估值为：
2 N S 1 2 2 x xi x N N ( N 1) i 1
数学期望的置信区间（置信水平为1- α）为：
S S x t , x t N N
2016/6/21
航空航天中的计算方法
Page 4
4.1 蒙特卡罗法的基本概念 基本思想： 对随机现象总体进行抽样，通过样本的仿真和分析， 对结果进行统计。 基本过程：
对影响系统的随机现象根据其概率分布产生随机样本 （抽样）； 用确定性方法对样本进行分析； 对随机样本分析结果进行统计，得到统计规律。


得：
E X N N
2
2
2016/6/21
航空航天中的计算方法
Page 8
4.2 蒙特卡罗法的理论依据 中心极限定理 如果随机变量X1，X2，…，XN独立同分布，且具有有限 非零的方差σ2，则当N充分大时，有： t 2 t t 2 2 P X N E( X ) e dt 1 0 N 2 其中α称为置信度，1- α称为置信水平。 X N N , N 蒙特卡罗方法（均值）的误差： N
2016/6/21 航空航天中的计算方法 Page 3
4.1 蒙特卡罗法的基本概念 4.1 蒙特卡罗法的基本概念 问题
干扰和偏差客观存在，无法消除： 风、推力误差、结构和装配误差 如何分析干扰和偏差对飞行过程的影响？ 随机因素对系统运动影响的统计规律？
蒙特卡罗（Monte Carlo）法
通过随机变量的统计试验（或随机模拟），求解问题近 似解的数值方法，又称统计试验法或随机模拟法。
F ( x ) f ( t )dt 1 e ax ,
0 x
x0
反函数：
1 F 1 ( r ) ln(1 r ), r 0,1 a
用随机数产生服从指数分布的随机变量X的抽样： 1 x ln r , R U 0,1 a
2016/6/21 航空航天中的计算方法 Page 13
2016/6/21
航空航天中的计算方法
Page 15
4.3 随机变量的抽样 例：正态分布随机变量的生成 设X，Y独立服从标准正态分布：
1 x2 y2 2 f ( x, y) e 2 1 2 2 f ( , ) e f ( ) f ( ) 2 1 2 2 f ( ) e , f ( ) 2
2016/6/21
航空航天中的计算方法
Page 5
4.1 蒙特卡罗法的基本概念 蒙特卡罗法的基本步骤 （1）建立系统数学模型；
（2）确定系统中的各种随机扰动因素及概率分布； （3）产生各随机扰动变量的抽样值； （4）将随机扰动变量的抽样值输入系统数学模型，进行 确定性仿真计算，得到扰动下系统运动过程； （5）重复（4），进行多次抽样仿真，得到系统运动过 程子样； （6）对子样仿真结果进行处理，得到系统运动的统计特 征量。