有没有numpy内置函数可以从列表中排除异常值？

Question

是否有numpy内置函数可以执行以下操作？即，获取列表d并返回一个列表filtered_d，其中根据假定的点分布从d中删除任何异常值元素。

import numpy as np

def reject_outliers(data):
    m = 2
    u = np.mean(data)
    s = np.std(data)
    filtered = [e for e in data if (u - 2 * s < e < u + 2 * s)]
    return filtered

>>> d = [2,4,5,1,6,5,40]
>>> filtered_d = reject_outliers(d)
>>> print filtered_d
[2,4,5,1,6,5]

我说“类似于”的原因是该函数可能允许在不同分布（泊松分布、高斯分布等）和这些分布内的不同离群值阈值下进行操作（例如我在这里使用的m）。

- aaren

---

相关: scipy.stats能否识别和屏蔽明显的异常值？，尽管那个问题似乎处理更复杂的情况。对于您描述的简单任务，使用外部包似乎有些过度。 - Sven Marnach

我在想，考虑到主要的numpy库中内置函数的数量，没有类似的功能似乎有些奇怪。这似乎是处理原始、嘈杂数据的常见需求。 - aaren

线性异常值可以通过numpy std函数找到，但是如果数据是非线性的，例如二次函数或三次函数，标准差将无法很好地处理任务，因为它需要回归来帮助解决异常值。 - Weilory

这就是我编写这个仓库的原因：outliers.py。 - Weilory

Answer 1

处理异常值时，重要的一点是应该尽可能使用鲁棒性更强的估计器。分布的均值会受到异常值的影响，但例如中位数则会少得多。

在 eumiro 的回答基础上继续建设：

def reject_outliers(data, m = 2.):
    d = np.abs(data - np.median(data))
    mdev = np.median(d)
    s = d/mdev if mdev else np.zeros(len(d))
    return data[s<m]

在这里，我用更健壮的中位数替换了平均值，并用中位数绝对距离替换了标准差。然后，我通过它们（再次）的中位数值来缩放距离，以便m处于合理的相对比例尺上。

请注意，要使data[s<m]语法起作用，data必须是一个numpy数组。

- Benjamin Bannier

---

这基本上是所提到的修正Z-score，但门槛不同。如果我的数学没错，他们推荐 m 为 3.5 / .6745 ~= 5.189（他们将s乘以0.6745，并指定m为3.5...还要取abs(s)）。有人可以解释一下选择m的原因吗？或者这是您从特定数据集中确定的东西？ - Charlie G

@BenjaminBannier：你能否提供一些具体的解释来选择 m 的值，而不是像“纯度和效率的相互作用”这样的模糊陈述吗？ - stackoverflowuser2010

@stackoverflowuser2010：就像我所说的，这取决于您的具体要求，即我们需要将信号样本清洁到什么程度（假阳性），或者我们可以扔掉多少信号测量值来保持信号清洁（假阴性）。至于某个特定用例的具体示例评估，请参见例如http://www.desy.de/~blist/notes/whyeffpur.ps.gz。 - Benjamin Bannier

当我使用一个浮点数列表调用函数时，出现以下错误： “TypeError: 只有整数标量数组可以转换为标量索引” - Vasilis

@Charlie，如果您观看图表https://www.itl.nist.gov/div898/handbook/eda/section3/eda356.htm#MAD，您会发现在处理标准差为1的正态分布时（实际上这并不是需要使用改进型z得分的情况），MAD约为0.68，这解释了缩放因子。因此，选择m = 3.5意味着您想要丢弃0.05％的数据。 - Fato39

Answer 2

这个方法与您的方法几乎相同，只是更针对Numpy（也仅适用于Numpy数组）:

def reject_outliers(data, m=2):
    return data[abs(data - np.mean(data)) < m * np.std(data)]

- eumiro

---

如果“m”足够大（例如“m=6”），那么这种方法就足够好，但对于小的“m”值来说，这种方法会受到均值和方差不是强鲁棒估计量的影响。 - Benjamin Bannier

这并不是对该方法的抱怨，而是对“异常值”这个模糊概念的抱怨。 - Eelco Hoogendoorn

你如何选择一个m？ - john k

我还没有让它正常工作。我一直收到一个错误：return data[abs(data - np.mean(data)) < m * np.std(data)] TypeError: 只有整数标量数组可以转换为标量索引，或者它会导致我的程序冻结。 - john k

@johnktejik 数据参数需要是一个numpy数组。 - Sander van Leeuwen

Answer 3

Benjamin Bannier的答案在中位数到距离中位数的距离为0时产生了一个传递效应，因此我发现下面这个修改后的版本对于下面示例中给出的情况更有帮助。

def reject_outliers_2(data, m=2.):
    d = np.abs(data - np.median(data))
    mdev = np.median(d)
    s = d / (mdev if mdev else 1.)
    return data[s < m]

例子：

data_points = np.array([10, 10, 10, 17, 10, 10])
print(reject_outliers(data_points))
print(reject_outliers_2(data_points))

给出:

[[10, 10, 10, 17, 10, 10]]  # 17 is not filtered
[10, 10, 10, 10, 10]  # 17 is filtered (it's distance, 7, is greater than m)

- Egal

---

你是什么意思，Oleg？难道简单的一维数组不起作用吗？ - Egal

你是什么意思，Oleg？难道简单的一维数组不起作用吗？ - undefined

好的，我简单明了地表达一下。由于保密协议，我不能分享我的代码，但是这个解决方案并没有给我预期的答案。我已经删除了我的回答，因为它可能会误导别人，在其他情况下，它可能有效。用户ankostis提供的下面的解决方案起了作用。 - Олег Місько

嗯，我简单明了地表达了。由于保密协议，我不能分享我的代码，但是这个解决方案并没有给我预期的答案。我删除了我的回答，因为它可能会误导他人，在他们的情况下可能有效。用户ankostis下面的解决方案解决了问题。 - undefined

Answer 4

基于Benjamin的工作，使用 pandas.Series，并用 IQR替换MAD：

def reject_outliers(sr, iq_range=0.5):
    pcnt = (1 - iq_range) / 2
    qlow, median, qhigh = sr.dropna().quantile([pcnt, 0.50, 1-pcnt])
    iqr = qhigh - qlow
    return sr[ (sr - median).abs() <= iqr]

例如，如果您设置iq_range=0.6，四分位距的百分位数将变为：0.20 <--> 0.80，因此将包括更多的异常值。

- ankostis

Answer 5

一种替代方法是对标准差进行稳健估计（假设符合高斯统计）。在查找在线计算器时，我发现90％的百分位数对应于1.2815σ，而95％的百分位数则是1.645σ (http://vassarstats.net/tabs.html?#z)。

举个简单例子：

import numpy as np

# Create some random numbers
x = np.random.normal(5, 2, 1000)

# Calculate the statistics
print("Mean= ", np.mean(x))
print("Median= ", np.median(x))
print("Max/Min=", x.max(), " ", x.min())
print("StdDev=", np.std(x))
print("90th Percentile", np.percentile(x, 90))

# Add a few large points
x[10] += 1000
x[20] += 2000
x[30] += 1500

# Recalculate the statistics
print()
print("Mean= ", np.mean(x))
print("Median= ", np.median(x))
print("Max/Min=", x.max(), " ", x.min())
print("StdDev=", np.std(x))
print("90th Percentile", np.percentile(x, 90))

# Measure the percentile intervals and then estimate Standard Deviation of the distribution, both from median to the 90th percentile and from the 10th to 90th percentile
p90 = np.percentile(x, 90)
p10 = np.percentile(x, 10)
p50 = np.median(x)
# p50 to p90 is 1.2815 sigma
rSig = (p90-p50)/1.2815
print("Robust Sigma=", rSig)

rSig = (p90-p10)/(2*1.2815)
print("Robust Sigma=", rSig)

我收到的输出是：

Mean=  4.99760520022
Median=  4.95395274981
Max/Min= 11.1226494654   -2.15388472011
Sigma= 1.976629928
90th Percentile 7.52065379649

Mean=  9.64760520022
Median=  4.95667658782
Max/Min= 2205.43861943   -2.15388472011
Sigma= 88.6263902244
90th Percentile 7.60646688694

Robust Sigma= 2.06772555531
Robust Sigma= 1.99878292462

这个结果接近于期望值2。

如果我们想要删除超过5个标准差的点（在1000个数据点中，预计有1个值> 3个标准差）：

y = x[abs(x - p50) < rSig*5]

# Print the statistics again
print("Mean= ", np.mean(y))
print("Median= ", np.median(y))
print("Max/Min=", y.max(), " ", y.min())
print("StdDev=", np.std(y))

这将会产生：

Mean=  4.99755359935
Median=  4.95213030447
Max/Min= 11.1226494654   -2.15388472011
StdDev= 1.97692712883

我不知道哪种方法更有效/更健壮

- Chris

Answer 6

我希望在这个答案中提供两种方法，一种基于 "z-score" ，另一种基于 "IQR"。

此答案提供的代码适用于单维度 numpy 数组和多维度 numpy 数组。

首先让我们导入一些模块。

import collections
import numpy as np
import scipy.stats as stat
from scipy.stats import iqr

基于z得分的方法

这种方法将检验数字是否超出了三个标准差。根据这个规则，如果值是异常值，则该方法返回true，否则返回false。

def sd_outlier(x, axis = None, bar = 3, side = 'both'):
    assert side in ['gt', 'lt', 'both'], 'Side should be `gt`, `lt` or `both`.'

    d_z = stat.zscore(x, axis = axis)

    if side == 'gt':
        return d_z > bar
    elif side == 'lt':
        return d_z < -bar
    elif side == 'both':
        return np.abs(d_z) > bar

基于四分位距的方法

该方法将检测值是否小于q1 - 1.5 * iqr或大于q3 + 1.5 * iqr，类似于SPSS的绘图方法。

def q1(x, axis = None):
    return np.percentile(x, 25, axis = axis)

def q3(x, axis = None):
    return np.percentile(x, 75, axis = axis)

def iqr_outlier(x, axis = None, bar = 1.5, side = 'both'):
    assert side in ['gt', 'lt', 'both'], 'Side should be `gt`, `lt` or `both`.'

    d_iqr = iqr(x, axis = axis)
    d_q1 = q1(x, axis = axis)
    d_q3 = q3(x, axis = axis)
    iqr_distance = np.multiply(d_iqr, bar)

    stat_shape = list(x.shape)

    if isinstance(axis, collections.Iterable):
        for single_axis in axis:
            stat_shape[single_axis] = 1
    else:
        stat_shape[axis] = 1

    if side in ['gt', 'both']:
        upper_range = d_q3 + iqr_distance
        upper_outlier = np.greater(x - upper_range.reshape(stat_shape), 0)
    if side in ['lt', 'both']:
        lower_range = d_q1 - iqr_distance
        lower_outlier = np.less(x - lower_range.reshape(stat_shape), 0)

    if side == 'gt':
        return upper_outlier
    if side == 'lt':
        return lower_outlier
    if side == 'both':
        return np.logical_or(upper_outlier, lower_outlier)

最后，如果你想过滤掉异常值，请使用 numpy 选择器。

祝您拥有美好的一天。

- Losses Don

Answer 7

我想做类似的事情，只是将数字设为NaN而不是从数据中删除它，因为如果你删除它，你会改变长度，这可能会破坏绘图（例如，如果你只从表格中的一列中删除异常值，但你需要让它保持与其他列相同，以便可以相互绘制）。 为此，我使用了numpy的掩码函数：

def reject_outliers(data, m=2):
    stdev = np.std(data)
    mean = np.mean(data)
    maskMin = mean - stdev * m
    maskMax = mean + stdev * m
    mask = np.ma.masked_outside(data, maskMin, maskMax)
    print('Masking values outside of {} and {}'.format(maskMin, maskMax))
    return mask

- Alex S

---

你也可以使用 np.clip 将它们剪裁到允许的最小和最大值以保持维度。 - Andi R

Answer 8

考虑到标准差由于极端值的存在变得非常大时，上述所有方法均会失败。

（与平均值计算类似，应该计算中位数。然而，平均值“更容易出现这种标准差误差”。）

您可以尝试迭代地应用算法或使用四分位距进行过滤： （此处的“因子”与n* sigma范围有关，但仅当数据遵循高斯分布时才如此）

import numpy as np

def sortoutOutliers(dataIn,factor):
    quant3, quant1 = np.percentile(dataIn, [75 ,25])
    iqr = quant3 - quant1
    iqrSigma = iqr/1.34896
    medData = np.median(dataIn)
    dataOut = [ x for x in dataIn if ( (x > medData - factor* iqrSigma) and (x < medData + factor* iqrSigma) ) ] 
    return(dataOut)

- K. Foe

---

抱歉，我忽略了上面已经有一个IQR建议。由于代码更短，我应该保留这个答案还是删除它？ - K. Foe

Answer 9

有很多答案，但我添加了一个新的可能对作者或其他用户有用的答案。

您可以使用Hampel滤波器。但您需要使用Series进行操作。

Hampel滤波器返回异常值索引，然后您可以从Series中删除它们，然后将其转换回List。

要使用Hampel滤波器，您可以轻松地使用pip安装软件包：

pip install hampel

使用方法：

# Imports
from hampel import hampel
import pandas as pd

list_d = [2, 4, 5, 1, 6, 5, 40]

# List to Series
time_series = pd.Series(list_d)

# Outlier detection with Hampel filter
# Returns the Outlier indices
outlier_indices = hampel(ts = time_series, window_size = 3)

# Drop Outliers indices from Series
filtered_d = time_series.drop(outlier_indices)

filtered_d.values.tolist()

print(f'filtered_d: {filtered_d.values.tolist()}')

输出结果如下：

filtered_d: [2, 4, 5, 1, 6, 5]

其中，ts 是一个 pandas 的 Series 对象，window_size 是一个总窗口大小，将被计算为 2 * window_size + 1。

对于这个序列，我将 window_size 设置为 3。

与 Series 工作的很酷的一件事是能够生成图形：

# Imports
import matplotlib.pyplot as plt

plt.style.use('seaborn-darkgrid')

# Plot Original Series
time_series.plot(style = 'k-')
plt.title('Original Series')
plt.show()
    
# Plot Cleaned Series
filtered_d.plot(style = 'k-')
plt.title('Cleaned Series (Without detected Outliers)')
plt.show()

输出结果将是：

如果想了解更多有关Hampel滤波器的内容，我推荐以下阅读材料：

• Python实现的Hampel滤波器
• Hampel滤波器进行异常值检测
• 使用Hampel滤波器清洗时间序列数据

- whoisraibolt

Answer 10

在这里，我找到了x中的异常值，并用它们附近的一组点（win）的中位数替换它们（参考Benjamin Bannier的中位偏差）

def outlier_smoother(x, m=3, win=3, plots=False):
    ''' finds outliers in x, points > m*mdev(x) [mdev:median deviation] 
    and replaces them with the median of win points around them '''
    x_corr = np.copy(x)
    d = np.abs(x - np.median(x))
    mdev = np.median(d)
    idxs_outliers = np.nonzero(d > m*mdev)[0]
    for i in idxs_outliers:
        if i-win < 0:
            x_corr[i] = np.median(np.append(x[0:i], x[i+1:i+win+1]))
        elif i+win+1 > len(x):
            x_corr[i] = np.median(np.append(x[i-win:i], x[i+1:len(x)]))
        else:
            x_corr[i] = np.median(np.append(x[i-win:i], x[i+1:i+win+1]))
    if plots:
        plt.figure('outlier_smoother', clear=True)
        plt.plot(x, label='orig.', lw=5)
        plt.plot(idxs_outliers, x[idxs_outliers], 'ro', label='outliers')                                                                                                                    
        plt.plot(x_corr, '-o', label='corrected')
        plt.legend()
    
    return x_corr

- scrx2