仪器精度实验

仪器精度理论实验报告

仪器精度理论实验报告

学 院： 专 业： 姓 名： 学 号：

仪器精度理论实验报告

1 Monte Carlo模拟误差分析原理

在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度，这种方法称为GUM 法。GUM 法适用于误差是线性的，然而在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。Monte Carlo(MCM)法就是较为常用的数学工具。

两种方法的不确定度统计模拟分析方法如下： 设有测量方程

y =f (x 1, x 2,..., x n ) (1-1)

以计算机模拟m 次测量的各误差分量δx 1j , δx 2j , …, δx nj , j =1,2,…, m 。则测量数据可写成

x 1j =X 1+δx

j 1j 2

x 2j =X 2+δx x nj =X n +δx nj

j =1,2,…, m (1-2)

式中，X 1, X 2, , X n 为各项分量的真值。将各项数据带入测量方程，得

y j =f (x 1j , x 2j , x nj )

=f (X 1+δx 1j , X 2+δx 2j , , X n +δx nj )

(1-3)

该结果也可写为

y j =Y +δy j

=f (X 1+δx 1j , X 2+δx 2j , , X n +δx nj )

(1-4)

式中，Y 为被测量的真值，y j 为m 个测量数据，j =1,2,…, m ; δy j 为由各项误差分量

δx 1j , δx 2j , …, δx nj , j =1,2,…, m 造成的测量结果y j 的误差。

按GUM 法，由m 个测量数据y j ，取其算数平均值得

1m

y =∑y j (1-5)

m j =1

计算其残差

1m

v j =y j -y =y j -∑y j (1-6)

m j =1

由Bessel 公式得测量数据的标准差为

仪器精度理论实验报告

按MCM 法，将y j 以适宜的子区间间隔绘制成直方图，得到频率分布。对于置信概率p ，若pm 为整数，设q =pm ，否则取q 的值为pm +2的整数部分。则⎡⎣y low , y high ⎤⎦为Y 的其中，对任意的r =1, , m -q ，y low =y (r ) ，y high =y (r +q ) 。对于对称的PDF ，100p %包含区间，

如GUM 法中使用的正态分布和缩放位移t 分布，概率对称100p%包含区间和最短100p%包含区间是相同的。如果PDF 不对称，可采用最短100p %包含区间。确定r *，使得

y (r *+q ) -y (r *) ≤y (r +q ) -y (r ) ，r =1, , M -q ，可得最短100p %包含区间⎡y (r ) , y (r +q ) ⎤。

⎣

*

*

⎦

在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。

用MATLAB 可以在短时间内获得大量的误差分量样本，可以很好的模拟统计实验。本次实验包含三个小实验：

（1）MCM 法与GUM 法的实现及比较； （2）自适应MCM 法的实现；

（3）基于最短包含区间的MCM 法的实现。

仪器精度理论实验报告

2 MCM法与GUM 法的比较实验

2.1 实验步骤

GUM 法与MCM 法的实现步骤如下：

(1) 利用MATLAB 的rand 函数生成[0，1]区间内的均匀分布的随机数ξ1和ξ2； (2) 给出误差分量的随机值：利用MATLAB ，由均匀分布随机数ξ1生成标准正态分布随机数η1，误差分量随机数可表示为

δ1=η1u 1=5η1mV (2-1)

同理得

δ2=η2u 2=7η2mV (2-2)

(3) 求和的随机数：误差和的随机数δ=δ1+δ2；

(4) 重复以上步骤，得误差和的随机数系列：δi =δ1i +δ2i i =1, 2, n ；

(5) 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m 组（m 尽可能大），每组间隔为∆，记数各区间的随机数的数目n j ，以∆为底，以

n j n ∆

为高作第j (j =1, 2 m ) 区间的矩形，最终的m 组矩形构

∑n

成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。以频率

j =1

k

j

n

=99.73%为

界划定区间，该区间半宽即为测量总误差的置信概率为99.73%的扩展不确定度。

(6) 合成的标准不确定度：总误差随机数平均值为

1n

δ=∑δi (2-3)

n i =1

各误差随机数的残差

v i =δi - (2-4)

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

u =s =

∑v

i =1

n

2i

n -1

(2-5)

实验流程图如图2-1所示：

仪器精度理论实验报告

图2-1 实验流程图

2.2 程序实现

用MATLAB 实现如下： tic;

clear;clc;close all; bdclose all; %%设定参数值%%

%%随机信号点数N ，均值为1，标准差u1,u2%% N=10^6; M=N/100; x=0:1:M; x_=[1:M]; u1=0.005; u2=0.007;

%%产生两个在(0,1)上服从均匀分布的, 种子为53, 每一次都相同的随机数X1和X2 rand('state',53); 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;

仪器精度理论实验报告

lta 为误差分布的间距

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;

%%正态分布误差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(abs(H))/N;

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积

仪器精度理论实验报告

s_1(i)=1*abs(floor(H(floor(M/2+i))))/N; s_2(i)=1*abs(floor(H(floor(M/2-i+1))))/N; 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_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;

种子设定为学号后两位53，程序运行后，得到合成误差delta 概率分布的三组曲线： (1) N=106，M=104； (2) N=106，M=103； (3) N=105，M=104。

分别如图2-2(a)，2-2(b)，2-2(c)所示。

仪器精度理论实验报告

-0.05

-0.04

-0.03

-0.02

-0.01

0.01

0.02

0.03

0.04

0.05

图2-2(a) N=10, M=10delta 概率分布

6

4

-0.05

-0.04-0.03-0.02-0.0100.010.020.030.040.05

图-2(b) N=10, M=10delta 概率分布

6

3

仪器精度理论实验报告

-0.04

-0.03-0.02-0.0100.010.020.030.04

图2-2(c) N=10, M=10delta 概率分布

5

4

三种情况下，算出的GUM 法和MCM 法得到的不确定度分别为： (1) 图(a)，delta_n_u =0.0084，s =0.0086； (2) 图(b)，delta_n_u =0.0073，s =0.0086； (3) 图(c)，delta_n_u =0.0084，s =0.0086。

2.3 问题讨论

(1) N(采样点数) ，M(划分的区间数) 与直方图的关系(平滑，Y 轴的高度) 。

答：通过对比图(a)与图(b)可以看出：划分的区间M 越小，Y 轴高度越高，曲线越平滑。因为当采样点数N 不变时，区间数M 减少导致单个区间的频数变大，所以区间的过渡更接近理想正态分布曲线。同时，通过对比图(a)与图(c)可以看出：采样点数N 越大，Y 轴高度越高，曲线越平滑。因为当划分的区间数M 不变时，采样点数N 增加导致单个区间的频数变大，所以区间的过渡更接近理想正态分布曲线。

(2)Bessel公式计算的标准不确定度(GUM)与99.73%直方图面积的扩展不确定度(MCM)两者之间会存在误差，这个误差与哪些因素有关(N,M,N与M 的关系) 。

答：通过对比由(1)和(2)得到的参数可知：取N=106，M=104时，u 和s 的误差不大；取N=106，M=103时，u 和s 的误差很大。再通过对比由(1)和(3)得到的参数可知：取N=106，M=104时，u 和s 的误差不大；取N=105，M=104时，u 和s 的误差更小。综上说明：采样点数N 对误差的影响不大，划分区间数M 对误差的影响很大，且N 与M 的比值越大，误差越大。

仪器精度理论实验报告

3 自适应MCM 法实验

3.1 实验步骤

自适应MCM 法的步骤如下：

(1) 基于前一个实验，构建衡量Monte Carlo 法和GUM 法计算得到标准不确定度差值的函数。

(2) 将随机数个数N ，分割区间数M 分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。

(3) 分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N ，M 组合。

3.2 程序实现

用MATLAB 实现如下： tic;

close all;clear;clc;

%根据三维网格曲线求最佳的N ，M 组合 %a,b模拟N 和M a=logspace(1,6); b=logspace(1,6);

%j代表N ，jj 代表M ，且N>M for j=1:max(size(a)); for jj=1:j;

Result1(j,jj)=shiyan(a(j),b(jj)); end end

%求Result1的非零最小值，并返回行号和列号 Result1(find(Result1==0))=NaN; f=min(min(Result1)); [n,m]=find(Result1==f);

[a,b]=meshgrid(logspace(1,6)); %求对应行列的N 、M M=a(n,m); N=b(n,m); %作三维网格曲线 figure(1),surfc(a,b,Result1); xlabel('M');

仪器精度理论实验报告

zlabel('Δ');

%根据三维曲线求最佳的N ，M 组合 %c、d 模拟N 、M

c=logspace(1,6); d=logspace(1,6);

for jjj=1:max(size(c));

Result2(jjj)=shiyan(c(jjj),d(jjj)); end

%求Result2的最小值，返回最小值和列号 [Y,nm]=min(Result2); NM=c(nm); %作三维曲线

figure(2),plot3(c,d,Result2);grid on; xlabel('M'); ylabel('N'); zlabel('Δ'); toc;

其中函数shiyan 的代码如下，以N 、M 为参量，以GUM 法与MCM 法所得不确定度之差作为返回值，种子数为53。 function y=shiyan(N,M) %%设定参数值%%

%%随机信号点数N ，均值为1，标准差u1,u2 Mu=1;

x=0:1:floor(M); x_=[1:floor(M)]; u1=0.005; u2=0.007;

%%产生两个在(0,1)上服从均匀分布的, 种子为53, 每一次都相同的随机数X1和X2%% rand('state',53); 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;

仪器精度理论实验报告

delta=delta1+delta2;

d_delta=(max(delta)-min(delta))./(floor(M)-1); %%d_delta为误差分布的间距 delta_n=[min(delta):d_delta:max(delta)]; %%delta_n为误差分布序列 H=hist(delta,floor(M)); %%计算直方图的面积%% S=sum(1.*abs(H))/N; %% 计算直方图的面积%%

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积 i=1:1:floor(M./2);

s_1(i)=1*abs(floor(H(floor(M./2+i))))/N; s_2(i)=1.*abs(floor(H(floor(M./2-i+1))))/N; 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

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)); %%计算GUM 法与MCM 法所得不确定度的差值 y=abs(delta_n_u-s);

仪器精度理论实验报告

出Δ的最小值及对应的M 、N 值，如图3-1(a)所示。

Δ

M

5

N

图3-1(a) 三维网格曲线图

三维曲线如图3-1(b)所示，X 轴和Y 轴分别为M 、N ，Z 轴为Δ=abs(delta_n_u-s)，并得到Δ的最小值及M 、N 。

Δ

M

5

N

图3-1(b) 三维曲线图

3.3 问题讨论

如何根据三维网格曲线和三维曲线获得最佳的N ，M 组合。

仪器精度理论实验报告

线，其误差Δ的最小值为3.6949×10-8，在矩阵中处于43行49列，对应的M 和N 分别为M=1.9307×105，N=7.9060×105。对于三维曲线，其误差Δ的最小值为1.6837×10-6，处于41行41列，对应的M 和N 为1.2068×105。

仪器精度理论实验报告

4 基于最短包含区间的MCM 法实验

4.1 实验步骤

基于最短包含区间的MCM 法的实现步骤如下： (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*对应的区间长度极为最短包含区间。

4.2 程序实现

用MATLAB 实现如下： tic;

clear;clc;close all; bdclose all; %%设定参数值%%

%%随机信号点数N ，均值1，标准差%% N=10^6; M=N/100; x=0:1:M; x_=[1:M]; u1=1; u2=1;

%%产生两个在(-1,1)内的均匀分布, 种子为53, 每一次都相同的随机数X1和X2 rand('state',53);

X1=unifrnd(-1,1,1,N); X2=unifrnd(-1,1,1,N);

%%产生正弦分布Y1和正态分布Y2 Y1=sin(X1);

Y2=sqrt(-2*log(abs(X1))).*cos(2*pi*X2); %% 为做直方图先定义好X 轴的坐标数据%% delta1=u1*Y1; delta2=u2*Y2;

delta=delta1+delta2;

%%d_delta为误差分布的间距

d_delta=(max(delta)-min(delta))/(M-1);

仪器精度理论实验报告

delta_n=[min(delta):d_delta:max(delta)]; %% 画直方图 %%正弦分布的直方图 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([-2,2,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(1*abs(H))/N; %% 计算直方图的面积%%

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积 i=1:1:M/2;

s_1(i)=1*abs(floor(H(floor(M/2+i))))/N; s_2(i)=1*abs(floor(H(floor(M/2-i+1))))/N; 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

仪器精度理论实验报告

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/2+iii))-delta_n(floor(M/2-iii+1)))/6; %% 计算最短包含区间面积 %s1为包含q 个采样点的概率和 P=0.9973; %置信概率

q=P*M; %q为最短包含区间应该包含的采样点个数 %求置信概率为99.73%时的所有包含区间 for j=1:M;

for jj=j:M;

area1(jj)=sum(H(j:jj))*M/N; if area1(jj)>q result(j)=jj; break; else ; end end end

%求所有包含区间的长度 for k=1:length(result);

result1(k)=result(k)-k; end

%求最短包含区间的长度及起始区间r [q,r]=min(result1); delta_m=q*d_delta;

%%理论计算标准不确定度%% delta_mean=mean(delta);

delta_cancha=delta-delta_mean;

s=sqrt((sum(delta_cancha.^2))/(N-1));

%%计算最短包含区间与MCM 法、GUM 法所得区间之差

仪器精度理论实验报告

delta_u2=abs(delta_m-6*s); toc;

程序运行后得到的合成误差频率分布图如图4-1所示，由图可知其频率分布为非正态分布，用最短包含区间法得到最短区间为[1654，5795]，即r=1654，q=9973。再与普通MCM 法、GUM 法的包含区间比较，分别得到两两差值为delta_u1= 0.3895，delta_u2= 0.7660。

-6

-4-20246

图4-1 非正态分布的频率分布

4.3 问题讨论

(1) 最短包含区间MCM 法与普通MCM 法有何区别。最短包含区间MCM 适合解决怎样的问题。

答：当PDF 对称时，最短包含区间MCM 法比普通MCM 法算得的包含区间相同；当PDF 不对称时，最短包含区间MCM 法算得的包含区间更小。最短包含区间MCM 法需要对所有满足置信概率的区间循迹，并找到区间的最小值，需要算多个区间；普通MCM 法只需从中心点向两侧计算得到满足置信概率的区间即可，只需算一次包含区间。最短包含区间MCM 法适用于PDF 不对称的情况，以获取输出量的最佳估计值、不确定度和包含区间。

(2) 最短包含区间法与GUM 法同时计算不确定度的时候，哪一个略显保守。 答：最短包含区间法得到的包含区间比GUM 法小，GUM 法略显保守。

仪器精度理论实验报告

仪器精度理论实验报告

学 院： 专 业： 姓 名： 学 号：

仪器精度理论实验报告

1 Monte Carlo模拟误差分析原理

在误差分析的过程中，常用的方法是通过测量方程推导出误差传递方程，再通过不确定度的合成公式获得间接测量量的标准不确定度和扩展不确定度，这种方法称为GUM 法。GUM 法适用于误差是线性的，然而在有些场合下，测量方程较难获得，在这种情况下研究误差的特性就需要借助于模拟统计的方式进行计算。Monte Carlo(MCM)法就是较为常用的数学工具。

两种方法的不确定度统计模拟分析方法如下： 设有测量方程

y =f (x 1, x 2,..., x n ) (1-1)

以计算机模拟m 次测量的各误差分量δx 1j , δx 2j , …, δx nj , j =1,2,…, m 。则测量数据可写成

x 1j =X 1+δx

j 1j 2

x 2j =X 2+δx x nj =X n +δx nj

j =1,2,…, m (1-2)

式中，X 1, X 2, , X n 为各项分量的真值。将各项数据带入测量方程，得

y j =f (x 1j , x 2j , x nj )

=f (X 1+δx 1j , X 2+δx 2j , , X n +δx nj )

(1-3)

该结果也可写为

y j =Y +δy j

=f (X 1+δx 1j , X 2+δx 2j , , X n +δx nj )

(1-4)

式中，Y 为被测量的真值，y j 为m 个测量数据，j =1,2,…, m ; δy j 为由各项误差分量

δx 1j , δx 2j , …, δx nj , j =1,2,…, m 造成的测量结果y j 的误差。

按GUM 法，由m 个测量数据y j ，取其算数平均值得

1m

y =∑y j (1-5)

m j =1

计算其残差

1m

v j =y j -y =y j -∑y j (1-6)

m j =1

由Bessel 公式得测量数据的标准差为

仪器精度理论实验报告

按MCM 法，将y j 以适宜的子区间间隔绘制成直方图，得到频率分布。对于置信概率p ，若pm 为整数，设q =pm ，否则取q 的值为pm +2的整数部分。则⎡⎣y low , y high ⎤⎦为Y 的其中，对任意的r =1, , m -q ，y low =y (r ) ，y high =y (r +q ) 。对于对称的PDF ，100p %包含区间，

如GUM 法中使用的正态分布和缩放位移t 分布，概率对称100p%包含区间和最短100p%包含区间是相同的。如果PDF 不对称，可采用最短100p %包含区间。确定r *，使得

y (r *+q ) -y (r *) ≤y (r +q ) -y (r ) ，r =1, , M -q ，可得最短100p %包含区间⎡y (r ) , y (r +q ) ⎤。

⎣

*

*

⎦

在执行自适应蒙特卡洛方法的基本过程中，蒙特卡洛试验次数不断增加，直至所需要的各种结果达到统计意义上的稳定。如果某结果的两倍标准偏差小于标准不确定度的数值容差时，则认定该数值结果稳定。

用MATLAB 可以在短时间内获得大量的误差分量样本，可以很好的模拟统计实验。本次实验包含三个小实验：

（1）MCM 法与GUM 法的实现及比较； （2）自适应MCM 法的实现；

（3）基于最短包含区间的MCM 法的实现。

仪器精度理论实验报告

2 MCM法与GUM 法的比较实验

2.1 实验步骤

GUM 法与MCM 法的实现步骤如下：

(1) 利用MATLAB 的rand 函数生成[0，1]区间内的均匀分布的随机数ξ1和ξ2； (2) 给出误差分量的随机值：利用MATLAB ，由均匀分布随机数ξ1生成标准正态分布随机数η1，误差分量随机数可表示为

δ1=η1u 1=5η1mV (2-1)

同理得

δ2=η2u 2=7η2mV (2-2)

(3) 求和的随机数：误差和的随机数δ=δ1+δ2；

(4) 重复以上步骤，得误差和的随机数系列：δi =δ1i +δ2i i =1, 2, n ；

(5) 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m 组（m 尽可能大），每组间隔为∆，记数各区间的随机数的数目n j ，以∆为底，以

n j n ∆

为高作第j (j =1, 2 m ) 区间的矩形，最终的m 组矩形构

∑n

成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。以频率

j =1

k

j

n

=99.73%为

界划定区间，该区间半宽即为测量总误差的置信概率为99.73%的扩展不确定度。

(6) 合成的标准不确定度：总误差随机数平均值为

1n

δ=∑δi (2-3)

n i =1

各误差随机数的残差

v i =δi - (2-4)

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

u =s =

∑v

i =1

n

2i

n -1

(2-5)

实验流程图如图2-1所示：

仪器精度理论实验报告

图2-1 实验流程图

2.2 程序实现

用MATLAB 实现如下： tic;

clear;clc;close all; bdclose all; %%设定参数值%%

%%随机信号点数N ，均值为1，标准差u1,u2%% N=10^6; M=N/100; x=0:1:M; x_=[1:M]; u1=0.005; u2=0.007;

%%产生两个在(0,1)上服从均匀分布的, 种子为53, 每一次都相同的随机数X1和X2 rand('state',53); 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;

仪器精度理论实验报告

lta 为误差分布的间距

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;

%%正态分布误差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(abs(H))/N;

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积

仪器精度理论实验报告

s_1(i)=1*abs(floor(H(floor(M/2+i))))/N; s_2(i)=1*abs(floor(H(floor(M/2-i+1))))/N; 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_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;

种子设定为学号后两位53，程序运行后，得到合成误差delta 概率分布的三组曲线： (1) N=106，M=104； (2) N=106，M=103； (3) N=105，M=104。

分别如图2-2(a)，2-2(b)，2-2(c)所示。

仪器精度理论实验报告

-0.05

-0.04

-0.03

-0.02

-0.01

0.01

0.02

0.03

0.04

0.05

图2-2(a) N=10, M=10delta 概率分布

6

4

-0.05

-0.04-0.03-0.02-0.0100.010.020.030.040.05

图-2(b) N=10, M=10delta 概率分布

6

3

仪器精度理论实验报告

-0.04

-0.03-0.02-0.0100.010.020.030.04

图2-2(c) N=10, M=10delta 概率分布

5

4

三种情况下，算出的GUM 法和MCM 法得到的不确定度分别为： (1) 图(a)，delta_n_u =0.0084，s =0.0086； (2) 图(b)，delta_n_u =0.0073，s =0.0086； (3) 图(c)，delta_n_u =0.0084，s =0.0086。

2.3 问题讨论

(1) N(采样点数) ，M(划分的区间数) 与直方图的关系(平滑，Y 轴的高度) 。

答：通过对比图(a)与图(b)可以看出：划分的区间M 越小，Y 轴高度越高，曲线越平滑。因为当采样点数N 不变时，区间数M 减少导致单个区间的频数变大，所以区间的过渡更接近理想正态分布曲线。同时，通过对比图(a)与图(c)可以看出：采样点数N 越大，Y 轴高度越高，曲线越平滑。因为当划分的区间数M 不变时，采样点数N 增加导致单个区间的频数变大，所以区间的过渡更接近理想正态分布曲线。

(2)Bessel公式计算的标准不确定度(GUM)与99.73%直方图面积的扩展不确定度(MCM)两者之间会存在误差，这个误差与哪些因素有关(N,M,N与M 的关系) 。

答：通过对比由(1)和(2)得到的参数可知：取N=106，M=104时，u 和s 的误差不大；取N=106，M=103时，u 和s 的误差很大。再通过对比由(1)和(3)得到的参数可知：取N=106，M=104时，u 和s 的误差不大；取N=105，M=104时，u 和s 的误差更小。综上说明：采样点数N 对误差的影响不大，划分区间数M 对误差的影响很大，且N 与M 的比值越大，误差越大。

仪器精度理论实验报告

3 自适应MCM 法实验

3.1 实验步骤

自适应MCM 法的步骤如下：

(1) 基于前一个实验，构建衡量Monte Carlo 法和GUM 法计算得到标准不确定度差值的函数。

(2) 将随机数个数N ，分割区间数M 分别作为该函数的自变量，定义自变量的取值范围，从而获得相应的函数值。

(3) 分别进行三维网格作图和三维曲线作图，通过观察曲线获得最佳的N ，M 组合。

3.2 程序实现

用MATLAB 实现如下： tic;

close all;clear;clc;

%根据三维网格曲线求最佳的N ，M 组合 %a,b模拟N 和M a=logspace(1,6); b=logspace(1,6);

%j代表N ，jj 代表M ，且N>M for j=1:max(size(a)); for jj=1:j;

Result1(j,jj)=shiyan(a(j),b(jj)); end end

%求Result1的非零最小值，并返回行号和列号 Result1(find(Result1==0))=NaN; f=min(min(Result1)); [n,m]=find(Result1==f);

[a,b]=meshgrid(logspace(1,6)); %求对应行列的N 、M M=a(n,m); N=b(n,m); %作三维网格曲线 figure(1),surfc(a,b,Result1); xlabel('M');

仪器精度理论实验报告

zlabel('Δ');

%根据三维曲线求最佳的N ，M 组合 %c、d 模拟N 、M

c=logspace(1,6); d=logspace(1,6);

for jjj=1:max(size(c));

Result2(jjj)=shiyan(c(jjj),d(jjj)); end

%求Result2的最小值，返回最小值和列号 [Y,nm]=min(Result2); NM=c(nm); %作三维曲线

figure(2),plot3(c,d,Result2);grid on; xlabel('M'); ylabel('N'); zlabel('Δ'); toc;

其中函数shiyan 的代码如下，以N 、M 为参量，以GUM 法与MCM 法所得不确定度之差作为返回值，种子数为53。 function y=shiyan(N,M) %%设定参数值%%

%%随机信号点数N ，均值为1，标准差u1,u2 Mu=1;

x=0:1:floor(M); x_=[1:floor(M)]; u1=0.005; u2=0.007;

%%产生两个在(0,1)上服从均匀分布的, 种子为53, 每一次都相同的随机数X1和X2%% rand('state',53); 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;

仪器精度理论实验报告

delta=delta1+delta2;

d_delta=(max(delta)-min(delta))./(floor(M)-1); %%d_delta为误差分布的间距 delta_n=[min(delta):d_delta:max(delta)]; %%delta_n为误差分布序列 H=hist(delta,floor(M)); %%计算直方图的面积%% S=sum(1.*abs(H))/N; %% 计算直方图的面积%%

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积 i=1:1:floor(M./2);

s_1(i)=1*abs(floor(H(floor(M./2+i))))/N; s_2(i)=1.*abs(floor(H(floor(M./2-i+1))))/N; 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

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)); %%计算GUM 法与MCM 法所得不确定度的差值 y=abs(delta_n_u-s);

仪器精度理论实验报告

出Δ的最小值及对应的M 、N 值，如图3-1(a)所示。

Δ

M

5

N

图3-1(a) 三维网格曲线图

三维曲线如图3-1(b)所示，X 轴和Y 轴分别为M 、N ，Z 轴为Δ=abs(delta_n_u-s)，并得到Δ的最小值及M 、N 。

Δ

M

5

N

图3-1(b) 三维曲线图

3.3 问题讨论

如何根据三维网格曲线和三维曲线获得最佳的N ，M 组合。

仪器精度理论实验报告

线，其误差Δ的最小值为3.6949×10-8，在矩阵中处于43行49列，对应的M 和N 分别为M=1.9307×105，N=7.9060×105。对于三维曲线，其误差Δ的最小值为1.6837×10-6，处于41行41列，对应的M 和N 为1.2068×105。

仪器精度理论实验报告

4 基于最短包含区间的MCM 法实验

4.1 实验步骤

基于最短包含区间的MCM 法的实现步骤如下： (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*对应的区间长度极为最短包含区间。

4.2 程序实现

用MATLAB 实现如下： tic;

clear;clc;close all; bdclose all; %%设定参数值%%

%%随机信号点数N ，均值1，标准差%% N=10^6; M=N/100; x=0:1:M; x_=[1:M]; u1=1; u2=1;

%%产生两个在(-1,1)内的均匀分布, 种子为53, 每一次都相同的随机数X1和X2 rand('state',53);

X1=unifrnd(-1,1,1,N); X2=unifrnd(-1,1,1,N);

%%产生正弦分布Y1和正态分布Y2 Y1=sin(X1);

Y2=sqrt(-2*log(abs(X1))).*cos(2*pi*X2); %% 为做直方图先定义好X 轴的坐标数据%% delta1=u1*Y1; delta2=u2*Y2;

delta=delta1+delta2;

%%d_delta为误差分布的间距

d_delta=(max(delta)-min(delta))/(M-1);

仪器精度理论实验报告

delta_n=[min(delta):d_delta:max(delta)]; %% 画直方图 %%正弦分布的直方图 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([-2,2,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(1*abs(H))/N; %% 计算直方图的面积%%

%%s_1表示正向直方图的每一个单元的面积 %%s_2表示反向直方图的每一个单元的面积 %%s_表示正反向两两对称每一对单元的面积 %%area表示以中心为对称轴的累加面积 i=1:1:M/2;

s_1(i)=1*abs(floor(H(floor(M/2+i))))/N; s_2(i)=1*abs(floor(H(floor(M/2-i+1))))/N; 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

仪器精度理论实验报告

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/2+iii))-delta_n(floor(M/2-iii+1)))/6; %% 计算最短包含区间面积 %s1为包含q 个采样点的概率和 P=0.9973; %置信概率

q=P*M; %q为最短包含区间应该包含的采样点个数 %求置信概率为99.73%时的所有包含区间 for j=1:M;

for jj=j:M;

area1(jj)=sum(H(j:jj))*M/N; if area1(jj)>q result(j)=jj; break; else ; end end end

%求所有包含区间的长度 for k=1:length(result);

result1(k)=result(k)-k; end

%求最短包含区间的长度及起始区间r [q,r]=min(result1); delta_m=q*d_delta;

%%理论计算标准不确定度%% delta_mean=mean(delta);

delta_cancha=delta-delta_mean;

s=sqrt((sum(delta_cancha.^2))/(N-1));

%%计算最短包含区间与MCM 法、GUM 法所得区间之差

仪器精度理论实验报告

delta_u2=abs(delta_m-6*s); toc;

程序运行后得到的合成误差频率分布图如图4-1所示，由图可知其频率分布为非正态分布，用最短包含区间法得到最短区间为[1654，5795]，即r=1654，q=9973。再与普通MCM 法、GUM 法的包含区间比较，分别得到两两差值为delta_u1= 0.3895，delta_u2= 0.7660。

-6

-4-20246

图4-1 非正态分布的频率分布

4.3 问题讨论

(1) 最短包含区间MCM 法与普通MCM 法有何区别。最短包含区间MCM 适合解决怎样的问题。

答：当PDF 对称时，最短包含区间MCM 法比普通MCM 法算得的包含区间相同；当PDF 不对称时，最短包含区间MCM 法算得的包含区间更小。最短包含区间MCM 法需要对所有满足置信概率的区间循迹，并找到区间的最小值，需要算多个区间；普通MCM 法只需从中心点向两侧计算得到满足置信概率的区间即可，只需算一次包含区间。最短包含区间MCM 法适用于PDF 不对称的情况，以获取输出量的最佳估计值、不确定度和包含区间。

(2) 最短包含区间法与GUM 法同时计算不确定度的时候，哪一个略显保守。 答：最短包含区间法得到的包含区间比GUM 法小，GUM 法略显保守。