哈工大误差分析课程设计MonteCarlo

1、Monte Carlo模拟误差分析课程设计1 实验目的1.1了解MATLAB软件的基本功能和使用1.2 学习不确定度的统计模拟分析方法1.3 研究误差概率密度函数和Bessel公式获得扩展不确定度的方法和影响因素2实验原理在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度(GUM)。在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。Monte Carlo(MCM)法就是较为常用的数学工具，具体原理相见相关资料。此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对

2、这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度(MCM)，并与理论上估计标准不确定度的Bessel公式、极差法作(GUM)比较，完成实验内容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。已知两项误差分量服从正态分布，标准不确定度分别为mV， mV，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率 P为99.73%的扩展不确定度)。3实验内容(1). 利用MATLAB软件生成0，1区间的均匀分布的随机数；(2). 给出误差分量的随机值：利用

3、MATLAB，由均匀分布随机数生成标准正态分布随机数，误差分量随机数可表示为mV；mV(3). 求和的随机数：误差和的随机数；(4). 重复以上步骤，得误差和的随机数系列： ；(5). 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为组（尽可能大），每组间隔为，记数各区间的随机数的数目，以为底，以为高作第()区间的矩形，最终构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。(6). 以频率为界划定区间，该区间半宽即为测量总误差的置信概率为95%的扩展不确定度。(7). 合成的标准不确定度：4.实验流程图：一实验1本实验中随机数种

4、子为014。并使分别取N为100000点和10000点两种情况下，得到M值分别为5*N, 2*N, N, N/2, N/5, N/10五种情况下的模拟图像。1实验1程序tic;clear;clc;close all;%设定参数值%随机信号点数N，均值为1，标准差u1,u2%N=105;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变换方法产生标准

5、正态分布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)p

6、lot(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;% 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,

7、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 on;%正态分布误差2幅度直方图%figure(4)axis(-1,1,0,N)hist(delta2,M);grid on;%合成误差幅度直方图%figure(5)axis(-1,1,0,N)H=hist(delta,M);hist(delta,M);grid on;%画包

8、络线%figure(6)HH=envelope(x_,H);plot(delta_n,HH,'b:');grid on;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)

9、+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; delta_n_u=(delta_n(floor(M/2+iii)-delta_n(floor(M/2-iii+1)/6;%理

10、论计算标准不确定度%delta_mean=mean(delta);delta_cancha=delta-delta_mean;s=sqrt(sum(delta_cancha.2)/(N-1);%toc;2. 实验1程序运行结果图（1）当M=N/10时Figure 1Figure 2Figure3Figure4Figure 5Figure 6（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如以下个图所示：图1.1 N=105, M=N*5，s=0.0086，detla_n_u=0.0087图1.2 N=105, M=N*2，s=0.0086，detla_n_u=0.0087图1.3 N

11、=105, M=N，s=0.0086，detla_n_u=0.0087图1.4 N=105, M=N/2，s=0.0086，detla_n_u=0.0087图1.5 N=105, M=N/5，s=0.0086，detla_n_u=0.0086图1.6 N=105, M=N/10，s=0.0086，detla_n_u=0.0085图1.7 N=104, M=N*5，s=0.0085，detla_n_u=0.0087图1.8 N=104, M=N*2，s=0.0085，detla_n_u=0.0087图1.9 N=104, M=N，s=0.0085，detla_n_u=0.0084图1.10 N=

12、104, M=N/2，s=0.0085，detla_n_u=0.0082图1.11 N=104, M=N/5，s=0.0085，detla_n_u=0.0078图1.12 N=104, M=N/10，s=0.0085，detla_n_u=0.0074表2 N=105时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M1/51/212510s0.00860.00860.00860.00860.00860.0086delta_n_u0.00870.00870.00870.00870.00860.0085|delta_n_us|0.00010.00010.00010.000100.

13、0001表2 N=104时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异k=N/M1/51/212510s0.00850.00850.00850.00850.00850.0085delta_n_u0.00870.00870.00840.00820.00780.0074|delta_n_us|0.00020.00020.00010.00030.00070.00113 实验需要讨论的问题(1). N(采样点数)，M(划分的区间数)与直方图的关系(平滑，Y轴的高度)。有图1.11.12可知：当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。其中，M<N时，Besse

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

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

16、); Result1(j,jj)=shiyan(a(j),b(jj); endendfigure(1),surfc(a,b,Result1);c=logspace(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;toc;2. 实验2程序运行结果图Figure 1 logspace(1,6)Figure 2 logspace(1,6)图2.1 logspace(1,5)图2.1 logspace(1,4)3 实

17、验需要讨论的问题如何根据三维网格曲线和三维曲线获得最佳的N，M组合。本实验中，N、M的取值，是由logspace(a,b,n)函数确定的。该函数的意义为生成从10的a次方到10的b次方之间按对数等分的n个元素的行向量。当logspace(a,b)时，则Matlab中默认n=50。 由以上实验可知，当M、N的取值越大，即logspace(a,b)中a与b的取值越大时，其误差越小，但是由于取值越大，要求的计算量也随之增加，所以N、M的取值使figure 2中蓝色图线随着横坐标的增加越快趋于平滑越好，且越接近原点处其斜率越大，波动越小越好。但是随着N、M值的增加，其计算量也会随着增加，这就要求计算机

18、达到一定的配置，所以要在计算机现有配置的基础上，选择最佳N、M的值。三实验3基于最短包含区间的MCM法如果PDF不对称，可采用最短包含区间。确定，使得，可得最短包含区间。(1). 先确定q(P=0.9973,M=1000000,q=PM=10000)(2). 重新计算包络线下的面积(不是对称的时候)(3). 根据算法：，计算r*(4). r*对应的区间长度极为最短包含区间1实验3程序%x为均匀分布，正态分布，反正弦分布%y=sin(x)为何种分布tic;clear;clc;close all; %设定参数值%随机信号点数N，均值1，标准差%N=106;M=N/10;x=0:1:M;x_=1:M

19、;u1=0.005;u2=0.007;%rand('state',014);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)h

20、ist(Y1,M);grid on;figure(2)axis(-1,1,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:');grid on;hold on;%计算直方图的面积%S=sum(d_delta*abs(H); % 计算直方图的面积%s_1表示正向直方图的每一个单元的面积%s_2表示反向直方图的每一个单元的面积%s_表示正反向两两对称每一对单

21、元的面积%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);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; delta_n_u1=(delta_n(floor(M/2+iii)