最近在学这门课,正好讲到粗大误差的处理,然后我发现这个处理实在难绷,数据稍微多一点计算器算残差和标准差就麻烦的要死,正巧有自学过部分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的