Monte-Carlo-模拟误差分析课程设计.doc

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. Monte Carlo模拟误差分析的实验原理
在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度(GUM)。在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。Monte Carlo(MCM)法就是较为常用的数学工具，具体原理相见相关资料。
此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度——(MCM)，并与理论上估计标准不确定度的Bessel公式、极差法作——(GUM)比较，完成实验内容。并以此作为基础，分析GUM法与MCM法的区别与联系，影响MCM法的参数，自适应MCM法和基于最短包含区间的MCM法。
已知两项误差分量服从正态分布，标准不确定度分别为mV， mV，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率 P为99.73%的扩展不确定度)。

3. Monte Carlo模拟误差分析的实验内容
3.1 MCM法与GUM法进行模拟误差分析和不确定度计算
(1). 利用MATLAB软件生成[0，1]区间的均匀分布的随机数；
(2). 给出误差分量的随机值：
利用MATLAB，由均匀分布随机数生成标准正态分布随机数，误差分量随机数可表示为
mV；
同理得
mV
(3). 求和的随机数：误差和的随机数；
(4). 重复以上步骤，得误差和的随机数系列： ；
(5). 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为组（尽可能大），每组间隔为，记数各区间的随机数的数目，以为底，以为高作第()区间的矩形，最终的组矩形构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。
(6). 以频率为界划定区间，该区间半宽即为测量总误差的置信概率为99.73%的扩展不确定度。
(7). 合成的标准不确定度：
A. 总误差随机数平均值

B. 各误差随机数的残差

C. 按照Bessel公式估计标准不确定度

实验流程图：

图6

实验说明：
本实验中随机数种子为31。以下为N分别为100000点和500000点两种情况下，M分别等于N/10、N/5、N/2、N、2N、5N六种情况下的模拟图像。

实验程序：
tic;
clear;clc;close all;
bdclose all;
%%设定参数值%%
%%随机信号点数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,每一次都相同的随机数X1和X2%%
rand('state',31);
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);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,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;
%%画包络线%
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)+s_2(i);

area(1)=s_(1);
for ii=1:M/2-1
area(ii+1)=area(ii)+s_(ii);
end

%% 计算99.73%的直方图面积
for iii=1:M/2;
area(iii);
if (area(iii)-0.9973*S)>=0;
break
end
end
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))/(N-1));
%%%%%%%%%%%%%%%
toc;

程序运行结果
（1）当N=100000，M=N/10时：

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

计算结果：s=0.0086，delta_n_u= 0.0089。

（2）当更改N与M不同的倍数关系时，可得到不同的计算结果，如下各图所示：

图3.1 N=100000,M=N/5; s=0.0086,delta_n_u= 0.0090

图3.2 N=100000,M=N/2; s=0.0086,delta_n_u= 0.0091

图3.3 N=100000, M=N; s=0.0086, delta_n_u= 0.0091

图3.4 N=100000,M=2N; s=0.0086,delta_n_u= 0.0091

图3.5 N=100000,M=5N; s=0.0086; delta_n_u= 0.0091

图3.6 N=500000,M=N/5; s=0.0086,delta_n_u= 0.0086

图3.7 N=500000,M=N/2; s=0.0086,delta_n_u= 0.0086

图3.8 N=500000,M=N; s=0.0086,delta_n_u= 0.0086

图3.9 N=500000,M=2N; s=0.0086,delta_n_u= 0.0086

图3.10 N=500000,M=5N; s=0.0086,delta_n_u= 0.0086

表1 N=100000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异
k=N/M
10
5
2
1
1/2
1/5
s
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
delta_n_u
0.0089
0.0090
0.0091
0.0091
0.0091
0.0091
|delta_n_u﹣s|
0.0003
0.0004
0.0005
0.0005
0.0005
0.0005

表2 N=500000时，N与M成不同倍数k时，直方图计算结果与理论计算结果的差异
k=N/M

10
5
2
1
1/2
1/5
s
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
delta_n_u
0.0086
0.0086
0.0086
0.0086
0.0086
0.0086
|delta_n_u﹣s|
0
0
0
0
0
0

实验需要讨论的问题：
（1）N(采样点数)，M(划分的区间数)与直方图的关系(平滑，Y轴的高度)。
根据以上各图分析知：当N固定的情况下，随着M值得增大直方图的平滑性变差，Y轴高度下降。其中，M图面积的扩展不确定度两者之间的误差随着M的增大而逐渐减小。当N改变时，即当N增大时可使直方图更为精细，且此时不改变直方图包络的基本形状。
（2） Bessel公式计算的标准不确定度与99.73%直方图面积的扩展不确定度两者之间会存在误差，这个误差与哪些因素有关(N,M,N>=M)
此误差的大小和M、N的相对大小值有关。当N>=M时，由于对离散的误差值统计运算存在舍入误差导致误差，此误差随着M的增大可消除此项舍入误差。当M>N时，增大M值，误差值稳定，但不能改善误差值。
3.2 自适应MCM法
在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。
(1). 基于前一个实验，构建衡量Monte Carlo法和GUM法计算得到标准不确定度差值的函数。
(2). 将随机数个数N，分割区间数M分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。
(3). 分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N，M组合。
实验示例程序：
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));
end
end
figure(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));
end
figure(2),plot3(c,d,Result2);grid on;
toc;
其中shiyan（）子程序为：
function y=shiyan(N,M)
%%clear;clc;close all;
%%bdclose all;
%%设定参数值%%
%%随机信号点数N，均值为1，标准差u1,u2%%
%N=10^5;
Mu=1;
%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*Y2;
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);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,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,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_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)-1
area(ii+1)=area(ii)+s_(ii);
end

%% 计算99.73%的直方图面积
for iii=1:floor(M./2);
area(iii);
if (area(iii)-0.9973*S)>=0;
break
end
end
%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));
y=abs(delta_n_u-s);
程序运行结果：

Figure 1

Figure 2
实验需要讨论的问题：
如何根据三维网格曲线和三维曲线获得最佳的N，M组合。
根据shiyan（）子函数知：程序返回值为y=abs(delta_n_u-s)；显然，当y=0时即可获得N，M的最佳组合，即三维网格曲线和三维曲线的Z坐标为0时的N，M。
3.3 基于最短包含区间的MCM法
如果PDF不对称，可采用最短包含区间。确定，使得，，可得最短包含区间。
(1). 先确定q(P=0.9973,M=1000000,q=PM=10000)
(2). 重新计算包络线下的面积(不是对称的时候)
(3). 根据算法：，，计算r*
(4). r*对应的区间长度极为最短包含区间
实验示例程序：
%%x为均匀分布，正态分布，反正弦分布
%%y=sin(x)为何种分布
tic;
clear;clc;close all;

%%设定参数值%%
%%随机信号点数N，均值1标准差%%
N=10^6;
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*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;
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);
plotdelta_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)+s_2(i);

area(1)=s_(1);
for ii=1:M/2-1
area(ii+1)=area(ii)+s_(ii);
end

%% 计算99.73%的直方图面积
for iii=1:M/2;
area(iii);
if (area(iii)-0.9973*S)>=0;
break
end
end
plot([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/2iii))-delta_n(floor(M/2-iii+1)))/6;

%% 计算最短包含区间面积
P=0.9973; %置信概率
q=P*M; % 参数q
iiii=1:M;
s2(iiii)=d_delta*abs((H(iiii)));
area2(1)=s2(1);
for j=1:M-1
area2(j+1)=area2(j)+s2(j);
end
%%%%%%
for jj=1:M-q
area2(jj);
for r=1:M/2
if area2(r+q)-area2(r)<=area2(jj+q)area2(jj);
break
end
end
end
plot([delta_n(r),delta_n(r+q)],[H(r),H(r+q)],'r*');grid on;
delta_n_u2=(delta_n(floor(r+q))-delta_n(floor(r)))/6
%%理论计算标准不确定度%%
delta_mean=mean(delta);
delta_cancha=delta-delta_mean;
s=sqrt((sum(delta_cancha.^2))/(N-1))
%%%%%%%%%%%%%%%
toc;
程序运行结果：

Figure 1

Figure 2

Figure 3

Figure 4
计算结果：delta_n_u2 =0.0122，s = 0.0071
实验需要讨论的问题：
(1) 最短包含区间MCM法与普通MCM法有何区别。最短包含区间MCM适合解决怎样的问题。
适合解决概率密度函数非对称的问题
(2) 最短包含区间与GUM法同时计算不确定度的时候，哪一个略显保守。
由计算得delta_n_u2 =0.0113，s = 0.0071，显然，在评定不确定度的时候，最短包含区间略显保守。

22