蒙特卡罗Monte Carlo模拟误差分析课程设计哈工大误差

1、Monte Carlo模拟误差分析课程设计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. MATLAB软件介绍实验内容2.1 介绍MATLAB软件的基本知识MATLAB名字由MATrix和LABoratory 两词的前三个字母组合而

2、成。20世纪七十年代，时任美国新墨西哥大学计算机科学系主任的Cleve Moller出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK矩阵软件工具包库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLABMATLAB语言的主要特点(1). 具有丰富的数学功能(2). 具有很好的图视系统(3). 可以直接处理声言和图形文件。(4). 具有若干功能强大的应用工具箱。(5). 使用方便，具有很好的扩张功能。(6). 具有很好的帮助功能演示内容：(1). MATLAB的数值计算功能在“命令行”Command提示窗口中键入：“A=eye(5,5);A=ze

3、ros(5,5);A=ones(5,5)”等命令生成各类矩阵；在“命令行”Command提示窗口中键入：“v,d=eig (A)”生成特征矩阵和特征向量；在“命令行”Command提示窗口中键入：“expm(A)”对矩阵A求幂；在“命令行”Command提示窗口中键入：x=1 3 5;y=2 4 6;z=conv(x,y)；显示结果：z = 2 10 28 38 30(2). MATLAB的符号计算功能在“命令行”Command提示窗口中键入：syms a x；f=sin(a*x); df=diff(f,x); dfa=diff(f,a)；Command提示窗口显示结果： df =cos(a*

4、x)*a； dfa =cos(a*x)*x；2.2 MATLAB软件画图特性(1). MATLAB二维绘图命令函数：plot参数：线型、颜色、多重线、网格和标记、画面窗口分割、其他方式、隐函数的描绘)(2). MATLAB三维画图曲面与网格图命令函数：mesh三维带阴影曲面图：surf三维曲线命令：plot3演示内容：(1). MATLAB的二维绘图功能在命令行Command提示窗口中键入：close all; x=linspace(0, 2*pi, 100); % 100个点的x座标y=sin(x); % 对应的y座标plot(x,y);得到如下的结果：10.80.60.40.2-0.2-0

5、.4-0.6-0.8-1图1在命令行Command提示窗口中键入：“plot(x, sin(x), x, cos(x);”得到如下的结果：10.80.60.40.2-0.2-0.4-0.6-0.8-101234567图2在命令行Command提示窗口中键入：plot(x, sin(x), 'co', x, cos(x), 'g*');得到如下的结果：图3在命令行Command提示窗口中键入：xlabel('Input Value'); % x轴注解ylabel('Function Value'); % y轴注解title('

6、;Two Trigonometric Functions'); % 图形标题 legend('y = sin(x)','y = cos(x)'); % 图形注解grid on; % 显示格线得到如下的结果：Two Trigonometric FunctionsFunction ValueInput Value 图4(2). MATLAB的多维绘图功能在命令行Command提示窗口中键入：X,Y = meshgrid(-3:0.125:3); % 生成二维网格点 Z = peaks(X,Y); % 生成某种内置函数 mesh(X,Y,Z);得到如下的结果-

7、4图5其他的演示功能详见“MATLAB画图文档”3. Monte Carlo模拟误差分析的实验原理在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度(GUM)。在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。Monte Carlo(MCM)法就是较为常用的数学工具，具体原理相见相关资料。此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度(MCM)，并与理论上估计标准不

8、确定度的Bessel公式、极差法作(GUM)比较，完成实验内容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。已知两项误差分量服从正态分布，标准不确定度分别为u1=5mV， u2=7mV，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率 P为99.73%的扩展不确定度)。4. Monte Carlo模拟误差分析的实验内容4.1 MCM法与GUM法进行模拟误差分析和不确定度计算(1). 利用MATLAB软件生成0，1区间的均匀分布的随机数；(2). 给出误差分量的随机值：利用MATLAB

9、，由均匀分布随机数1生成标准正态分布随机数1，误差分量随机数可表示为1=1u1=51mV；同理得2=2u2=72mV(3). 求和的随机数：误差和的随机数=1+2；(4). 重复以上步骤，得误差和的随机数系列：i=1i+2i i=1,2, n；(5). 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m组（m尽可能大），每组间隔为，记数各区间的随机数的数目nj，以为底，以njn为高作第j(j=1,2 m)区间的矩形，最终的m组矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。n(6). 以频率j=1kjn=99.73%为界划

10、定区间，该区间半宽即为测量总误差的置信概率为99.73%的扩展不确定度。(7). 合成的标准不确定度：A. 总误差随机数平均值1n=i ni=1B. 各误差随机数的残差vi=i-C. 按照Bessel公式估计标准不确定度u=s=实验流程图：vi=1n2in-1图6实验1本实验中随机数种子为28。以下为N分别为100000点和500000点两种情况下，M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。1）程序tic;clear;clc;close all;%设定参数值%随机信号点数N，均值为1，标准差u1,u2%N=100000;M=N/10;x=0:1:M;x_=1:M;

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

12、elta)/(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);grid on;title('学号：13S101028 姓名：李明');figure(2),axis(0,N,-max(5*Y2),max(5*Y2)plot(Y2);grid on;title('学号：13S101028 姓名：李明');% hold on% plot(x,0,

13、'k');grid on;% plot(x,1,'r-');grid on;% plot(x,-1,'r-');grid on;% hold on%变换为任意均值和方差的正态分布%Z1=Sigma*Y1+Mu;%作图%高斯随机信号% subplot(2,2,2)% axis(0,N,-6,6)% plot(Z1);grid on;% hold on% plot(x,Mu,'k');% plot(x,Mu+Sigma,'r-');grid on;% plot(x,Mu-Sigma,'r-');grid

14、 on;% hold on%正态分布误差1幅度直方图%figure(3)axis(-1,1,0,N)hist(delta1,M);grid on;title('学号：13S101028 姓名：李明'); %正态分布误差2幅度直方图%figure(4)axis(-1,1,0,N)hist(delta2,M);grid on;title('学号：13S101028 姓名：李明'); %合成误差幅度直方图%figure(5)axis(-1,1,0,N)H=hist(delta,M);hist(delta,M);grid on;title('学号：13S1010

15、28 姓名：李明'); %画包络线%figure(6)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;title('学号：13S101028 姓名：李明'); hold on;%计算直方图的面积%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

16、);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);for ii=1:M/2-1area(ii+1)=area(ii)+s_(ii);end% 计算99.73%的直方图面积for iii=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');grid on; title(

17、9;学号：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不同的倍数关系时，可得到不同的计算结果，如下各图所

18、示：图7 N=100000, M=N/5; s=0.0086, delta_n_u= 0.0086图8 N=100000, M=N/2; s=0.0086,delta_n_u= 0.0086图9 N=100000, M=N; s=0.0086, delta_n_u= 0.0086图10 N=100000, M=2N; s=0.0086,delta_n_u= 0.0086图11 N=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

19、.0088。（4）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示：图18 N=500000, M=N/5; s=0.0086,delta_n_u= 0.0088图19 N=500000, M=N/2; s=0.0086, delta_n_u= 0.0088图20 N=500000, M=N; s=0.0086,delta_n_u= 0.0088图21 N=500000,M=2N; s=0.0086, delta_n_u= 0.0088图22 N=500000, M=5N; s=0.0086, delta_n_u=0.0088表1 N=100000时，N与M成不同倍数k时，直方

20、图计算结果与理论计算结果的差异表2 N=500000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异3）实验需要讨论的问题（1）N(采样点数)，M(划分的区间数)与直方图的关系(平滑，Y轴的高度)。根据以上各图分析知：当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。其中，M<N时，Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间的误差随着M的增大而逐渐减小。当N改变时，即当N增大时可使直方图更为精细，且此时不改变直方图包络的基本形状。（2） Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间会存在

21、误差，这个误差与哪些因素有关(N,M,N>=M)此误差的大小和M、N的相对大小值有关。当N>=M时，由于对离散的误差值统计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。当M>N时，增大M值，误差值稳定，但不能改善误差值。 实验2自适应MCM法在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。(1). 基于前一个实验，构建衡量Monte Carlo法和GUM法计算得到标准不确定度差值的函数。(2). 将随机数个数N，分割区间数M分

22、别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。(3). 分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N，M组合。1）程序tic;warning off;a,b=meshgrid(logspace(1,6);for j=1:max(size(a);for jj=1:max(size(b);Result1(j,jj)=shiyan(a(j),b(jj);endendfigure(1),surfc(a,b,Result1);title('学号：13S101028 姓名：李明'); c=logspace(1,6);d=logspace(1,6);for

23、 jjj=1:max(size(c);Result2(jjj)=shiyan(c(jjj),d(jjj);endfigure(2),plot3(c,d,Result2);grid on;title('学号：13S101028 姓名：李明'); toc;其中shiyan（）子程序为：function y=shiyan(N,M)%clear;clc;close all;%bdclose all;%设定参数值%随机信号点数N，均值为1，标准差u1,u2%N=105;Mu=1;%M=N/10;x=0:1:floor(M);x_=1:floor(M);u1=0.005;u2=0.007;

24、%产生两个在(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).

25、/(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);grid on;% figure(2),% axis(0,N,-max(5*Y2),max(5*Y2)% plot(Y2);grid on;% % hold on% % plot(x,0,'k');grid on;% % plot(x,1,'r-');grid on

26、;% % plot(x,-1,'r-');grid on;% % hold on% %变换为任意均值和方差的正态分布% % %Z1=Sigma*Y1+Mu;% %作图% %高斯随机信号% % subplot(2,2,2)% % axis(0,N,-6,6)% % plot(Z1);grid on;% % hold on% % plot(x,Mu,'k');% % plot(x,Mu+Sigma,'r-');grid on; % % plot(x,Mu-Sigma,'r-');grid on;% % hold on% %正态分布误差

27、1幅度直方图% % figure(3)% axis(-1,1,0,N)% hist(delta1,M);grid on;% %正态分布误差2幅度直方图% % figure(4)% axis(-1,1,0,N)% hist(delta2,M);grid on;% %合成误差幅度直方图% figure(5)% axis(-1,1,0,N)H=hist(delta,floor(M);% hist(delta,M);grid on;% %画包络线% figure(6)% HH=envelope(x_,H);% plot(delta_n,HH,'b:');grid on;% hold o

28、n;%计算直方图的面积%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);for ii=1:floor(M./2)-1area(ii+1)=area(i

29、i)+s_(ii);end% 计算99.73%的直方图面积for iii=1:floor(M./2);area(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');grid on; delta_n_u=(delta_n(floor(M./2+iii)-delta_n(floor(M./2-iii+1)./6;%理论计算标准不确定度% delta_mean=mean(delta);delta_canch

30、a=delta-delta_mean;s=sqrt(sum(delta_cancha.2)./(floor(N)-1); y=abs(delta_n_u-s);2）程序运行结果图23图243）实验需要讨论的问题：如何根据三维网格曲线和三维曲线获得最佳的N，M组合。根据shiyan（）子函数知：程序返回值为y=abs(delta_n_u-s)；显然，当y=0时即可获得N，M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N，M。实验3基于最短包含区间的MCM法如果PDF不对称，可采用最短100p%包含区间。确定r*，使得y(r*+q)-y(r*)y(r+q)-y(r)，r=1,M-q，可得最

31、短100p%包含区间y(r),y(r+q)。 *(1). 先确定q(P=0.9973,M=1000000,q=PM=10000)(2). 重新计算包络线下的面积(不是对称的时候)(3). 根据算法：y(r*+q)-y(r*)y(r+q)-y(r)，r=1,M-q，计算r*(4). r*对应的区间长度极为最短包含区间1）程序%x为均匀分布，正态分布，反正弦分布%y=sin(x)为何种分布tic;clear;clc;close all;%设定参数值%随机信号点数N，均值1，标准差%N=106;M=N/10;x=0:1:M;x_=1:M;u1=0.005;u2=0.007;%rand('st

32、ate',28);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);grid on;title('学号

33、：13S101028 姓名：李明');figure(2)axis(-1,1,0,N)hist(Y2,M);grid on;title('学号：13S101028 姓名：李明');figure(3);axis(-1,1,0,N)hist(delta,M);grid on;title('学号：13S101028 姓名：李明');H=hist(delta,M);%画包络线%figure(4)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;title('学号：13S101028 姓名：李明');hold on;%计算直方图的面积%S=sum(d_delta*abs(H);% 计算直方图的面积%s_1表示正向直方图的每一个单元的面积%s_2表示反向直方图的每一个单元的面积%s_表示正反向两两对称每一对单元的面积%area表示以中心为对称轴的累加面积i=1:1:M/2;s_1(i)=d_del