Monte_Carlo 模拟误差分析课程设计

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

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

3、法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率P为99.73%的扩展不确定度)。3. MonteCarlo模拟误差分析的实验内容3.1MCM法与GUM法进行模拟误差分析和不确定度计算(1) .利用MATLAB软件生成0,1区间的均匀分布的随机数g；(2) .给出误差分量的随机值：利用MATLAB，由均匀分布随机数g生成标准正态分布随机数耳，误差分量随机数可表示为8=qu=5耳mV；1111同理得8=qu=7耳mV2222(3) .求和的随机数：误差和的随机数8=8+8；12(4) .重复以上步骤，得误差和的随机数系列：8=8+8i=1,2,n；i1i2i(5) .作

4、误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m组(m尽可能大)，每组间隔为A，记数各区间n的随机数的数目n.，以A为底，以j为高作第j(j=1,2m)区间的矩形，最终的m组jnA矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。丈n(6) .以频率4-=99.73%为界划定区间，该区间半宽即为测量总误差的置信概率为n99.73%的扩展不确定度。(7) .合成的标准不确定度：A.总误差随机数平均值nii=1B.各误差随机数的残差v=88iiC.按照Bessel公式估计标准不确定度i实验流程图：图6实验说明：本实验中随机数种子

5、为31。以下为N分别为100000点和500000点两种情况下M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。实验程序：tic;clear;clc;closeall;bdcloseall;%设定参数值%随机信号点数N,均值Mu,标准差Sigma%N=10人5;Mu=1;Sigma=2;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%产生两个在(0,1)上服从均匀分布的,种子为31,每一次都相同的随机数XI和X2%rand('state',31);X1=rand(1,N);X2=rand(1,N);%按照Box-Muell

6、er变换方法产生标准正态分布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

7、(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

8、9;);%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=envelo

9、pe(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);fo

10、rii=1:M/2-1area(ii+1)=area(ii)+s_(ii);end%计算99.73%的直方图面积foriii=1: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),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);d

11、elta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.A2)/(N-l);%toc;程序运行结果(1)当N=1OOOOO,M=N/1O时:FigurelFigure2Figure3Figure5Figure6计算结果：s=O.OO86,delta_n_u=0.0089。(2)当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示:35302520151017-0.04-0.03-0.02-0.0100.010.020.030.040.05图3.1N=100000,M=N/5;s=0.0086,delta_n_u=0.0090图3.2

12、N=100000,M=N/2;s=0.0086,deltanu=0.0091图3.3N=100000,M=N;s=0.0086,delta_n_u=0.0091图3.4N=100000,M=2N;s=0.0086,delta_n_u=0.0091图3.5N=100000,M=5N;s=0.0086;deltanu=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=

13、500000,M=2N;s=0.0086,deltanu=0.008604oLjKOMLT!fWWWVW5n;BI:BI:BI;s!;j888888883SSWSRSSSSyi3Ioo.G02o01o02-03Io-0.140图3.10N=500000,M=5N;s=0.0086,delta_n_u=0.0086表1N=100000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086deltanu0.00890.00900.00910.00910.00910.0091ldeltan

14、u-si0.00030.00040.00050.00050.00050.0005表2N=500000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086deltanu0.00860.00860.00860.00860.00860.0086|deltanu-s|000000实验需要讨论的问题：(1)N(采样点数)，M(划分的区间数)与直方图的关系(平滑，Y轴的高度)。根据以上各图分析知：当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。其中，M<N时，Besse

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

16、不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。(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);Result

17、1(j,jj)=shiyan(a(j),b(jj);endendfigure(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);endfigure(2),plot3(c,d,Result2);gridon;toc;其中shiyan()子程序为：functiony=shiyan(N,M)%clear;clc;closeall;%bdcloseall;%设定参数值%随机信号点数N,均值为1,标准差ul,u2%N=10人5;Mu=1;

18、%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*Y

19、2;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;%

20、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

21、(-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表示

22、正向直方图的每一个单元的面积%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);ar

23、ea(iii);if(area(iii)-0.9973*S)>=0;breakendend%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)

24、-l);y=abs(delta_n_u-s);程序运行结果：l23实验需要讨论的问题：如何根据三维网格曲线和三维曲线获得最佳的N，M组合。根据shiyan()子函数知：程序返回值为y=abs(delta_n_u-s)；显然，当y=0时即可获得N,M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N,M。使得3.3基于最短包含区间的MCM法如果PDF不对称，可采用最短100p%包含区间。确定r*，yy5yy，r=1,Mq，可得最短100p%包含区间y,y(r*+q)(r*)(r+Q)(r)L(”*)(+q)-(1) .先确定q(P=0.9973,M=1000000,q=PM=10000)(2

25、) .重新计算包络线下的面积(不是对称的时候)(3) .根据算法：yy<yy，r=1,Mq，计算r*(r*+q)(r*)(r+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*l

26、og(X1).*cos(2*pi*X2);delta1=u1*Y1;%d_delta为误差分布的间距%delta_n为误差分布序列delta2=u2*Y2;delta=delta1+delta2;d_delta=(max(delta)-min(delta)/(M-1);delta_n=min(delta):d_delta:max(delta);%画直方图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);gr

27、idon;H=hist(delta,M);%画包络线%figure(4)HH=envelope(x_,H);plot(delta_n,HH,'b:');gridon;holdon;%计算直方图的面积%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+l);s_(i)=s_l(i)+s_2(i);area(l)=s_(l);forii=l:M/2-larea(ii+l)=area(ii)+s_(ii);end%计算99.73%的直方图面积foriii=l:M/