是否有一个 numpy 内置函数可以拒绝列表中的异常值

问题描述：

是否有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我在这里使用的一样）。

解决方案 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 数组。

解决方案 2：

此方法与您的方法几乎相同，只是更多的 numpyst（也仅适用于 numpy 数组）：

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

解决方案 3：

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

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)

解决方案 4：

在本杰明 (Benjamin) 的基础上，使用，并用 IQRpandas.Series替换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，因此将包含更多异常值。

解决方案 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

我不知道哪种方法更高效/更稳健

解决方案 6：

我想做一些类似的事情，除了将数字设置为 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

解决方案 7：

我想在这个答案中提供两种方法，基于“z 分数”的解决方案和基于“IQR”的解决方案。

此答案中提供的代码适用于单个 dimnumpy数组和多个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

基于IQR的方法

此方法将测试值是否小于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选择器。

祝你今天过得愉快。

解决方案 8：

请考虑一下，当您的标准差由于巨大的异常值而变得非常大时，上述所有方法都会失败。

（类似于平均值计算失败，而应该计算中位数。尽管平均值“更容易出现像 stdDv 这样的错误”。）

您可以尝试迭代应用您的算法或使用四分位距进行过滤：（此处“因子”与*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)

解决方案 9：

在这里，我找到了中的异常值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

解决方案 10：

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

您可以使用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()}')

输出将是：

已筛选：[2, 4, 5, 1, 6, 5]

其中，ts是一个 pandasSeries对象，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 滤波器的更多信息，我推荐以下读物：

• Hampel 滤波器的 Python 实现

• 使用 Hampel 滤波器进行异常值检测

• 使用 Hampel 过滤器清理时间序列数据

解决方案 11：

我的解决方案是让异常值等于之前的值。

test_data = [2,4,5,1,6,5,40, 3]
def reject_outliers(data, m=2):
    mean = np.mean(data)
    std = np.std(data)
    for i in range(len(data)) :
        if np.abs(data[i] -mean) > m*std :
            data[i] = data[i-1]
    return data
reject_outliers(test_data)

输出：

[2, 4, 5, 1, 6, 5, 5, 3]

解决方案 12：

如果您想获得异常值的索引位置，idx_list它将返回它。

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

data_points = np.array([8, 10, 35, 17, 73, 77])  
print(reject_outliers(data_points))

after rejection: [ 8 10 35 17], index positions of outliers: [4 5]

解决方案 13：

对于一组图像（每个图像有 3 个维度），我想拒绝我使用的每个像素的异常值：

mean = np.mean(imgs, axis=0)
std = np.std(imgs, axis=0)
mask = np.greater(0.5 * std + 1, np.abs(imgs - mean))
masked = np.multiply(imgs, mask)

然后可以计算平均值：

masked_mean = np.divide(np.sum(masked, axis=0), np.sum(mask, axis=0))

（我用它来进行背景减法）

解决方案 14：

沿轴修剪 NumPy 数组中的异常值，并将其替换为该轴上的最小值或最大值（取两者中较接近的值）。阈值为 z 值：

def np_z_trim(x, threshold=10, axis=0):
    """ Replace outliers in numpy ndarray along axis with min or max values
    within the threshold along this axis, whichever is closer."""
    mean = np.mean(x, axis=axis, keepdims=True)
    std = np.std(x, axis=axis, keepdims=True)
    masked = np.where(np.abs(x - mean) < threshold * std, x, np.nan)
    min = np.nanmin(masked, axis=axis, keepdims=True)
    max = np.nanmax(masked, axis=axis, keepdims=True)
    repl = np.where(np.abs(x - max) < np.abs(x - min), max, min)
    return np.where(np.isnan(masked), repl, masked)

解决方案 15：

我的解决方案是删除顶部和底部的百分位数，保留等于边界的值：

def remove_percentile_outliers(data, percent_to_drop=0.001):
    low, high = data.quantile([percent_to_drop / 2, 1-percent_to_drop / 2])
    return data[(data >= low )&(data <= high)]