m_grubbs.m

一组仪器性能分析程序

function [m,s,n,id]=M_grubbs(x,xid,n,id)
%输入数组x,xid，粗大误差点的初始个数n=0，初始id=0（必须输入0或值为0的变量）
%输出根据Grubbs准则剔除粗大误差后的平均值和标准差以及被剔除的粗大误差点的id和个数

%如果在此基础上修改，输出部分添加x,xid则得到剔除粗大误差后的数组和对应编号
L=length(x)
if L~=length(xid)           
    return                              %如果数组和对应id长度不同则退出
end
m=mean(x)                               %平均值
s=std(x)                                %标准差

%figure
%hold on
%plot(xid,x,'*k')
%plot([xid(1),xid(L)],[m+3*s,m+3*s],'-b')
%plot([xid(1),xid(L)],[m-3*s,m-3*s],'-b')
%hold off

flag=false
for i=1:1:L
    if abs(x(i)-m)>3*s      
        flag=true                       %当找到粗大误差点则跳出循环，这时i为该点的下标
        break
    end
end
if flag 
    n=n+1                               %粗大误差点个数+1
    id(n)=xid(i)                        %粗大误差点的编号记录入数组id
    x=[x(1:i-1) x(i+1:L)]               %剔除粗大误差点（已验证能剔除任意点）
    xid=[xid(1:i-1) xid(i+1:L)]
    [m,s,n,id]=M_grubbs(x,xid,n,id)           %递归调用Grubbs以在新数组中剔除粗大误差
end
%这里递归调用Grubbs函数有两点
%（1）位置放最后，以免从最里层返回时进入上一层继续计算，有可能发生错误（比如数组下标越界）
%（2）调用时需要同步返回值，否则返回值保持最顶层的值