哈工大误差分析课程设计--Monte-Carlo

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

1 实验目的

1.1了解MATLAB 软件的基本功能和使用 1.2 学习不确定度的统计模拟分析方法

1.3 研究误差概率密度函数和Bessel 公式获得扩展不确定度的方法和影响因素

2实验原理

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

此次课程设计中按照实验要求产生的随机数可以模拟测量误差，通过对这些随机数的概率密度分布函数的面积、包络线和概率特征点的求取，可以获得随机误差的标准不确定度——(MCM)，并与理论上估计标准不确定度的Bessel 公式、极差法作——(GUM)比较，完成实验内容。并以此作为基础，分析GUM 法与MCM 法的区别与联系，影响MCM 法的参数，自适应MCM 法和基于最短包含区间的MCM 法。

已知两项误差分量服从正态分布，标准不确定度分别为51=u mV ， 72=u mV ，用统计模拟分析法给出两项误差和的分布(误差分布的统计直方图，合成的标准差，合成的置信概率 P 为99.73%的扩展不确定度)。

3实验内容

(1). 利用MATLAB 软件生成[0，1]区间的均匀分布的随机数ξ； (2). 给出误差分量的随机值：

利用MATLAB ，由均匀分布随机数1ξ生成标准正态分布随机数1η，误差分量随机数可表示为11115ηηδ==u mV ；22227ηηδ==u mV (3). 求和的随机数：误差和的随机数21δδδ+=；

(4). 重复以上步骤，得误差和的随机数系列：i i i 21δδδ+= n i ,2,1=； (5). 作误差和的统计直方图：以误差数值为横坐标，以频率为纵坐标作图。作图区间应包含所有数据，按数值将区间等分为m 组（m 尽可能大），每组间隔为∆，记数各区间的随机数的数目j n ，以∆为底，以

∆

n n j 为高作第j (m j 2,1=)区间的

矩形，最终构成误差和的分布直方图，该图包络线线即为实验的误差分布曲线。

(6). 以频率

%951

=∑=n

n

k

j j

为界划定区间，该区间半宽即为测量总误差的置信概率

为95%的扩展不确定度。

(7). 合成的标准不确定度：1

1

2

-=

=∑=n v

s u n

i i

4.实验流程图：

一．实验1

本实验中随机数种子为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=10^5;

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变换方法产生标准正态分布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;