蒙特卡罗MonteCarlo模拟误差分析课程设计哈工大误差原理课程设计报告书

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提示窗口

3、中键入： “A=eye(5,5);A=zeros(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 = 210283830(2). MATLAB 的符号计算功能在“命令行 ”Command提示窗口中键入：symsax；f=sin(a*x); df=diff(f,x); dfa=diff(

4、f,a)；Command 提示窗口显示结果：df =cos(a*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); %

5、 对应的 y 座标plot(x,y);得到如下的结果：10.80.60.40.20-0.2-0.4-0.6-0.8-101234567图 1在命令行 Command 提示窗口中键入：“plot(x, sin(x), x, cos(x); ”得到如下的结果：10.80.60.40.20-0.2-0.4-0.6-0.8-101234567图 2在命令行 Command 提示窗口中键入：plot(x, sin(x), co, x, cos(x), g*);得到如下的结果：10.80.60.40.20-0.2-0.4-0.6-0.8-101234567图 3在命令行 Command 提示窗口中键入：x

6、label(Input Value); % x 轴注解ylabel(Function Value); % y 轴注解title(Two Trigonometric Functions); %图形标题legend(y = sin(x),y = cos(x); % 图形注解grid on; % 显示格线得到如下的结果：eulaVnoticnuFTwo Trigonometric Functions1y = sin(x)0.8y = cos(x)0.60.40.20-0.2-0.4-0.6-0.8-101234567Input Value图 4(2). MATLAB 的多维绘图功能在命令行 Comm

7、and 提示窗口中键入：X,Y = meshgrid(-3:0.125:3); %生成二维网格点Z = peaks(X,Y);% 生成某种置函数mesh(X,Y,Z);得到如下的结果1050-5-10424020-2-2-4-4图 5其他的演示功能详见 “MATLAB画图文档 ”3. Monte Carlo 模拟误差分析的实验原理在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度(GUM) 。在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。 Monte Carlo(MC

8、M) 法就是较为常用的数学工具，具体原理相见相关资料。此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度 (MCM) ，并与理论上估计标准不确定度的Bessel公式、极差法作 (GUM)比较，完成实验容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。已知两项误差分量服从正态分布，标准不确定度分别为u15 mV ，u27 mV ，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率P 为

9、99.73%的扩展不确定度 )。4. Monte Carlo 模拟误差分析的实验容4.1 MCM 法与 GUM 法进行模拟误差分析和不确定度计算(1). 利用 MATLAB 软件生成 0 ，1区间的均匀分布的随机数；(2). 给出误差分量的随机值：利用 MATLAB ，由均匀分布随机数1 生成标准正态分布随机数1 ，误差分量随机数可表示为11u15 1 mV ；同理得22 u27 2 mV(3). 求和的随机数：误差和的随机数12 ；(4). 重复以上步骤，得误差和的随机数系列：i1i2ii 1,2,n ；(5). 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所

10、有数据，按数值将区间等分为m 组（ m 尽可能大），每组间隔为，记数各区间的随机数的数目n j ，以为底，以 n j 为高作第 j ( j1,2m )区间的矩形，最终的 m 组n矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。knj(6). 以频率 j 199.73% 为界划定区间，该区间半宽即为测量总误差的置信概率为n99.73%的扩展不确定度。(7). 合成的标准不确定度：A. 总误差随机数平均值1nnii1B. 各误差随机数的残差viiC. 按照 Bessel公式估计标准不确定度nvi2i1usn1实验流程图：生成 0,1 区间均匀分布的随机数利用 Box-Muller

11、方法生成正态分布的随机数根据各项误差的标准不确定度与正态分布的随机数生成各项误差分量的随机值图 6计算 99.73%的包罗面积对应的区间长度即为误差和的标准不确定度并于理论计算值比较，分析参数的变化对计算结果的影响获得误差和的直方图获得误差和的随机值实验 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

12、;x_=1:M;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

13、)-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);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,k);grid on;% pl

14、ot(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 on;% hold on%正态分布误差 1幅度直方图 %figure(3)axis(-1,1,0,N)hist(delta1,M);grid o

15、n;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( 学号： 13S101028 ：明 );%画包络线 %figure(6)HH=envelope(x_,H);plot(delta_n,HH,b:);grid on;title( 学号： 13S101028 ：明

16、 );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);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);

17、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(号： 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-d

18、elta_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 不同的倍数关系时，可得到不同的计算结果，如下各图所示：图 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

19、_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.0088。（ 4）当更改 N 与 M 不同的倍数关系时，可得到不同的计算结果，如下各图所示：图 18 N=500000, M=N/5; s=0.0086, delta_n_u= 0.0088图 19 N=500000, M=N

20、/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 时，直方图计算结果与理论计算结果的差异k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086delta_n_u0.00850.00860.00860

21、.00860.00860.0086|delta_n_u s|0.000100000表 2 N=500000 时， N 与 M 成不同倍数k 时，直方图计算结果与理论计算结果的差异k=N/M105211/21/5s0.00860.00860.00860.00860.00860.0086delta_n_u0.00880.00880.00880.00880.00880.0088|delta_n_u s|0.00020.00020.00020.00020.00020.00023）实验需要讨论的问题（ 1）N( 采样点数 )，M( 划分的区间数 )与直方图的关系 (平滑， Y 轴的高度 )。根据以上各图

22、分析知： 当 N 固定的情况下， 随着 M 值得增大直方图的平滑性变差，Y 轴高度下降。其中， M=M)此误差的大小和 M 、 N 的相对大小值有关。当N=M 时，由于对离散的误差值统计运算存在舍入误差导致误差， 此误差随着 M 的增大可消除此项舍入误差。 当 MN 时，增大 M 值，误差值稳定，但不能改善误差值。实验 2 自适应 MCM 法在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。(1). 基于前一个实验，构建衡量Monte Carlo 法和 GUM 法

23、计算得到标准不确定度差值的函数。(2). 将随机数个数 N，分割区间数 M 分别作为该函数的自变量，定义自变量的取值围，从而获得相应的函数值。(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=logspa

24、ce(1,6);d=logspace(1,6);for 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

25、(M);u1=0.005;u2=0.007;%产生两个在 (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(d

26、elta)-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);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;%

27、% 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% %正态分布误差 1 幅度直方图 % figure(3)% axis(-1,1,0,N)% hist(

28、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 on;%计算直方图的面积 %S=sum(d_delta.*abs(H);% 计算直方图的面积 %s

29、_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(ii)+s_(ii);end% 计算 99.73%的直方图面积for iii=1:floo

30、r(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_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.2)./(floor(N)-1

31、);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 不对称，可采用最短100 p%包含区间。确定r * ，使得y(r*q )y( r* )y( r q)y(r )， r1, Mq ，可得最短100p% 包含区间y( r * ) , y( r*q )。(1). 先确定

32、 q(P=0.9973,M=1000000,q=PM=10000)(2). 重新计算包络线下的面积(不是对称的时候 )(3). 根据算法： y( r* q) y( r* )y(r q ) y( r ) ， r 1, , Mq ，计算 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(state,28);X1=rand(

33、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( 学号： 13S101028 ：明 );figure(2)axi

34、s(-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_表示正反