哈工大误差分析课程设计-Monte_Carlo

MonteCarlo模拟误差分析课程设计1实验目的1.1了解MATLAB软件的基本功能和使用1.2学习不确定度的统计模拟分析方法1.3研究误差概率密度函数和Bessel公式获得扩展不确定度的方法和影响因素2实验原理在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度(GUM)。在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。MonteCarlo(MCM)法就是较为常用的数学工具，具体原理相见相关资料。此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度——(MCM)，并与理论上估计标准不确定度的Bessel公式、极差法作——(GUM)比较，完成实验内容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。已知两项误差分量服从正态分布，标准不确定度分别为51umV，72umV，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率P为99.73%的扩展不确定度)。3实验内容(1).利用MATLAB软件生成[0，1]区间的均匀分布的随机数；(2).给出误差分量的随机值：利用MATLAB，由均匀分布随机数1生成标准正态分布随机数1，误差分量随机数可表示为11115umV；22227umV(3).求和的随机数：误差和的随机数21；(4).重复以上步骤，得误差和的随机数系列：iii21ni,2,1；(5).作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m组（m尽可能大），每组间隔为，记数各区间的随机数的数目jn，以为底，以nnj为高作第j(mj2,1)区间的矩形，最终构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。(6).以频率%951nnkjj为界划定区间，该区间半宽即为测量总误差的置信概率为95%的扩展不确定度。(7).合成的标准不确定度：112nvsunii4.实验流程图：生成[0,1]区间均匀分布的随机数利用Box-Muller方法生成正态分布的随机数根据各项误差的标准不确定度与正态分布的随机数生成各项误差分量的随机值获得误差和的随机值获得误差和的直方图计算95%的包罗面积对应的区间长度即为误差和的标准不确定度并于理论计算值比较，分析参数的变化对计算结果的影响一．实验1本实验中随机数种子为014。并使分别取N为100000点和10000点两种情况下，得到M值分别为5*N,2*N,N,N/2,N/5,N/10五种情况下的模拟图像。1．实验1程序tic;clear;clc;closeall;%%设定参数值%%%%随机信号点数N，均值为1，标准差u1,u2%%N=10^5;M=N/10;x=0:1:M;x_=[1:M];u1=0.005;u2=0.007;%%产生两个在(0,1)上服从均匀分布的,种子为0,每一次都相同的随机数X1和X2%%rand('state',014);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-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+1),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))/(N-1));%%%%%%%%%%%%%%%toc;2.实验1程序运行结果图（1）当M=N/10时Figure1Figure2Figure3Figure4Figure5Figure6（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如以下个图所示：图1.1N=10^5,M=N*5，s=0.0086，detla_n_u=0.0087图1.2N=10^5,M=N*2，s=0.0086，detla_n_u=0.0087图1.3N=10^5,M=N，s=0.0086，detla_n_u=0.0087图1.4N=10^5,M=N/2，s=0.0086，detla_n_u=0.0087图1.5N=10^5,M=N/5，s=0.0086，detla_n_u=0.0086图1.6N=10^5,M=N/10，s=0.0086，detla_n_u=0.0085图1.7N=10^4,M=N*5，s=0.0085，detla_n_u=0.0087图1.8N=10^4,M=N*2，s=0.0085，detla_n_u=0.0087图1.9N=10^4,M=N，s=0.0085，detla_n_u=0.0084图1.10N=10^4,M=N/2，s=0.0085，detla_n_u=0.0082图1.11N=10^4,M=N/5，s=0.0085，detla_n_u=0.0078图1.12N=10^4,M=N/10，s=0.0085，detla_n_u=0.0074表2N=10^5时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M1/51/212510s0.00860.00860.00860.00860.00860.0086delta_n_u0.00870.00870.00870.00870.00860.0085|delta_n_u﹣s|0.00010.00010.00010.000100.0001表2N=10^4时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M1/51/212510s0.00850.00850.00850.00850.00850.0085delta_n_u0.00870.00870.00840.00820.00780.0074|delta_n_u﹣s|0.00020.00020.00010.00030.00070.00113实验需要讨论的问题(1).N(采样点数)，M(划分的区间数)与直方图的关系(平滑，Y轴的高度)。有图1.1~1.12可知：当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。其中，MN时，Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间的误差随着M的增大而逐渐减小。当N改变时，即当N增大时可使直方图更为精细，且此时不改变直方图包络的基本形状。(2).Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间会存在误差，这个误差与哪些因素有关(N,M,N=M)此误差的大小和M、N的相对大小值有关。当N=M时，由于对离散的误差值统计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。当MN时，增大M值，误差值稳定，且不能改善误差值。二．实验2——自适应MCM法在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。(1).基于前一个实验，构建衡量MonteCarlo法和GUM法计算得到标准不确定度差值的函数。(2).将随机数个数N，分割区间数M分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。(3).分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N，M组合。1．实验2程序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));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;2.实验2程序运行结果图Figure1logspace(1,6)Figure2logspace(1,6)图2.1logspace(1,5)图2.1logspace(1,4)3实验需要讨论的问题如何根据三