计算百分位数以去除异常值的快速算法

Posted 2023-03-12

【中文标题】计算百分位数以去除异常值的快速算法

【英文标题】：Fast Algorithm for computing percentiles to remove outliers

【发布时间】：2011-04-16 08:06:39

【问题描述】：

我有一个程序需要重复计算数据集的近似百分位数（顺序统计），以便在进一步处理之前删除异常值。我目前正在通过对值数组进行排序并选择适当的元素来做到这一点；这是可行的，但它在配置文件中是一个明显的亮点，尽管它只是该程序的一个相当小的部分。

更多信息：

该数据集包含多达 100000 个浮点数，并假定“合理”分布 - 特定值附近的密度不太可能重复或出现巨大峰值；如果由于某种奇怪的原因分布是奇怪的，那么近似值不太准确是可以的，因为数据可能无论如何都搞砸了，进一步的处理也很可疑。但是，数据不一定是均匀分布或正态分布的；它不太可能退化。 一个近似的解决方案会很好，但我确实需要了解该近似如何引入错误以确保它是有效的。 由于目标是消除异常值，我一直在计算相同数据的两个百分位数：例如一个为 95%，一个为 5%。 该应用程序使用 C# 语言，在 C++ 中进行了一些繁重的工作；伪代码或预先存在的库都可以。 一种完全不同的去除异常值的方法也可以，只要它是合理的。

更新：

看来我正在寻找一个近似的selection algorithm。

虽然这一切都是在一个循环中完成的，但数据每次都（略有）不同，因此像 for this question 那样重用数据结构并不容易。

实施的解决方案

使用 Gronim 建议的***选择算法将这部分运行时间减少了大约 20 倍。

由于我找不到 C# 实现，这就是我想出的。即使对于小输入，它也比 Array.Sort 更快；在 1000 个元素时，它的速度提高了 25 倍。

public static double QuickSelect(double[] list, int k) 
    return QuickSelect(list, k, 0, list.Length);

public static double QuickSelect(double[] list, int k, int startI, int endI) 
    while (true) 
        // Assume startI <= k < endI
        int pivotI = (startI + endI) / 2; //arbitrary, but good if sorted
        int splitI = partition(list, startI, endI, pivotI);
        if (k < splitI)
            endI = splitI;
        else if (k > splitI)
            startI = splitI + 1;
        else //if (k == splitI)
            return list[k];
    
    //when this returns, all elements of list[i] <= list[k] iif i <= k

static int partition(double[] list, int startI, int endI, int pivotI) 
    double pivotValue = list[pivotI];
    list[pivotI] = list[startI];
    list[startI] = pivotValue;

    int storeI = startI + 1;//no need to store @ pivot item, it's good already.
    //Invariant: startI < storeI <= endI
    while (storeI < endI && list[storeI] <= pivotValue) ++storeI; //fast if sorted
    //now storeI == endI || list[storeI] > pivotValue
    //so elem @storeI is either irrelevant or too large.
    for (int i = storeI + 1; i < endI; ++i)
        if (list[i] <= pivotValue) 
            list.swap_elems(i, storeI);
            ++storeI;
        
    int newPivotI = storeI - 1;
    list[startI] = list[newPivotI];
    list[newPivotI] = pivotValue;
    //now [startI, newPivotI] are <= to pivotValue && list[newPivotI] == pivotValue.
    return newPivotI;

static void swap_elems(this double[] list, int i, int j) 
    double tmp = list[i];
    list[i] = list[j];
    list[j] = tmp;

感谢 Gronim，为我指明了正确的方向！

【问题讨论】：

我知道这是旧的，但我已经实现了这个，但是我如何实际执行它以获得第 5 个和第 95 个百分位数？

如果你想要第 5 个百分位，QuickSelect(values, (int)(values.Length*0.05+0.5))。如果您想要第 95 个百分位数，QuickSelect(values, (int)(values.Length*0.95+0.5)) - 请注意，您必须将 0.05 和 0.95 的小数索引四舍五入到整个索引，至少除非您的列表长度是 20 的倍数。如果您的列表很短，您可以考虑插值而不是仅仅选择一个索引，但对于大多数用法我怀疑它是否重要 - 如果你不关心只选择 5%，你可能不在乎它实际上是 4.8% 还是其他 - 无论如何都没有确切的第 5 个百分位数，一般来说。

【参考方案1】：

将数据的最小值和最大值之间的间隔划分为（例如）1000 个 bin 并计算直方图。然后构建部分总和，看看它们首先超过 5000 或 95000 的位置。

【讨论】：

不错...快速排序，并切断了顶部和底部的5000。不知道分布不知道如何做得更好。

桶排序更适合这个。

这听起来非常实用，尽管并不总是有效。一些极端的异常值真的会扭曲你的垃圾箱......

【参考方案2】：

不是专家，但我的记忆表明：

要准确确定您需要排序和计数的百分位数 从数据中抽取样本并计算百分位值听起来像是一个不错的近似方案，如果你能得到一个好的样本的话 如果不是，按照 Henrik 的建议，如果您对桶进行计数并计算它们，则可以避免完全排序

【参考方案3】：

您可以仅从数据集的一部分（例如前几千个点）估计百分位数。

Glivenko–Cantelli theorem 确保这是一个相当不错的估计，如果您可以假设您的数据点是独立的。

【讨论】：

不幸的是，数据点不是独立的，它们是按外部标准排序的——但我可以按随机顺序迭代。我不明白链接定理实际上如何让我估计百分位数 - 你能举个例子吗？为正态分布？

@Eamon：链接定理简单地说，经验分布函数（根据数据计算百分位数时会隐含使用）是对实际分布的良好估计。你不必实际使用它 =)

啊，好吧，我明白你的意思了:-)

【参考方案4】：

我能想到几个基本的方法。首先是计算范围（通过找到最高和最低值），将每个元素投影到一个百分位数（（x - min）/范围），并丢弃任何低于 0.05 或高于 0.95 的元素。

第二个是计算均值和标准差。距平均值 2 个标准差的跨度（在两个方向上）将包含 95% 的正态分布样本空间，这意味着您的异常值将在 97.5 个百分位数内。计算一个系列的平均值是线性的，标准 dev 也是线性的（每个元素的差值与平均值之和的平方根）。然后，从平均值中减去 2 个 sigma，然后在平均值上加上 2 个 sigma，就得到了异常值限制。

这两个都将在大致线性的时间内计算；第一个需要两次通过，第二个需要三个（一旦你有你的限制，你仍然必须丢弃异常值）。由于这是一个基于列表的操作，我认为你不会找到任何具有对数或恒定复杂度的东西。任何进一步的性能提升都需要优化迭代和计算，或者通过对子样本（例如每三个元素）执行计算来引入错误。

【讨论】：

第一个建议不是扔掉外面的第 5 个百分位数，而是根据最极端的异常值做一些事情，这是非常不稳定的。第二个建议假设数据是正态分布的，但显然不是。

【参考方案5】：

我曾经通过计算standard deviation 来识别异常值。距离大于平均值标准偏差的 2（或 3）倍的所有事物都是异常值。 2 次 = 约 95%。

由于您正在计算平均值，因此计算标准差也非常容易，速度非常快。

您也可以只使用数据的一个子集来计算数字。

【讨论】：

数据不是正态分布的。

【参考方案6】：

Henrik 的直方图解决方案将起作用。您还可以使用选择算法在 O(n) 中有效地找到包含 n 个元素的数组中的 k 个最大或最小元素。要将其用于第 95 个百分位数，请设置 k=0.05n 并找到 k 个最大的元素。

参考：

http://en.wikipedia.org/wiki/Selection_algorithm#Selecting_k_smallest_or_largest_elements

【讨论】：

对，这就是我要找的——选择算法！

【参考方案7】：

According 对其创建者SoftHeap 可用于：

计算精确或近似中位数 和最佳百分比。也是 用于近似排序...

【讨论】：

@Eamon SoftHeap 背后的整个想法及其应用程序真的很酷。

@EugenConstantinDinca：谢谢你的好主意！是否在某处有实际的实现，或者论文/维基是唯一的来源？

@Legend 我已经找到了它的一些不同语言（从 C++ 到 Haskell）的实现，但没有使用过，所以我不知道它们有多大用处。

@EugenConstantinDinca：哦，我明白了。感谢您的信息。

【参考方案8】：

您的问题的一个很好的一般答案似乎是RANSAC。 给定一个模型和一些噪声数据，该算法有效地恢复了模型的参数。 您将不得不选择一个可以映射您的数据的简单模型。任何光滑的东西都应该没问题。假设是几个高斯的混合体。 RANSAC 将设置模型的参数并同时估计一组内联。然后扔掉任何不适合模型的东西。

【讨论】：

我有一组数字 - 不是一些复杂的模型 - RANSAC 看起来会很慢而且容易出错，而且对于这样一个简单的情况，存在更好的解决方案。

【参考方案9】：

一组 100k 个元素的数据几乎不需要时间来排序，所以我假设你必须重复这样做。如果数据集是相同的数据集，只是略有更新，最好构建一棵树 (O(N log N))，然后在新点进入时删除和添加新点（O(K log N) 其中K 是更改的点数）。否则，已经提到的kth 最大元素解决方案会为每个数据集提供O(N)。

【参考方案10】：

即使数据不是正态分布的，您也可以过滤掉 2 或 3 个标准差；至少，它将以一致的方式完成，这应该很重要。

当您删除异常值时，std dev 会发生变化，您可以循环执行此操作，直到 std dev 的变化最小。您是否要这样做取决于您为什么要以这种方式处理数据。一些统计学家对去除异常值持重大保留意见。但是有些人会删除异常值以证明数据是相当正态分布的。

【讨论】：

如果数据大部分位于极端情况下——即与正常情况相反，如果你愿意的话——那么这种方法可能会删除大量数据。我真的不想删除超过一小部分的数据，最好只在这些是异常值时删除。我正在抑制异常值，因为它们会分散注意力 - 它们只是从可视化中裁剪出来的，而不是从实际数据中裁剪出来的。

根据定义，只有一小部分数据可能处于极端状态。根据切比雪夫不等式，只有 1/9 的分布可以超过 3 个标准差；只有 1/16 可以相差 4 个偏差。而这些限制只有在你的分布只有两个尖峰的退化情况下才会达到。因此，计算 O(N) 中的偏差是过滤异常值的一种有效且有效的方法。

@MSalters：（是的，回复了 3 年的评论）：切比雪夫不等式不够精确，无法实用。要裁剪到至少 95% 的数据集，我需要执行 4.5 sigma；但如果数据恰好是正常的，我会显示 99.999% 的数据 - 与目标相去甚远。换句话说，我会被缩小 2.25 倍，即显示的区域比必要的多 5 倍，从而使有趣的部分变得很小。如果数据比正常数据高，那就更糟了。所以，当然，这可能是一个绝对的最低限度，但它不是一个很好的近似值。