MonteCarlo模拟误差分析课程设计.docx

1、MonteCarlo模拟误差分析课程设计Monte Carlo模拟误差分析课程设计1.实验目的1.1学习并掌握MATLAB软件的基本功能和使用。1.2学习并掌握基于Monte Carlo Method(MCM)分析的不确定度计算方法。1.3 研究 Guide to the expression of Uncertainty in Measurement(GUM)法与 MCM 法的区别与联系和影响因素，自适应MCM方法，基于最短包含区间的MCM法。2.Monte Carlo模拟误差分析的实验原理在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不 确定度的合成公式获得间接测量

2、量的标准不确定度和扩展不确定度(GUM)o在有些场合 下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进 行计算。Monte Carlo(MCM)法就是较为常用的数学工具，具体原理相见相关资料。此次课程设汁中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数 的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不 确定度一一(MCM),并与理论上佔计标准不确定度的Bessel公式、极差法作一一(GUM) 比较，完成实验内容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响 MCM法的参数，自适应MCM法和基于最短包含区间的MCM法

3、。已知两项误差分量服从正态分布，标准不确定度分别为z/.=5mV, i/2=7mV,用统汁模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成 的置信概率P为99.73%的扩展不确定度)。3Monte Carlo模拟误差分析的实验内容3.1 MCM法与GUM法进行模拟误差分析和不确定度计算(1).利用MATLAB软件生成0, 1区间的均匀分布的随机数；(2).给出误差分量的随机值：利用MATLAB,山均匀分布随机数生成标准正态分布随机数卩，误差分量随机数可表示为Q =5】mV；同理得=f72ll2 =7%mV(3).求和的随机数：误差和的随机数二+；(4).重复以上步骤，得

4、误差和的随机数系列：J, =Jlr+J2/ / = 1,2,- n：(5).作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应 包含所有数据，按数值将区间等分为加组(加尽可能大)，每组间隔为4,记数各区间 的随机数的数目心，以/为底，以巴为高作第八丿1,2加)区间的矩形，最终的加组nA矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。k(6).以频率上 = 99.73%为界划定区间，该区间半宽即为测量总误差的置信概率为n99.73%的扩展不确定度。(7).合成的标准不确定度：A.总误差随机数平均值一 1 71B.各误差随机数的残差v； = 6C.按照Bess

5、el公式估计标准不确定度实验流程图:实验说明：本实验中随机数种子为31。以下为N分别为100000点和500000点两种悄况下,M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。实验程序：tic;clcar:clc;close all;bdclose all;%设定参数值%随机信号点数N,均值Mu,标准差Sigma%N=10A5;Mu=l;Sigma=2;M=N/10;x=0:l:M:x_=l:M;u 1=0.005;u2=0.007;%产生两个在(0,1)上服从均匀分布的,种子为31,每一次都相同的随机数XI和X2% rand(*state3r);Xl=rand(LN)

6、;X2=rand(LN);%按照Box-Mueller变换方法产生标准正态分布Y1和Y2%Yl=sqrt(-2*log(XI).*cos(2*pi*X2);Y2=sqrt(-2*log(X 1 ).*sin(2*pi*X2);%为做直方图先立义好X轴的坐标数据delta l=ul*Yl;delta2=u2*Y2;delta=delta 1 +delta2;d_delta=(max(delta)-min(delta)/(M-l); %d_delta 为误差分布的间距delta_n=min(delta):d_delta:max(delta); %delta_n 为误差分布序列%作图%髙斯随机信号f

7、igured),axis(O,Nrmax(5*Yl)jnax(5*Yl)plot(Yl );grid on;figure ，axis(0,Nrmax(5*Y2),max(5*Y2)J)plot(Y2);grid on;% hold on% plot(x,0,k);grid on;% plot(x J；r-);grid on;% plot(XrI/r-);grid on;% hold on%变换为任意均值和方差的正态分布%Zl=Sigma*Y 1 +Mu;%作图%高斯随机信号% subplot(2,2,2)% axis(0,Nr6,6)% plot(Zl);grid on;% hold on%

8、plot(x,Mu,k);% plot(x,Mu+Sigma,r-,);grid on;% pIot(x,Mu-Sigma/r-r);grid on;% hold on%正态分布误差1幅度直方图%figure axis(-l JQN)hist(deltal,M);grid on;%正态分布误差2幅度直方图%figure(4)axis(-l J AN) hist(delta2,M);grid on;%合成误差幅度直方图figure(5)axis(-lQN)H=hist(deIta.M);hist(de!ta.M);grid on;%画包络线%figure(6)HH=enveIope(x_,H);

9、plot(dclta_nHHb:);grid on;hold on;%i|-算直方图的而积%S=sum(d_dclta*abs( H);%计算直方图的而积%s_l表示正向直方图的每一个单元的而积 %s_2表示反向直方图的每一个单元的而积 %s_表示正反向两两对称每一对单元的而积 %area表示以中心为对称轴的累加面积 i=l:l:M/2;s_l(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_l(i)+s_2(i);area =s_ ;for ii=l:M/2-l a

10、rea(ii+l )=area(ii)+s_(ii);end%计算99.73%的直方图面积 for iii=l:M/2;area(iii);if(area(iii)-0.9973*S)=0;breakplot(delta_n(M/2-iii+l),delta_n(M/2+iii)JH(M/2-iii).H(M/2+iii)Jro,);grid on; delta_n_u=(delta_n(fioor(M/2+iii)-delta_n(floor(M/2-iii+l )/6;%理论计算标准不确宦度% delta_mean=mean(delta); delta_cancha=delta-delta

11、_niean: s=sqrt(sum(delta_cancha.A2)/(N-l);% toe;程序运行结果(1)当 N =100000, M=N/10 时:Figure 2Figure 4Figure 5Figure 6计算结果：s=O.OO86, delta_n_u= 0.0089o（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示:图 3.1 N=100000.M=N/5; s=O.OO86.delta_n_u= 0.009020 Q I 厂! | I I ! I -0 04 -G 03 -0 02 -0 01 0 0 01 0 02 0 03 0 04 0 05图

12、3.2 N=100000.M=N/2; s=0.0086.delta n u= 0.00915图 3.3 N= 1OOOOO. M=N; s=0.0086, delta_n_u= 0.0091 1 G1 1 1 1 G 1 M -0.03 -0.02 -0.01 0 0.01 0.02 0.03 D.04 0.06图 3.5 N=1 OOOOO,M=5N; s=0.0086; delta_n_u= 0.009125-0.06 -0.04 -0.03 -0.02 -0.01 0 0.01 0 02 0.03 0.04图 3.7 N=5OOOOO.M=N/2; s=O.OO86.delta n u

13、= 0.008618图 3.8 N=500000.M=N： s=0.0086.delta_n_u= 0.0086图 3.9 N=5OOOOO.M=2N; s=0.0086.delta n u= 0.0086图 3.10 N=500000.M=5N; s=0.0086.delta n u= 0.0086表1 N=100000时，N与M成不同倍数k时，直方图讣算结果与理论汁算结果的差异k=N/M105211/21/5S0.00860.00860.00860.00860.00860.0086delta n u0.00890.00900.00910.00910.00910.0091ldelta n u

14、 si0.00030.00040.00050.00050.00050.0005表2 N=500000时，N与M成不同倍数k时，直方图讣算结果与理论汁算结果的差异k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086delta n u0.00860.00860.00860.00860.00860.0086ldelta n u si000000实验需要讨论的问题：（1）N（采样点数），M（划分的区间数）与直方图的关系（平滑，Y轴的高度）。根据以上各图分析知：当N固定的悄况下，随着M值得增大直方图的平滑性变差, Y轴高度下降。其中，MvN时，Bes

15、sel公式计算的标准不确定度与99.73%直方图面积 的扩展不确定度两者之间的误差随着M的增大而逐渐减小。当N改变时，即当N增大 时可使直方图更为精细，且此时不改变直方图包络的基本形状。（2） Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间 会存在误差，这个误差与哪些因素有关（N,M,N=M）此误差的大小和M、N的相对大小值有关。当N=M时，由于对离散的误差值统 计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。当MN时, 增大M值，误差值稳定，但不能改善误差值。3.2自适应MCM法在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加

16、，直至所需 要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的 数值容差时，则认定该数值结果稳定。（1） .基于前一个实验，构建衡量Monte Carlo法和GUM法计算得到标准不确定度差值 的函数。（2） .将随机数个数N,分割区间数M分别作为该函数的自变量，定义自变量的取值范圉, 从而获得相应的函数值。（3） .分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N, M组合。 实验示例程序：tic;warning off;a,b=meshgrid(logspace( 1,6);for j= 1:max(size(a);for jj= 1 :max(size

17、(b);Result l(j,jj)=shiyan(a(j),b(jj);endendfigure(l ),surfic(a,b,Result 1);c=logspace(l,6);d=logspace(l,6);for jjj= 1:max(size(c);Result2(jjj)=shiyan(c(jjj),d(jjj);endfigure(2),pIot3(c.d,Result2);grid on;toe;其中shiyan ()子程序为：function y=shiyan(N,M)%clear;clc;close all;%bdclose all;%设定参数值%随机信号点数N,均值为1,

18、标准差ul,u2%N=10A5;Mu=l;%M=N/10;x=0:l:floor(M);x_=l:floor(M);ul=0.005;u2=0.007;%产生两个在(0,1)上服从均匀分布的,种子为39；每一次都相同的随机数XI和X2% rand(*state39);X1 =rand( 1 ,floor(N);X2=rand( 1 ,floor(N);%按照Box-Mueller变换方法产生标准正态分布Y1和Y2%Y1 =sqrt(-2*log(X l).*cos(2*pi*X2);Y2=sqrt(-2*log(X l).*sin(2*pi*X2);%为做直方图先定义好X轴的坐标数据%delt

19、al=uI*Yl;delta2=u2*Y2;delta=delta 1 +delta2;d_delta=(max(delta)-min(delta)./(floor(M)-1); %d_delta 为误差分布的间距delta_n=niin(delta):d_delta:max(delta); %delta_n 为误差分布序列%作图%高斯随机信号% figure(l),% axis(0,N,-max(5*Y1 ),max(5*Y 1)% plot(Yl);grid on;% figure(2),% axis(0,Nrmax(5*Y2),max(5*Y2)% plot(Y2);grid on;%

20、% hold on% % plot(x,0,k);grid on;% % plot(xj/r*);grid on;% % plot(x, Wgrid on;% % hold on%变换为任意均值和方差的正态分布% %Zl=Sigma*Yl+Mu;%作图% %ISJ斯随机信号% % subplot(222)% % axis(0,N,-6,6)% % plot(Zl);grid on;% % hold on% % plot(x,Mu,k);% % plot(x,Mu+Sigma,);grid on;% % plot(x,Mu-Sigma/rrid on;% % hold on% %正态分布误差1幅

21、度直方图% figure(3)%axis(-l,lAN)% hist(delta 1 ,M);grid on;% %正态分布误差2幅度直方图% figure(4)% axis(卜 1,1,0,N)% hist(delta2,M);grid on;%合成误差幅度直方图% figure(5)%axis(-lJAN)H=hist(delta,floor(M);% hist(delta,M);grid on;% %画包络线% figure(6)% HH=envelope(x-,H);% plot(delta_n,HH/b:*);grid on;% hold on;；15%ik算直方图的面积%S=sum

22、(d_delta.*abs(H);%计算直方图的面积%s_l表示正向直方图的每一个单元的面积 %s_2表示反向直方图的每一个单元的面积 %s_表示正反向两两对称每一对单元的面积 %area表示以中心为对称轴的累加面积 i=l:l:floor(M./2);s_l(i)=d_deIta.*abs(floor(H(floor(M./2+i);s_2(i)=d_delta.*abs(floor(H(floor(M./2-i+1); s_(i)=s_l(i)+s_2(i);area(l)=s_(l);for ii= 1:floor(M./2)-larea(ii+ l)=area(ii)+s_(ii);e

23、nd%计算99.73%的直方图面积for iii=l:floor(M./2);area(iii);if (area(iii)-0.9973*S)=0;break endend %plot(delta.n(M/24ii+lXdelta_n(M/2+iii)JH(M/2-iii).H(M/2+iii)；ro,);grid on;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;10，17x怕

24、5实验需要讨论的问题：如何根据三维网格曲线和三维曲线获得最佳的N, M组合。根据shiyan ()子函数知：程序返回值为y=abs(delta_n_u-s)；显然，当y=0时即可获得N, M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N, Mo3.3基于最短包含区间的MCM法如果PDF不对称，可釆用最短100p%包含区间。确定,使得人如-)畑-治，r = l,M-q,可得最短100p%包含区间(爲。(1).先确定 q(P=0.9973,M= 1 OOOOOO,q=PM= 10000)(2).重新计算包络线下的面积(不是对称的时候)(3).根据算法：= 计算 r 龙(4). f对应的区间

25、长度极为最短包含区间实验示例程序：%x为均匀分布，正态分布，反正弦分布%y=sin(x)为何种分布tic;clear;clc;close all;%设定参数值%随机信号点数N,均值1,标准差%N=10A6;M=N/10;x=0:l:M;x_=l:MJ;ul =0.005;u2=0.007;%randfstate31);Xl=rand(l,N);X2=rand(l,N);Yl=sin(Xl);Y 2=sqrt(-2*log(X 1). *cos(2*pi*X2);deltal=ul*Yl;delta2=u2*Y2;delta 二 delta 1 +delta2;d_delta=(max(delt

26、a)-min(delta)/(M-1); %d_delta 为误差分布的间距delta_n= niin(delta):d_delta: max(delta); %delta_n 为误差分布序列%画直方图figure(l)axis(-lJ,O,N)hist(Yl,M);grid on;figure(2)axis(-l,l,0,N)hist(Y2,M);grid on;figure(3);axis(卜 1,1,0,N)hist(delta,M);grid on;H=hist(delta,M);%画包络线%figure(4)HH=envelope(x_,H);plot(delta_n,HH,b:);

27、grid on;hold on;%i|-算直方图的面积%S=sum(d_delta*abs(H);%计算直方图的面积%s_l表示正向直方图的每一个单元的面积%s_2表示反向直方图的每一个单元的面积%s_表示正反向两两对称每一对单元的面积%area表示以中心为对称轴的累加面积 i=l:l:M/2;s_l(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_l(i)+s_2(i);area(l)=s_(l);for ii=l:M/2-larea(ii+ l)=area(ii)+s

28、_(ii);end%计算99.73%的直方图面积for iii=l:M/2;area(iii);if (area(iii)-0.9973*S)=0;breakendendplot(delta_n(M/2-iii+1 ),delta_n(M/2+iii),H(M/2-iii+1 ),H(M/2+iii)；ro,);grid on; delta_n_u 1 =(delta_n(floor(M/2+iii)-delta_n(floor(M/2-iii+1 )/6;%讣算最短包含区间面积P=09973; %置信概率q=P*M; % 参数qiiii=l:M;s2(iiii)=d_delta*abs(H(

29、iiii);area2(l)=s2(l);for j=l:M-larea2(j+ I)=area2(j)+s2(j);end%for jj=l:M-qarea2(jj);for r=l:M/2if area2(r+q)-area2(r)=area2(jj+q)-area2(jj); breakendendend plot(delta_n(r),delta_n(r+q),H(r),H(r+q),r*,);grid on;delta_n_u2=(delta_n(floor(r+q)-delta_n(floor(r)/6 %理论计算标准不确定度% delta_mean=mean(delta); delta_cancha=delta-delta_mean;s=sqrt