Pythonic 检测一维观测数据中异常值的方法

Posted 2023-03-12

【中文标题】Pythonic 检测一维观测数据中异常值的方法

【英文标题】：Pythonic way of detecting outliers in one dimensional observation data

【发布时间】：2014-04-16 17:44:47

【问题描述】：

对于给定的数据，我想将异常值（由 95% 置信水平或 95% 分位数函数或任何所需的值定义）设置为 nan 值。以下是我现在正在使用的数据和代码。如果有人能进一步解释我，我会很高兴。

import numpy as np, matplotlib.pyplot as plt

data = np.random.rand(1000)+5.0

plt.plot(data)
plt.xlabel('observation number')
plt.ylabel('recorded value')
plt.show()

【问题讨论】：

您更了解您的数据，但我认为 Winsorising 会比删除更好。此外，如果将这些数据设置为 nan，则必须处理它。看看np.percentile函数。

仅供参考：Detect and exclude outliers in Pandas dataframe

【参考方案1】：

按照@Martin 的建议使用np.percentile：

percentiles = np.percentile(data, [2.5, 97.5])

# or =>, <= for within 95%
data[(percentiles[0]<data) & (percentiles[1]>data)]

# set the outliners to np.nan
data[(percentiles[0]>data) | (percentiles[1]<data)] = np.nan

【讨论】：

使用数据的百分位数作为异常值测试是合理的第一步，但并不理想。问题是 1）您将删除一些数据，即使它不是异常值，以及 2）异常值严重影响方差，因此影响百分位值。最常见的异常值测试使用“中值绝对偏差”，它对异常值的存在不太敏感。

@Joe Kington 如果您能使用 python 代码实现您的方式，我将不胜感激。

@Joe Kington 我看到了你回答的链接。但是，没有更简单的方法可以使其主要使用 numpy 中可用的功能

@julie - 该函数广泛使用 numpy（它需要一个 numpy 数组作为输入并输出一个 numpy 数组）。异常值测试远远超出numpy 的范围。 （numpy 本身只包含核心数据结构和一些基本操作。它故意很小。）你可以说scipy.stats 将是一个合理的异常值测试位置，但其中有很多，而且有没有单一的最佳测试。因此，目前没有单功能异常值测试。

Statsmodels 在sm.robust.mad 中有一个中值绝对偏差函数。我不确定是否有用于单变量异常值测试的工具，但在回归框架中有影响/异常值的工具。将看到添加一些用于单变量异常值检测的工具。

【参考方案2】：

一个简单的解决方案也可以是，删除超出 2 个标准差（或 1.96）的东西：

import random
def outliers(tmp):
    """tmp is a list of numbers"""
    outs = []
    mean = sum(tmp)/(1.0*len(tmp))
    var = sum((tmp[i] - mean)**2 for i in range(0, len(tmp)))/(1.0*len(tmp))
    std = var**0.5
    outs = [tmp[i] for i in range(0, len(tmp)) if abs(tmp[i]-mean) > 1.96*std]
    return outs


lst = [random.randrange(-10, 55) for _ in range(40)]
print lst
print outliers(lst)

【讨论】：

这是用于 python 2 的吗？

在 python 3 中我应该用什么代替xrange？

python 2 中的 xrange 与 python 3 中的 range 相同。python 3 中不再有 xrange。

【参考方案3】：

一维数据中异常值的检测取决于其分布

1- 正态分布：

数据值几乎均匀分布在预期范围内： 在这种情况下，您可以轻松使用包括均值在内的所有方法，例如对于正态分布数据（中心极限定理和样本均值的抽样分布）的 3 或 2 个标准差（95% 或 99.7%）的置信区间。I 是一种非常有效的方法。 可汗学院统计和概率 - 抽样分布库中对此进行了解释。

如果您想要数据点的置信区间而不是平均值，另一种方法是预测区间。

数据值随机分布在一个范围内： 平均值可能不是数据的公平表示，因为平均值很容易受到异常值的影响（数据集中非常小的或非常大的值，不典型） 中位数是衡量数值数据集中心的另一种方法。

中值绝对偏差 - 一种以中值距离衡量所有点与中值距离的方法 http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h.htm - 有一个很好的解释，正如上面乔金顿的回答所解释的那样

2 - 对称分布：如果 z-score 计算和阈值相应地改变，中值绝对偏差也是一个很好的方法

说明： http://eurekastatistics.com/using-the-median-absolute-deviation-to-find-outliers/

3 - 非对称分布：双 MAD - 双中值绝对偏差 上面附加链接中的说明

附上我的python代码供参考：

 def is_outlier_doubleMAD(self,points):
    """
    FOR ASSYMMETRIC DISTRIBUTION
    Returns : filtered array excluding the outliers

    Parameters : the actual data Points array

    Calculates median to divide data into 2 halves.(skew conditions handled)
    Then those two halves are treated as separate data with calculation same as for symmetric distribution.(first answer) 
    Only difference being , the thresholds are now the median distance of the right and left median with the actual data median
    """

    if len(points.shape) == 1:
        points = points[:,None]
    median = np.median(points, axis=0)
    medianIndex = (points.size/2)

    leftData = np.copy(points[0:medianIndex])
    rightData = np.copy(points[medianIndex:points.size])

    median1 = np.median(leftData, axis=0)
    diff1 = np.sum((leftData - median1)**2, axis=-1)
    diff1 = np.sqrt(diff1)

    median2 = np.median(rightData, axis=0)
    diff2 = np.sum((rightData - median2)**2, axis=-1)
    diff2 = np.sqrt(diff2)

    med_abs_deviation1 = max(np.median(diff1),0.000001)
    med_abs_deviation2 = max(np.median(diff2),0.000001)

    threshold1 = ((median-median1)/med_abs_deviation1)*3
    threshold2 = ((median2-median)/med_abs_deviation2)*3

    #if any threshold is 0 -> no outliers
    if threshold1==0:
        threshold1 = sys.maxint
    if threshold2==0:
        threshold2 = sys.maxint
    #multiplied by a factor so that only the outermost points are removed
    modified_z_score1 = 0.6745 * diff1 / med_abs_deviation1
    modified_z_score2 = 0.6745 * diff2 / med_abs_deviation2

    filtered1 = []
    i = 0
    for data in modified_z_score1:
        if data < threshold1:
            filtered1.append(leftData[i])
        i += 1
    i = 0
    filtered2 = []
    for data in modified_z_score2:
        if data < threshold2:
            filtered2.append(rightData[i])
        i += 1

    filtered = filtered1 + filtered2
    return filtered

【讨论】：

在 Python 3 中，它应该是 medianIndex = int(points.size/2)。此外，如果我运行代码并将阈值设置为零，它会崩溃并显示消息name 'sys' is not defined。 Laslty，函数调用中的self 永远不会被使用。

也可以使用：medianIndex = points.size//2 来避免浮动值

【参考方案4】：

我已经修改了http://eurekastatistics.com/using-the-median-absolute-deviation-to-find-outliers 的代码，它给出了与 Joe Kington 相同的结果，但使用 L1 距离而不是 L2 距离，并且支持非对称分布。原始的 R 代码没有 Joe 的 0.6745 乘数，所以我还添加了它以在这个线程中保持一致性。不是 100% 确定是否有必要，但可以进行比较。

def doubleMADsfromMedian(y,thresh=3.5):
    # warning: this function does not check for NAs
    # nor does it address issues when 
    # more than 50% of your data have identical values
    m = np.median(y)
    abs_dev = np.abs(y - m)
    left_mad = np.median(abs_dev[y <= m])
    right_mad = np.median(abs_dev[y >= m])
    y_mad = left_mad * np.ones(len(y))
    y_mad[y > m] = right_mad
    modified_z_score = 0.6745 * abs_dev / y_mad
    modified_z_score[y == m] = 0
    return modified_z_score > thresh

【讨论】：

如何在多元数据上使用基于 MAD 的方法？您提到的文章很棒，但我猜它适用于单维数据。我想知道修改它以使其也适用于多变量数据的最简单方法。

对于多变量数据，没有简单的方法可以做到这一点。一种简单的方法是一次只将该方法应用于一个变量，并查看某些样本是否在任何维度上都是异常值。

@sergeyf 我们如何选择阈值？通读原帖，那里也挖不出来。

@ekta 通常和不幸的答案是用你的实际数据尝试一堆不同的阈值，看看有什么遗漏。由于这是一维的，因此您可以像 Joe Kington 的回答那样进行可视化。如果它更容易，您可以将阈值视为类似于标准偏差的数量。所以3.5很多。但我之前使用过接近 6 的数字 - 仅取决于您的数据。

@TheAG 啊，我明白你在说什么。我之所以说绝对是因为我们不在乎异常值是在右尾还是左尾。但是，如果您在乎，那么删除绝对值是有道理的！

【参考方案5】：

使用percentile 的问题在于，被识别为异常值的点是您的样本量的函数。

测试异常值的方法有很多种，您应该考虑如何对它们进行分类。理想情况下，您应该使用先验信息（例如，“任何高于/低于此值的东西都是不现实的，因为......”）

然而，一个常见的、不太合理的异常值测试是根据“中值绝对偏差”删除点。

这是 N 维情况的实现（来自本文的一些代码：https://github.com/joferkington/oost_paper_code/blob/master/utilities.py）：

def is_outlier(points, thresh=3.5):
    """
    Returns a boolean array with True if points are outliers and False 
    otherwise.

    Parameters:
    -----------
        points : An numobservations by numdimensions array of observations
        thresh : The modified z-score to use as a threshold. Observations with
            a modified z-score (based on the median absolute deviation) greater
            than this value will be classified as outliers.

    Returns:
    --------
        mask : A numobservations-length boolean array.

    References:
    ----------
        Boris Iglewicz and David Hoaglin (1993), "Volume 16: How to Detect and
        Handle Outliers", The ASQC Basic References in Quality Control:
        Statistical Techniques, Edward F. Mykytka, Ph.D., Editor. 
    """
    if len(points.shape) == 1:
        points = points[:,None]
    median = np.median(points, axis=0)
    diff = np.sum((points - median)**2, axis=-1)
    diff = np.sqrt(diff)
    med_abs_deviation = np.median(diff)

    modified_z_score = 0.6745 * diff / med_abs_deviation

    return modified_z_score > thresh

这和one of my previous answers很像，不过我想详细说明一下样本量的影响。

让我们将基于百分位数的异常值检验（类似于@CTZhu 的答案）与针对各种不同样本大小的中值绝对偏差 (MAD) 检验进行比较：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def main():
    for num in [10, 50, 100, 1000]:
        # Generate some data
        x = np.random.normal(0, 0.5, num-3)

        # Add three outliers...
        x = np.r_[x, -3, -10, 12]
        plot(x)

    plt.show()

def mad_based_outlier(points, thresh=3.5):
    if len(points.shape) == 1:
        points = points[:,None]
    median = np.median(points, axis=0)
    diff = np.sum((points - median)**2, axis=-1)
    diff = np.sqrt(diff)
    med_abs_deviation = np.median(diff)

    modified_z_score = 0.6745 * diff / med_abs_deviation

    return modified_z_score > thresh

def percentile_based_outlier(data, threshold=95):
    diff = (100 - threshold) / 2.0
    minval, maxval = np.percentile(data, [diff, 100 - diff])
    return (data < minval) | (data > maxval)

def plot(x):
    fig, axes = plt.subplots(nrows=2)
    for ax, func in zip(axes, [percentile_based_outlier, mad_based_outlier]):
        sns.distplot(x, ax=ax, rug=True, hist=False)
        outliers = x[func(x)]
        ax.plot(outliers, np.zeros_like(outliers), 'ro', clip_on=False)

    kwargs = dict(y=0.95, x=0.05, ha='left', va='top')
    axes[0].set_title('Percentile-based Outliers', **kwargs)
    axes[1].set_title('MAD-based Outliers', **kwargs)
    fig.suptitle('Comparing Outlier Tests with n='.format(len(x)), size=14)

main()

---

---

---

---

请注意，无论样本量如何，基于 MAD 的分类器都能正常工作，而基于百分位数的分类器分类的点越多，样本量越大，无论它们是否实际上是异常值。

【讨论】：

乔，+1，这是一个很好的答案。虽然我想知道，如果 OP 的数据总是一致地受到干扰（random.rand()），或者大多数时候可能总是遵循其他分布。如果数据总是受到一致的干扰，我不确定是否使用MAD。

@CTZhu - 很好，特别是如果 OP 的数据是对数正态分布的。对于模糊对称分布，与正态分布的偏差应该没有太大关系，但对于强不对称分布（例如对数正态），MAD 不是一个好的选择。 （尽管你总是可以将它应用到日志空间来解决这个问题。）所有这些只是为了强调一点，你应该对你选择的任何异常值测试进行一些思考。

@JoeKington 您使用中位数的每个地方，但 diff 计算为 L2 范数（**2）；中位数是最小化 L1 范数的值，而在 L2 范数中，“均值”是中心；我期待着如果你从 L1 标准的中位数开始。你有什么理由 **2 在计算 diff 时比绝对值更好吗？

在您基于 mad 的异常值检测中，您如何设置值 0.6745 和 3.5 用于什么目的，如何确定这些值？这很混乱

@JoeKington 的 PDF 论文 pdf-archive.com/2016/07/29/outlier-methods-external/… 的备用镜像