|
from http://blog.csdn.net/zhzht19861011/article/details/4620896 程序实现的功能：
disp('#################################################################')
disp('## 基本数据处理程序 ##')
disp('####请选择： ####')
disp('## 0:求平均值 ##')
disp('## 1:求残差 ##')
disp('## 2:求标准差 ##')
disp('## 3:算术平均值的标准差 ##')
disp('## 4:极限误差 ##')
disp('## 5:判断粗大误差 ##')
disp('## 6:判断系统误差 ##')
disp('## 7:一键分析 ##')
disp('## 8:退出 ##')
disp('##################################################################')
源程序：
clear
clc
%以下为罗曼诺夫斯基准则中的判别表 a=0.05
cd_va=[4.97 3.56 3.30 2.78 2.62 2.51 2.43 2.37 2.33 2.29 2.26 2.24 2.22 2.20 2.18 2.17 2.26 2.15 2.14 2.13 2.12 2.11 2.10 2.10 2.09 2.09 2.08];
%以下为t分布表，a=0.01
t=[4.03 3.71 3.50 3.36 3.25 3.17 3.11 3.05 3.01 2.98 2.95 2.92 2.90 2.88 2.86];
%以下为检验周期性系统误差需要的表 p=0.95
zhouqi=[0.39 0.41 0.445 0.468 0.492 0.512 0.531 0.548 0.564 0.578 0.591 0.603 0.614 0.624];
%以下数据为罗曼诺夫斯基判别粗大误差所用数据 取a=0.01
luo=[11.46 6.53 5.04 4.36 3.96 3.71 3.53 3.41 3.31 3.23 3.17 3.12 3.08 3.04 3.01 3.00 2.95 2.93];
%=0; %显著度
while 1
%x=[20.42 20.43 20.40 20.43 20.42 20.43 20.39 20.30 20.40 20.43 20.42 20.41 20.39 20.39 20.40];
disp('请在下面输入测量列，格式如右所示：[n1 n2 n3 n4 ... nx]')
x=input('')
l=length(x);
mean_x=0; %平均值
for i=1:l
mean_x=mean_x+x(1,i)/l;
end
v=0;
for i=1:l
v(1,i)=x(1,i)-mean_x; %v为残差
end
y=0; %中间变量
for i=1:l
y=y+v(1,i)*v(1,i)/(l-1);
end
q=sqrt(y); % q为标准差
qx=q/sqrt(l); %算术平均的标准差
dx=qx*t(1,(l-4)); %极限误差
%------------------------------------------------------
%以下用罗曼夫斯基准则判别粗大误差
%------------------------------------------------------
for i=1:l %对残差求绝对值
s(1,i)=abs(v(1,i));
end
max_=max(s); %求残差绝对值中的最大值
n=1; %改粗大误差值对应的列号
for i=1:l
if s(1,i)==max_
break;
else
n=n+1;
end
end
qq=1; %剔除算法
pp=1;
for i=1:l
if i==n
pp=pp+1;
else
x1(1,qq)=x(1,pp); %剩余测量列
qq=qq+1;
pp=pp+1;
end
end
%计算剩余测量列的平均值、标准差、等
len=length(x1);
mean_x1=0; %平均值
for i=1:len
mean_x1=mean_x1+x1(1,i)/len;
end
v1=0;
for i=1:len
v1(1,i)=x1(1,i)-mean_x1; %v1为残差
end
y1=0; %中间变量
for i=1:len
y1=y1+v1(1,i)*v1(1,i)/(len-1);
end
q1=sqrt(y1); % q1为标准差
qx1=q1/sqrt(len); %算术平均的标准差
dx1=qx1*t(1,(len-4)); %极限误差
kl=luo(1,(l-3));
kl=kl*q1;
tl=x(1,n)-x1;
tl=abs(tl);
if tl>kl
lyt_f=0;
else
lyt_f=1;
end
for i=1:len %对残差求绝对值
s1(1,i)=abs(v1(1,i));
end
max_2=max(s1); %求残差绝对值中的最大值
n2=1;
for i=1:len
if s1(1,i)==max_2
break;
else
n2=n2+1;
end
end
% 第二次剔除可疑值
qq1=1;
pp1=1;
for i=1:len
if i==n2
pp1=pp1+1;
else
h1(1,qq1)=x1(1,pp1);
qq1=qq1+1;
pp1=pp1+1;
end
end
len2=length(x1);
mean_x2=0; %平均值
for i=1:len2
mean_x2=mean_x2+x1(1,i)/len2;
end
v2=0;
for i=1:len2
v2(1,i)=x1(1,i)-mean_x2; %v为残差
end
y2=0;
for i=1:len2
y2=y2+v2(1,i)*v2(1,i)/(len2-1);
end
q2=sqrt(y2); % q为标准差
k2=q2*luo(1,(l-3));
can1=x1(1,n2)-mean_x2;
can1=abs(can1);
if can1>k2
lyt_f1=0; %有粗大误差
else
lyt_f1=1;
end
%残余误差校核法校核系统误差
if lyt_f==1
for i=1:l
v2(1,i)=v(1,I);
end
else
for i=1:len
v2(1,i)=v1(1,i);
end
end
%校核线性误差
len1=length(v2);
k=(len1+1)/2;
k=fix(k); %截尾取整
sum1=0;
for i=1:k
sum1=sum1+v2(1,i);
end
sum2=0;
for i=(k+1):len1
sum2=sum2+v2(1,i);
end
sub=sum1-sum2;
sub=abs(sub);
if sub<0.1 %差值小于0.1断定含有线性误差,wucha=0
wucha=1;
else
wucha=0;
end
% 校核周期线性误差
bz=0;
for i=1:(len1-1)
bz=bz+(v2(1,i)-v2(1,(i+1)))*(v2(1,i)-v2(1,(i+1)));
end
ba=0;
for i=1:len1
bz=ba+v2(1,i)*v2(1,i);
end
xx1=bz/(2*ba);
xx2=zhouqi(1,(len1-3)); %xx1>xx2不存在周期性误差
while 1
clc
disp('#################################################################')
disp('## 基本数据处理程序 ##')
disp('####请选择： ####')
disp('## 0:求平均值 ##')
disp('## 1:求残差 ##')
disp('## 2:求标准差 ##')
disp('## 3:算术平均值的标准差 ##')
disp('## 4:极限误差 ##')
disp('## 5:判断粗大误差 ##')
disp('## 6:判断系统误差 ##')
disp('## 7:一键分析 ##')
disp('## 8:退出 ##')
disp('##################################################################')
disp('输入你的选项：')
option=input('')
clc
switch option
case 0
if lyt_f==1
disp('该测量列的平均值：')
mean_x
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列中存在粗大误差，是否剔除？ 0：剔除，1：不剔除')
pdf=input('')
if pdf==0
clc
disp('剔除粗大误差后的平均值为：')
mean_x1
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有粗大误差，你未剔除改粗大误差，未剔除粗大误差测量列的平均值为：')
mean_x
disp('按回车键继续……')
zhu=input('')
end
end
case 1
if lyt_f==1
disp('该测量列的残差：')
v
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列中存在粗大误差，是否剔除？ 0：剔除，1：不剔除')
pdf=input('')
clc
if pdf==0
disp('剔除粗大误差后的残差为：')
v1
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有粗大误差，你未剔除改粗大误差，未剔除粗大误差测量列的平均值为：')
v
disp('按回车键继续……')
zhu=input('')
end
end
case 2
if lyt_f==1
disp('该测量列的标准差：')
q
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列中存在粗大误差，是否剔除？ 0：剔除，1：不剔除')
pdf=input('')
clc
if pdf==0
disp('剔除粗大误差后的标准差为：')
q1
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有粗大误差，你未剔除改粗大误差，未剔除粗大误差测量列的标准差值为：')
q
disp('按回车键继续……')
zhu=input('')
end
end
case 3
if lyt_f==1
disp('该测量列算术平均值的标准差：')
qx
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列中存在粗大误差，是否剔除？ 0：剔除，1：不剔除')
pdf=input('')
clc
if pdf==0
disp('剔除粗大误差后算术平均值的标准差为：')
qx1
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有粗大误差，你未剔除改粗大误差，未剔除粗大误差测量列算术平均值的标准差为：')
qx
disp('按回车键继续……')
zhu=input('')
end
end
case 4
if lyt_f==1
disp('按t分布、a=0.01计算测量列的极限误差：')
dx
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列中存在粗大误差，是否剔除？ 0：剔除，1：不剔除')
pdf=input('')
clc
if pdf==0
disp('剔除粗大误差后，按t分布、a=0.01 计算测量列的极限误差为：')
dx1
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有粗大误差，你未剔除改粗大误差！！未剔除粗大误差的测量列，按t分布、a=0.01计算测量列的列极限误差为：')
dx
disp('按回车键继续……')
zhu=input('')
end
end
case 5
if lyt_f==1
disp('使用罗曼诺夫斯基判别法判别粗大误差，a=0.01')
disp('该测量列无粗大误差')
disp('按回车键继续……')
zhu=input('')
else
disp('使用罗曼诺夫斯基判别法判别粗大误差，a=0.01')
disp('该测量列中含有粗大误差')
disp('该粗大误差所在的测量列列值为：')
n
disp('该值为：')
x(1,n)
disp('是否剔除该值？ 0：剔除，1不剔除')
jj=input('')
clc
disp('再次使用罗曼诺夫斯基判别法判别剩余测量列的粗大误差，a=0.01')
disp('计算剔除后是否还有粗大误差……')
disp(' ')
disp('计算完毕：')
if lyt_f1==1
disp(' ')
disp('剔除后的测量列不含粗大误差')
disp('按回车键继续……')
zhu=input('')
else
disp(' ')
disp('使用罗曼诺夫斯基判别法判别粗大误差，a=0.01')
disp('剔除后的测量列仍含粗大误差')
disp('该粗大误差所在的测量列列值为：')
n1
disp('该值为：')
x1(1,n2)
disp('按回车键继续……')
zhu=input('')
end
end
case 6
disp('')
disp('使用残余误差校核法校核该测量列的系统误差')
if wucha==1
disp('该测量列含有线性系统误差')
else
disp('该测量列不含线性系统误差')
end
disp('')
disp('使用小样本序校核法校核该测量列的周期性系统误差')
if xx1>xx2
disp('')
disp('该测量列不含有周期性系统误差')
disp('按回车键继续……')
zhu=input('')
else
disp('该测量列含有周期性误差')
disp('按回车键继续……')
zhu=input('')
end
case 7
disp(['该测量列为',x ])
x
if lyt_f==1
disp('使用罗曼诺夫斯基判别法判别粗大误差，a=0.01')
disp(['该测量列无粗大误差'])
disp(['该测量列的算术平均值',mean_x])
mean_x
disp(['该测量列的标准差',q])
q
disp(['该测量列算术平均值的标准差',qx])
qx
disp(['该测量列的极限误差dx',dx])
else
disp('使用罗曼诺夫斯基判别法判别粗大误差，a=0.01')
disp('该测量列含有粗大误差')
disp('该粗大误差在的列值为：')
n
disp('该值为：')
x(1,n)
if lyt_f1==1
disp('将该值剔除……')
disp('再次使用罗曼诺夫斯基判别法判别剩余测量列的粗大误差，a=0.01')
disp('判别中……')
disp('判别完毕：')
disp('剔除粗大误差后的测量列无粗大误差')
else
disp('剔除粗大误差后的测量列仍有粗大误差')
end
disp(' ')
disp(['剔除粗大误差后的测量列的算术平均值',mean_x1])
mean_x1
disp(['剔除粗大误差后的测量列的标准差',q1])
q1
disp(['剔除粗大误差后的测量列算术平均值的标准差',qx1])
qx1
disp(['剔除粗大误差后的测量列的极限误差',dx1])
dx1
end
disp('使用残余误差校核法校核该测量列的线性系统误差:')
if wucha==1
disp('该测量列含有线性系统误差')
else
disp('该测量列不含线性系统误差')
end
disp(' ')
disp('使用小样本序校核法校核该测量列的周期性系统误差:')
if xx1>xx2
disp('该测量列不含有周期性系统误差')
else
disp('该测量列含有周期性系统误差')
end
disp('按回车键继续……')
zhu=input('')
case 8
break;
otherwise
disp('你的输入不符合规则，请重新输入。')
disp('按回车键继续……')
zhu=input('')
end
end
disp('确认退出程序？ 0：否，继续计算下组数据 / 1：是，退出程序')
flag=input('')
clc
if flag==1
break
end
end
|
