MonteCarlo模拟误差分析课程设计.docx
《MonteCarlo模拟误差分析课程设计.docx》由会员分享，可在线阅读，更多相关《MonteCarlo模拟误差分析课程设计.docx（31页珍藏版）》请在冰豆网上搜索。
MonteCarlo模拟误差分析课程设计
MonteCarlo模拟误差分析课程设计
1.实验目的
1.1学习并掌握MATLAB软件的基本功能和使用。
1.2学习并掌握基于MonteCarloMethod（MCM）分析的不确定度计算方法。
1.3研究GuidetotheexpressionofUncertaintyinMeasurement（GUM法与MCM法的
区别与联系和影响因素，自适应MCM方法，基于最短包含区间的MCM法。
2.MonteCarlo模拟误差分析的实验原理
在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不
确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度（GUM）。
在有些场合
下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。
MonteCarlo（MCM）法就是较为常用的数学工具，具体原理相见相关资料。
此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数
的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度一一（MCM），并与理论上估计标准不确定度的Bessel公式、极差法作一一（GUM）比较，完成实验内容。
并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。
已知两项误差分量服从正态分布，标准不确定度分别为5=5mV，U2=7mV，用
统计模拟分析法给出两项误差和的分布（误差分布的统计直方图，合成的标准差，合成的置信概率P为99.73%的扩展不确定度）。
3.MonteCarlo模拟误差分析的实验内容
3.1MCM法与GUM法进行模拟误差分析和不确定度计算
（1）.利用MATLAB软件生成［0，1］区间的均匀分布的随机数；
（2）.给出误差分量的随机值：
利用MATLAB，由均匀分布随机数打生成标准正态分布随机数1，误差分量随机
数可表示为
同理得
、2二2u2=72mV
（3）.求和的随机数：
误差和的随机数「•二:
.2；
（4）.重复以上步骤，得误差和的随机数系列：
、：
i-ii=2ii=1,2-n;
（5）.作误差和的统计直方图：
以误差数值为横坐标，以频率为纵坐标作图。
作图区间应包含所有数据，按数值将区间等分为m组（m尽可能大），每组间隔为.=，记数各区间
n
的随机数的数目nj，以厶为底，以j为高作第j（j=1,2…m）区间的矩形，最终的m组nA
矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。
k
（6）.以频率公99.73%为界划定区间，该区间半宽即为测量总误差的置信概率为
n
99.73%的扩展不确定度。
（7）.合成的标准不确定度：
A.总误差随机数平均值
B.各误差随机数的残差
Vj=「--■
C.按照Bessel公式估计标准不确定度
In
亡Vi2
u=s=彳一\n—1
实验流程图:
实验说明：
本实验中随机数种子为31。
以下为N分别为100000点和500000点两种情况下,
M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。
实验程序:
tic;
clear;clc;closeall;
bdcloseall;
%%设定参数值%%
%%随机信号点数N，均值Mu，标准差Sigma%%
N=10A5;
Mu=1;
Sigma=2;
M=N/10;
x=0:
1:
M;
x_=[1:
M];
u1=0.005;
u2=0.007;
%%产生两个在（0,1）上服从均匀分布的,种子为31,每一次都相同的随机数X1和X2%%rand（'state',31）;
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;
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,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:
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）,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.A2））/（N-1））;
%%%%%%%%%%%%%%%
toc;
程序运行结果
（1）当N=100000,M=N/10时:
012345678910
xIO4
Figure1
i-
12346676&
Figure2
Figure4
-0.04^）.03-0.02^0.0100.D1i0.020.030.040.05
Figure5
Figure6
计算结果：
s=0.0086,delta_n_u=0.0089
（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示:
图3.1N=100000,M=N/5;s=0.0086,delta_n_u=0.0090
图3.2N=100000,M=N/2;s=0.0086,deltanu=0.0091
图3.3N=100000,M=N;s=0.0086,delta_n_u=0.0091
10
IIIIX■II
04-0.03-0.020ooiW03OM005
图3.5N=100000,M=5N;s=0.0086;delta_n_u=0.0091
图3.6N=500000,M=N/5;s=0.0086,delta_n_u=0.0086
图3.7N=500000,M=N/2;s=0.0086,delta_n_u=0.0086
图3.8N=500000,M=N;s=0.0086,delta_n_u=0.0086
图3.9N=500000,M=2N;s=0.0086,deltanu=0.0086
图3.10N=500000,M=5N;s=0.0086,deltanu=0.0086
表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
deltanu
0.0089
0.0090
0.0091
0.0091
0.0091
0.0091
|deltanu-s|
0.0003
0.0004
0.0005
0.0005
0.0005
0.0005
表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
deltanu
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
|deltanu-s|
0
0
0
0
0
0
实验需要讨论的问题：
（1）N（采样点数），M（划分的区间数）与直方图的关系（平滑，丫轴的高度）。
根据以上各图分析知：
当N固定的情况下，随着M值得增大直方图的平滑性变差，丫轴高度下降。
其中，M当N改变时，即当N增大时可使直方图更为精细，且此时不改变直方图包络的基本形状。
（2）Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间会存在误差，这个误差与哪些因素有关（N,M,N>=M）
此误差的大小和M、N的相对大小值有关。
当N>=M时，由于对离散的误差值统计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。
当M>N时,增大M值，误差值稳定，但不能改善误差值。
3.2自适应MCM法
在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。
如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。
（1）.基于前一个实验，构建衡量MonteCarlo法和GUM法计算得到标准不确定度差值的函数。
（2）.将随机数个数N，分割区间数M分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。
（3）.分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N,M组合。
实验示例程序：
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
endfigure
（1）,surfc（a,b,Result1）;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;
toc;
其中shiyan（）子程序为：
functiony=shiyan（N,M）
%%clear;clc;closeall;
%%bdcloseall;
%%设定参数值%%
%%随机信号点数N，均值为1,标准差u1,u2%%
%N=10A5;
Mu=1;
%M=N/10;
x=0:
1:
floor（M）;
x_=[1:
floor（M）];
u1=0.005;
u2=0.007;
%%产生两个在（0,1）上服从均匀分布的,种子为39,每一次都相同的随机数X1和X2%%rand（'state',39）;
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）-1area（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_cancha42））./（floor（N）-1））;y=abs（delta_n_u-s）;
程序运行结果：
Figure1
00
Figure2
-
o
-
-
-
5
-CJ-
o
oo
1—
实验需要讨论的问题：
如何根据三维网格曲线和三维曲线获得最佳的N,M组合。
根据shiyan（）子函数知：
程序返回值为y=abs（delta_n_u-s；显然，当y=0时即可
获得N,M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N,M。
3.3基于最短包含区间的MCM法
如果PDF不对称，可采用最短100p%包含区间。
确定r*，使得
y（r*.q）-y（r*）空y（rq）-「=1,川，M_q，可得最短100p%包含区间口）。
（1）.先确定q（P=0.9973,M=1000000,q=PM=10000）
（2）.重新计算包络线下的面积（不是对称的时候）
（3）.根据算法：
y（r*q）-y（r*）-y（rq）-y（r），r=h川,M_q，计算r
（4）.r*对应的区间长度极为最短包含区间
实验示例程序：
%%x为均匀分布，正态分布，反正弦分布
%%y=sin（x）为何种分布
tic;
clear;clc;closeall;
%%设定参数值%%
%%随机信号点数N，均值1,标准差%%
N=10A6;
M=N/10;
x=0:
1:
M;
x_=[1:
M];
u1=0.005;
u2=0.007;
%%
rand（'state',31）;
X1=rand（1,N）;
X2=rand（1,N）;
Y1=sin（X1）;
Y2=sqrt（-2*log（X1））.*cos（2*pi*X2）;
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（[-1,1,0,N]）hist（Y1,M）;gridon;
figure
（2）axis（[-1,1,0,N]）hist（Y2,M）;gridon;
figure（3）;axis（[-1,1,0,N]）hist（delta,M）;gridon;H=hist（delta,M）;
%%画包络线%%figure（4）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:
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;delta_n_u1=（delta_n（floor（M/2+iii））-delta_n（floor（M/2-iii+1）））/6;
%%计算最短包含区间面积
P=0.9973;%置信概率q=P*M;%参数q
iiii=1:
M;
s2（iiii）=d_delta*abs（（H（iiii）））;area2
（1）=s2
（1）;
forj=1:
M-1area2（j+1）=ar
升级会员 