蒙特卡罗MonteCarlo模拟误差分析课程设计哈工大误差原理课程设计.docx
《蒙特卡罗MonteCarlo模拟误差分析课程设计哈工大误差原理课程设计.docx》由会员分享，可在线阅读，更多相关《蒙特卡罗MonteCarlo模拟误差分析课程设计哈工大误差原理课程设计.docx（40页珍藏版）》请在冰点文库上搜索。
蒙特卡罗MonteCarlo模拟误差分析课程设计哈工大误差原理课程设计
MonteCarlo模拟误差分析课程设计
MonteCarlo模拟误差分析课程设计
1.实验目的
1.1学习并掌握MATLAB软件的基本功能和使用。
1.2学习并掌握基于MonteCarloMethod（MCM）分析的不确定度计算方法。
1.3研究GuidetotheexpressionofUncertaintyinMeasurement（GUM）法与MCM法的区别与联系和影响因素，自适应MCM方法，基于最短包含区间的MCM法。
2.MATLAB软件介绍实验内容
2.1介绍MATLAB软件的基本知识
MATLAB名字由MATrix和LABoratory两词的前三个字母组合而成。
20世纪七十年代，时任美国新墨西哥大学计算机科学系主任的CleveMoller出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK矩阵软件工具包库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB
MATLAB语言的主要特点
（1）.具有丰富的数学功能
（2）.具有很好的图视系统
（3）.可以直接处理声言和图形文件。
（4）.具有若干功能强大的应用工具箱。
（5）.使用方便，具有很好的扩张功能。
（6）.具有很好的帮助功能
演示内容：
（1）.MATLAB的数值计算功能
在“命令行”Command提示窗口中键入：
“A=eye（5,5）;A=zeros（5,5）;A=ones（5,5）”等命令生成各类矩阵；在“命令行”Command提示窗口中键入：
“[v,d]=eig（A）”生成特征矩阵和特征向量；在“命令行”Command提示窗口中键入：
“expm（A）”对矩阵A求幂；在“命令行”Command提示窗口中键入：
x=[135];y=[246];z=conv（x,y）；显示结果：
z=210283830
（2）.MATLAB的符号计算功能
在“命令行”Command提示窗口中键入：
symsax；f=sin（a*x）;df=diff（f,x）;dfa=diff（f,a）；
Command提示窗口显示结果：
df=cos（a*x）*a；dfa=cos（a*x）*x；
2.2MATLAB软件画图特性
（1）.MATLAB二维绘图
命令函数：
plot
参数：
线型、颜色、多重线、网格和标记、画面窗口分割、其他方式、隐函数的描绘）
（2）.MATLAB三维画图
曲面与网格图命令函数：
mesh
三维带阴影曲面图：
surf
三维曲线命令：
plot3
演示内容：
（1）.MATLAB的二维绘图功能
在命令行Command提示窗口中键入：
closeall;x=linspace（0,2*pi,100）;%100个点的x座标
y=sin（x）;%对应的y座标
plot（x,y）;
得到如下的结果：
图1
在命令行Command提示窗口中键入：
“plot（x,sin（x）,x,cos（x））;”
得到如下的结果：
图2
在命令行Command提示窗口中键入：
plot（x,sin（x）,'co',x,cos（x）,'g*'）;
得到如下的结果：
图3
在命令行Command提示窗口中键入：
xlabel（'InputValue'）;%x轴注解
ylabel（'FunctionValue'）;%y轴注解
title（'TwoTrigonometricFunctions'）;%图形标题
legend（'y=sin（x）','y=cos（x）'）;%图形注解
gridon;%显示格线
得到如下的结果：
图4
（2）.MATLAB的多维绘图功能
在命令行Command提示窗口中键入：
[X,Y]=meshgrid（-3:
0.125:
3）;%生成二维网格点
Z=peaks（X,Y）;%生成某种内置函数
mesh（X,Y,Z）;
得到如下的结果
图5
其他的演示功能详见“MATLAB画图文档”
3.MonteCarlo模拟误差分析的实验原理
在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度（GUM）。
在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。
MonteCarlo（MCM）法就是较为常用的数学工具，具体原理相见相关资料。
此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度——（MCM），并与理论上估计标准不确定度的Bessel公式、极差法作——（GUM）比较，完成实验内容。
并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。
已知两项误差分量服从正态分布，标准不确定度分别为
mV，
mV，用统计模拟分析法给出两项误差和的分布（误差分布的统计直方图，合成的标准差，合成的置信概率P为99.73%的扩展不确定度）。
4.MonteCarlo模拟误差分析的实验内容
4.1MCM法与GUM法进行模拟误差分析和不确定度计算
（1）.利用MATLAB软件生成[0，1]区间的均匀分布的随机数
；
（2）.给出误差分量的随机值：
利用MATLAB，由均匀分布随机数
生成标准正态分布随机数
，误差分量随机数可表示为
mV；
同理得
mV
（3）.求和的随机数：
误差和的随机数
；
（4）.重复以上步骤，得误差和的随机数系列：
；
（5）.作误差和的统计直方图：
以误差数值为横坐标，以频率为纵坐标作图。
作图区间应包含所有数据，按数值将区间等分为
组（
尽可能大），每组间隔为
，记数各区间的随机数的数目
，以
为底，以
为高作第
（
）区间的矩形，最终的
组矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。
（6）.以频率
为界划定区间，该区间半宽即为测量总误差的置信概率为99.73%的扩展不确定度。
（7）.合成的标准不确定度：
A.总误差随机数平均值
B.各误差随机数的残差
C.按照Bessel公式估计标准不确定度
实验流程图：
图6
实验1
本实验中随机数种子为28。
以下为N分别为100000点和500000点两种情况下，M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。
1）程序
tic;
clear;clc;closeall;
%%设定参数值%%
%%随机信号点数N，均值为1，标准差u1,u2%%
N=100000;
M=N/10;
x=0:
1:
M;
x_=[1:
M];
u1=0.005;
u2=0.007;
%%产生两个在（0,1）上服从均匀分布的,种子为28,每一次都相同的随机数X1和X2%%
rand（'state',28）;
X1=rand（1,N）;
X2=rand（1,N）;
%%按照Box-Mueller变换方法产生标准正态分布Y1和Y2%%
Y1=sqrt（-2*log（X1））.*cos（2*pi*X2）;
Y2=sqrt（-2*log（X1））.*sin（2*pi*X2）;
%%为做直方图先定义好X轴的坐标数据%%
delta1=u1*Y1;
delta2=u2*Y2;
delta=delta1+delta2;
d_delta=（max（delta）-min（delta））/（M-1）;%%d_delta为误差分布的间距
delta_n=[min（delta）:
d_delta:
max（delta）];%%delta_n为误差分布序列
%%作图%%
%%高斯随机信号%%
figure
（1）,
axis（[0,N,-max（5*Y1）,max（5*Y1）]）
plot（Y1）;gridon;title（'学号：
13S101028姓名：
李明'）;
figure
（2）,
axis（[0,N,-max（5*Y2）,max（5*Y2）]）
plot（Y2）;gridon;title（'学号：
13S101028姓名：
李明'）;
%holdon
%plot（x,0,'k'）;gridon;
%plot（x,1,'r--'）;gridon;
%plot（x,-1,'r--'）;gridon;
%holdon
%%变换为任意均值和方差的正态分布%%
%Z1=Sigma*Y1+Mu;
%%作图%%
%%高斯随机信号%%
%subplot（2,2,2）
%axis（[0,N,-6,6]）
%plot（Z1）;gridon;
%holdon
%plot（x,Mu,'k'）;
%plot（x,Mu+Sigma,'r--'）;gridon;
%plot（x,Mu-Sigma,'r--'）;gridon;
%holdon
%%正态分布误差1幅度直方图%%
figure（3）
axis（[-1,1,0,N]）
hist（delta1,M）;gridon;title（'学号：
13S101028姓名：
李明'）;
%%正态分布误差2幅度直方图%%
figure（4）
axis（[-1,1,0,N]）
hist（delta2,M）;gridon;title（'学号：
13S101028姓名：
李明'）;
%%合成误差幅度直方图%%
figure（5）
axis（[-1,1,0,N]）
H=hist（delta,M）;
hist（delta,M）;gridon;title（'学号：
13S101028姓名：
李明'）;
%%画包络线%%
figure（6）
HH=envelope（x_,H）;
plot（delta_n,HH,'b:
'）;gridon;title（'学号：
13S101028姓名：
李明'）;
holdon;
%%计算直方图的面积%%
S=sum（d_delta*abs（H））;
%%计算直方图的面积%%
%%s_1表示正向直方图的每一个单元的面积
%%s_2表示反向直方图的每一个单元的面积
%%s_表示正反向两两对称每一对单元的面积
%%area表示以中心为对称轴的累加面积
i=1:
1:
M/2;
s_1（i）=d_delta*abs（floor（H（floor（M/2+i））））;
s_2（i）=d_delta*abs（floor（H（floor（M/2-i+1））））;
s_（i）=s_1（i）+s_2（i）;
area
（1）=s_
（1）;
forii=1:
M/2-1
area（ii+1）=area（ii）+s_（ii）;
end
%%计算99.73%的直方图面积
foriii=1:
M/2;
area（iii）;
if（area（iii）-0.9973*S）>=0;
break
end
end
plot（[delta_n（M/2-iii+1）,delta_n（M/2+iii）],[H（M/2-iii+1）,H（M/2+iii）],'ro'）;gridon;title（'学号：
13S101028姓名：
李明'）;
delta_n_u=（delta_n（floor（M/2+iii））-delta_n（floor（M/2-iii+1）））/6
%%理论计算标准不确定度%%
delta_mean=mean（delta）;
delta_cancha=delta-delta_mean;
s=sqrt（（sum（delta_cancha.^2））/（N-1））
toc;
2）程序运行结果图
（1）当N=100000，M=N/10时：
图1
图2
图3
图4
图5
图6
计算结果：
s=0.0086，delta_n_u=0.0085。
（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示：
图7N=100000,M=N/5;s=0.0086,delta_n_u=0.0086
图8N=100000,M=N/2;s=0.0086,delta_n_u=0.0086
图9N=100000,M=N;s=0.0086,delta_n_u=0.0086
图10N=100000,M=2N;s=0.0086,delta_n_u=0.0086
图11N=100000,M=5N;s=0.0086;delta_n_u=0.0086
（3）当N=500000时，计算结果如下所示：
图12
图13
图14
图15
图16
图17
计算结果：
s=0.0086，delta_n_u=0.0088。
（4）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示：
图18N=500000,M=N/5;s=0.0086,delta_n_u=0.0088
图19N=500000,M=N/2;s=0.0086,delta_n_u=0.0088
图20N=500000,M=N;s=0.0086,delta_n_u=0.0088
图21N=500000,M=2N;s=0.0086,delta_n_u=0.0088
图22N=500000,M=5N;s=0.0086,delta_n_u=0.0088
表1N=100000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异
k=N/M
10
5
2
1
1/2
1/5
s
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
delta_n_u
0.0085
0.0086
0.0086
0.0086
0.0086
0.0086
|delta_n_u﹣s|
0.0001
0
0
0
0
0
表2N=500000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异
k=N/M
10
5
2
1
1/2
1/5
s
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
delta_n_u
0.0088
0.0088
0.0088
0.0088
0.0088
0.0088
|delta_n_u﹣s|
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
3）实验需要讨论的问题
（1）N（采样点数），M（划分的区间数）与直方图的关系（平滑，Y轴的高度）。
根据以上各图分析知：
当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。
其中，M当N改变时，即当N增大时可使直方图更为精细，且此时不改变直方图包络的基本形状。
（2）Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间会存在误差，这个误差与哪些因素有关（N,M,N>=M）
此误差的大小和M、N的相对大小值有关。
当N>=M时，由于对离散的误差值统计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。
当M>N时，增大M值，误差值稳定，但不能改善误差值。
实验2——自适应MCM法
在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。
如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。
（1）.基于前一个实验，构建衡量MonteCarlo法和GUM法计算得到标准不确定度差值的函数。
（2）.将随机数个数N，分割区间数M分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。
（3）.分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N，M组合。
1）程序
tic;
warningoff;
[a,b]=meshgrid（logspace（1,6））;
forj=1:
max（size（a））;
forjj=1:
max（size（b））;
Result1（j,jj）=shiyan（a（j）,b（jj））;
end
end
figure
（1）,surfc（a,b,Result1）;title（'学号：
13S101028姓名：
李明'）;
c=logspace（1,6）;
d=logspace（1,6）;
forjjj=1:
max（size（c））;
Result2（jjj）=shiyan（c（jjj）,d（jjj））;
end
figure
（2）,plot3（c,d,Result2）;gridon;title（'学号：
13S101028姓名：
李明'）;
toc;
其中shiyan（）子程序为：
functiony=shiyan（N,M）
%%clear;clc;closeall;
%%bdcloseall;
%%设定参数值%%
%%随机信号点数N，均值为1，标准差u1,u2%%
%N=10^5;
Mu=1;
%M=N/10;
x=0:
1:
floor（M）;
x_=[1:
floor（M）];
u1=0.005;
u2=0.007;
%%产生两个在（0,1）上服从均匀分布的,种子为28,每一次都相同的随机数X1和X2%%
rand（'state',28）;
X1=rand（1,floor（N））;
X2=rand（1,floor（N））;
%%按照Box-Mueller变换方法产生标准正态分布Y1和Y2%%
Y1=sqrt（-2*log（X1））.*cos（2*pi*X2）;
Y2=sqrt（-2*log（X1））.*sin（2*pi*X2）;
%%为做直方图先定义好X轴的坐标数据%%
delta1=u1*Y1;
delta2=u2*Y2;
delta=delta1+delta2;
d_delta=（max（delta）-min（delta））./（floor（M）-1）;%%d_delta为误差分布的间距
delta_n=[min（delta）:
d_delta:
max（delta）];%%delta_n为误差分布序列
%%作图%%
%%高斯随机信号%%
%figure
（1）,
%axis（[0,N,-max（5*Y1）,max（5*Y1）]）
%plot（Y1）;gridon;
%
%figure
（2）,
%axis（[0,N,-max（5*Y2）,max（5*Y2）]）
%plot（Y2）;gridon;
%%holdon
%%plot（x,0,'k'）;gridon;
%%plot（x,1,'r--'）;gridon;
%%plot（x,-1,'r--'）;gridon;
%%holdon
%
%%%变换为任意均值和方差的正态分布%%
%%Z1=Sigma*Y1+Mu;
%
%%%作图%%
%%%高斯随机信号%%
%%subplot（2,2,2）
%%axis（[0,N,-6,6]）
%%plot（Z1）;gridon;
%%holdon
%%plot（x,Mu,'k'）;
%%plot（x,Mu+Sigma,'r--'）;gridon;
%%plot（x,Mu-Sigma,'r--'）;gridon;
%%holdon
%%%正态分布误差1幅度直方图%%
%figure（3）
%axis（[-1,1,0,N]）
%hist（delta1,M）;gridon;
%%%正态分布误差2幅度直方图%%
%figure（4）
%axis（[-1,1,0,N]）
%hist（delta2,M）;gridon;
%%%合成误差幅度直方图%%
%figure（5）
%axis（[-1,1,0,N]）
H=hist（delta,floor（M））;
%hist（delta,M）;gridon;
%%%画包络线%%
%figure（6）
%HH=envelope（x_,H）;
%plot（delta_n,HH,'b:
'）;gridon;
%holdon;
%%计算直方图的面积%%
S=sum（d_delta.*abs（H））;
%%计算直方图的面积%%
%%s_1表示正向直方图的每一个单元的面积
%%s_2表示反向直方图的每一个单元的面积
%%s_表示正反向两两对称每一对单元的面积
%%area表示以中心为对称轴的累加面积
i=1:
1:
floor（M./2）;
s_1（i）=d_delta.*abs（floor（H（floor（M./2+i））））;
s_2（i）=d_delta.*abs（floor（H（floor（M./2-i+1））））;
s_（i）=s_1（i）+s_2（i）;
area
（1）=s_
（1）;
forii=1:
floor（M./2）-1
area（ii+1）=area（ii）+s_（ii）;
end
%%计算99.73%的直方图面积
foriii=1:
floor（M./2）;
area（iii）;
if（area（iii）-0.9973*S）>=0;
break
end
end
%plot（[delta_n（M/2-iii+1）,delta_n（M/2+iii）],[H（M/2-iii）,H（M/2+iii）],'ro'）;gridon;
delta_n_u=（delta_n（floor（M./2+iii））-delta_n（floor（M./2-iii+1）））./6;
%%理论计算标准不确定度%%
delta_mean=mean（delta）;
delta_cancha=delta-delta_mean;
s=sqrt（（sum（delta_cancha.^2））./（floor（N）-1））;
y=abs（delta_n_u-s）;
2）程序运行结果
图23
图24
3）实验需要讨论的问题：
如何根据三维网格曲线和三维曲线获得最佳的N，M组合。
根据shiyan（）子函数知：
程序返回值为y=abs（delta_n_u-s）；显然，当y=0时即可获得N，M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N，M。
实验3——基于最短包含区间的MCM法
如果PDF不对称，可采用最短
包含区间。
确定
，使得
，
，可得最短
升级会员 