[误差理论与数据处理]使用Python与常见判别准则判别处理粗大误差

最近在学这门课，正好讲到粗大误差的处理，然后我发现这个处理实在难绷，数据稍微多一点计算器算残差和标准差就麻烦的要死，正巧有自学过部分Python内容，于是用Python实现了一遍作业里面要求的和书上列出过的方法，现在把源代码放在这里供后面学这门课的同志们参考。本程序的特点是可以一路生成完整的计算过程……打印的计算过程也算计算过程。

误差处理部分是有国家标准的，标准号GB/T 4883-2008，查表的时候查不到可以到这找。

第零部分：包引用

首先引用几个可能有点用但是实际上没啥大用处的包：

scipy.stats.t.interval：可以生成T分布表，我不想输入了所以引用一下。这个函数返回值是俩，一个正一个负，用的都是正值，而且书里的alpha和函数的alpha加起来是1（要反一下）。

statistics.stdev/fmain:stdev是标准差，fmean是快速浮点平均值，我不想做实现了而且这个模块是自带的所以干脆用一下

prettytable：打印残差的表格，懒得弄，用一下

bpr()函数：因为经常需要打印一串字符串+一串数字，干脆包装一下，顺便统一浮点数取整

from scipy.stats import t
from statistics import stdev,fmean
import prettytable as pt

def bpr(word,data):
    print(f"{word}{round(data,2)}")

第一部分：3σ准则（莱以特准则）

这个部分应该是最简单的部分，毕竟其实就一个标准差一个平均值就没了。stdev做的那个标准差是样本标准差（pstdev是全体样本标准差），下面除的是n-1，问题应该不大。实在不行自己实现一个就是了

def ThreeSigmaMethod(data,alpha=None):
    if len(data)<10:
        raise Warning("数据个数小于10个时，3σ原则失效。")
    
    #计算平均值和标准差
    avg=fmean(data);stderr=stdev(data)
    
    #虽然我们用不上，但是过程要求计算残差和残差平方
    vi=[x-avg for x in data]
    vi2=[x**2 for x in vi]
    
    #判断是否存在粗大误差
    #输出初始化：[正常数据],[粗大误差]
    out=[[],[]]
    for x in data:
        if not abs(x-avg)<=3*stderr:
            out[1].append(x)
        else:
            out[0].append(x)
    
    
    #过程生成部分
    bpr("数据平均值：",avg)
    tb=pt.PrettyTable(["原始数据","残差","残差平方"])
    for i in range(len(data)):
        tb.add_row([data[i],round(vi[i],6),round(vi2[i],6)])
    print(tb)
    
    return out

我自己写东西的时候不大喜欢开长行标记（45字符那个），而且编码习惯也不大好，别学我

这个方法输入一组data，输出一个嵌套了的列表，第一个元素是一个这样的列表[[非粗大误差数据],[粗大误差数据]]，后面几个元素是平均值、标准差和残差列表（未经平方）

第二部分：罗曼诺夫斯基准则

这个查的好像是T分布表，我就干脆拿scipy生成了

不想生成也可以自己把表打一个出来。

def RomanowskiMethod(data,alpha=0.95):
    #计算平均值、残差、残差平方
    avg=fmean(data)
    vi=[x-avg for x in data]
    vi2=[x**2 for x in vi]
    
    #保存原始数据
    sdata=data[:]
    
    #根据方法要求，直接删除可疑目标后再判断
    sus=data.pop(vi2.index(max(vi2)))
    
    #计算新的平均值和标准差
    avgn=fmean(data)
    stderr=stdev(data)
    susflag=False
    #根据方法，判断刚刚删除的是否为粗大误差
    if abs(sus-avgn)>round(abs(t.interval(alpha,len(data)+1)[1]),2)*stderr: #this dirty "if" is the actual core of this func.
        susflag=True

    
    #生成“计算过程”
    bpr("数据平均值：",avg)
    tb=pt.PrettyTable(["原始数据","残差","残差平方"])
    for i in range(len(sdata)):
        tb.add_row([sdata[i],round(vi[i],6),round(vi2[i],6)])
    print(tb)
    bpr("可疑数据项：",sus)
    bpr("排除可疑数据后的平均值：",avgn)
    bpr("排除可疑数据后的标准差：",stderr)
    if susflag:
        print("经验证，数据确实为粗大误差。")
    else:
        print("经验证，数据并非粗大误差")
    
    
    #返回结果
    if susflag:
        return [data,[sus]]
    data.append(sus)
    return [data,[]]

第三部分：格罗布斯准则

这俩方法都要查表，数据基本都是我从GB/T 4883-2008上手动录下来的，这个表真长啊

另外国标管这个叫格拉布斯准则

def GrubbsMethod(data,alpha=0.95):
    
    def _getGrubbsMethod(datascale,alpha):
        datascale-=3
        if alpha==0.95:
            return [1.15,1.46,1.67,1.82,1.94,2.03,2.11,\
            2.18,2.23,2.28,2.33,2.37,2.41,2.44,2.48,2.50,\
            2.53,2.56,2.58,2.60,2.62,2.64,2.66,2.74,2.81,2.87,2.96,3.17][datascale]
        if alpha==0.99:
            return [1.16,1.49,1.75,1.94,2.10,2.22,2.32,\
            2.41,2.48,2.55,2.61,2.66,2.70,2.75,2.78,2.82,\
            2.85,2.88,2.91,2.94,2.96,2.99,3.01,3.10,3.18,3.24,3.34,3.59][datascale]
    if len(data)<3:
        #表里没有
        raise ValueError("Data scale must be above 3.")
    if alpha not in [0.90,0.95,0.99,0.995]:
        #表里确实没有
        raise ValueError(f"Alpha {alpha} doesn't match any of the one listed : 0.95,0.99.")
    
    #排序数据
    sdata=data[:]
    sdata.sort()
    
    #计算平均值，标准差
    avg=fmean(data);stderr=stdev(data)
    
    #输出过程
    bpr("数据平均值：",avg)
    bpr("数据标准差：",stderr)
    bpr("数据最大值：",sdata[-1])
    bpr("数据最小值：",sdata[0])
    
    #获取当前阈值
    _GrubbusThreshold=_getGrubbsMethod(len(data),alpha)
    bpr("当前方法判断极限阈值：",_GrubbusThreshold)
    
    if avg-sdata[0]>=sdata[-1]-avg:
        #x(1)更值得怀疑的情况
        print("x(1)更值得怀疑")
        g1=(avg-sdata[0])/stderr
        if g1>_GrubbusThreshold:
            sus=sdata.pop(0)
            print(f"校核发现x(1):{sus}为粗大误差")
            return [sdata,[sus]]
        else:
            print("校核未发现粗大误差")
            return [sdata,[]]
    else:
        #x(n)更值得怀疑的情况
        print(f"x({len(data)})更值得怀疑")
        gn=(sdata[-1]-avg)/stderr
        if gn>_GrubbusThreshold:
            sus=sdata.pop(-1)
            print(f"校核发现x({len(data)}):{sus}为粗大误差")
            return [sdata,[sus]]
        else:
            print("校核未发现粗大误差")
            return [sdata,[]]

第四部分：狄克松准则

同上，手动输入的单边狄克松临界值表，我自己看了一下应该是没啥问题。GB/T 4883-2008也没给这个表咋算出来的，所以不能处理所有的alpha置信度情况，只有0.90/0.95/0.99/0.995四组，而且大于30的太长了我也没录入，凑合着过吧。

注意这里的alpha是没有用1减的，0.90就是90%，不是书上那个0.01（99%）和0.05（95%）凹

def DixonMethod(data,alpha=0.95):
    
    #查数据表动作的函数实现
    def _getDixonValue(datascale,alpha):
        datascale=-3
        if alpha==0.90:
            return [0.885,0.679,0.557,0.484,0.434,0.479,0.441,0.410,0.517,\
            0.490,0.467,0.491,0.470,0.453,0.437,0.424,0.412,0.401,0.391,\
            0.382,0.374,0.367,0.360,0.353,0.347,0.341,0.337,0.332][datascale]
        if alpha==0.95:
            return [0.941,0.765,0.642,0.562,0.507,0.554,0.512,0.477,\
            0.575,0.546,0.521,0.546,0.524,0.505,0.489,0.475,0.462,\
            0.450,0.440,0.431,0.422,0.413,0.406,0.399,0.393,0.387,0.381,0.376][datascale]
        if alpha==0.99:
            return [0.988,0.889,0.782,0.698,0.637,0.681,0.635,0.597\
            ,0.674,0.642,0.617,0.640,0.618,0.597,0.580,0.564,0.550,\
            0.538,0.526,0.516,0.507,0.497,0.489,0.482,0.474,0.468,0.462,0.456][datascale]
        if alpha==0.995:
            return [0.994,0.920,0.823,0.744,0.680,0.723,0.676,0.638,\
            0.707,0.675,0.649,0.672,0.649,0.629,0.611,0.595,0.580,\
            0.568,0.556,0.545,0.536,0.526,0.519,0.510,0.503,0.496,0.489,0.484][datascale]
    
    if len(data)<3:
        #表里没有
        raise ValueError("Data scale must be above 3.")
    if alpha not in [0.90,0.95,0.99,0.995]:
        #表里确实没有
        raise ValueError(f"Alpha {alpha} doesn't match any of the one listed in GB/T 4883-2008: 0.90,0.95,0.99,0.995.")
    
    #排序统计量
    sdata=data[:]
    sdata.sort()
    
    #统计量全算一遍
    Hr10=(sdata[-1]-sdata[-2])/(sdata[-1]-sdata[0])
    Hr11=(sdata[-1]-sdata[-2])/(sdata[-1]-sdata[1])
    Hr21=(sdata[-1]-sdata[-3])/(sdata[-1]-sdata[1])
    Hr22=(sdata[-1]-sdata[-3])/(sdata[-1]-sdata[2])
    Lr10=(sdata[0]-sdata[1])/(sdata[0]-sdata[-1])
    Lr11=(sdata[0]-sdata[1])/(sdata[0]-sdata[-2])
    Lr21=(sdata[0]-sdata[2])/(sdata[0]-sdata[-2])
    Lr22=(sdata[0]-sdata[2])/(sdata[0]-sdata[-3])
    
    #根据数据规模选择对应的统计量
    if len(data)<7:
        CurDH=Hr10;CurDL=Lr10
    elif len(data)<10:
        CurDH=Hr11;CurDL=Lr11
    elif len(data)<13:
        CurDH=Hr21;CurDL=Lr21
    elif len(data)<30:
        CurDH=Hr22;CurDL=Lr22
    else:
        raise ValueError(f"Data Scale {len(data)} is currently too big for this function. Consult GB/T 4883-2008 for further information.")
    out=[sdata[:],[]]
    _DixonThreshold=_getDixonValue(len(sdata),alpha)
    if CurDH>_DixonThreshold:
        out[1].append(out[0].pop(-1))
    if CurDL>_DixonThreshold:
        out[1].append(out[0].pop(0))
    
    bpr("当前使用的统计量（更大的）：",CurDH)
    bpr("当前使用的统计量（更小的）：",CurDL)
    bpr("当前设置下的极限值：",_DixonThreshold)
    if len(out[1])==0:
        print("经校核未发现问题。")
    else:
        bpr("经校核发现粗大误差：",out[1][0])
    return out

这个函数输入是一组data和一个置信度alpha，输出还是一个嵌套列表，[[非粗大误差数据],[粗大误差数据]]

第五部分：通用处理方法包装器

其实就是循环使用全部方法直到输出的结果中没有粗大误差数据

def globalMethodWrapper(methodname,data,alpha=0.95):
    CurMethod=methodname
    CurData=data[:]
    SusData=[]
    i=1
    
    while True:
        print(f"开始第{i}轮粗大误差判别：")
        out=CurMethod(CurData,alpha)
        if len(out[1])==0:
            print("数据已不含有粗大误差，结束。")
            return
        else:
            SusData.extend(out[1])
            CurData=out[0][:]
            print(f"当前判别的粗大误差为：{out[1][0]}")
            print(f"已经判别的所有粗大误差为：{'  '.join([str(x) for x in SusData])}")
            print()

这个包装器的输入直接用函数名就行（不要带括号进去），使用例：

globalMethodWrapper(ThreeSigmaMethod,Data)

至于前面的3sigma原则也加了个alpha是为了匹配后面几个有alpha的