估算统计中位数、众数、偏度和峰度的“在线”（迭代器）算法？

Question

有没有一种算法可以估计一组值的中位数、众数、偏度和/或峰度，但不需要一次性将所有值存储在内存中?

我想计算基本统计量：

• 平均值: 算术平均值
• 方差: 平均偏差的平方
• 标准差: 方差的平方根
• 中位数: 将数字中较大的一半与较小的一半分开的值
• 众数: 集合中出现最频繁的值
• 偏度: tl; dr
• 峰度: tl; dr

计算任何一个这些基本统计量的公式都很简单，并且我知道它们。许多统计库也实现了它们。

我的问题是我处理的集合中有大量（数十亿）的值：在Python中工作时，我不能只是创建一个包含数十亿个元素的列表或哈希表。即使我用C编写，十亿个元素的数组也不太实用。

数据没有排序。它是由其他进程随机生成的。每个集合的大小高度可变，并且大小不会事先知道。

我已经发现如何很好地处理平均值和方差，可以按任意顺序迭代每个集合中的值。（实际上，在我的情况下，我按它们生成的顺序来取）这是我正在使用的算法，由http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm提供：

• 初始化三个变量：count、sum和sum_of_squares。
• 对于每个值：
  • 增加计数。
  • 将该值添加到总和中。
  • 将该值的平方添加到平方总和中。
• 将总和除以计数，存储为变量mean。

将平方和除以计数，存储为变量mean_of_squares。

对均值进行平方，存储为square_of_mean。

从mean_of_squares中减去square_of_mean，存储为方差variance。

输出平均值和方差。

这种"在线"算法存在缺陷（例如，当sum_of_squares快速增长到超过整数范围或浮点精度时会出现精度问题），但它基本上可以给我需要的东西，而不必存储每个集合中的每个值。

但我不知道是否存在类似的技术来估计其他统计信息（中位数、众数、偏度、峰度）。我可以接受有偏估算器，甚至是在一定程度上牺牲准确性的方法，只要处理N个值所需的内存远远小于O（N）。

如果该库具有用于“在线”计算一个或多个这些操作的函数，则指向现有的统计库也将有所帮助。

- Ryan B. Lynch

---

数据是否会按顺序传递，而且您是否提前知道输入数量？ - chillysapien

Stephan：有些数据集是整数，有些是浮点数。人口分布非常接近正态（高斯）分布，因此我们可以建立置信区间，但没有硬性的范围边界（除了在某些情况下 x > 0）。 - Ryan B. Lynch

@Ryan B. Lynch：啊，你的文本描述了朴素的单遍方差算法，所以我认为你不知道Welform的方法。而且我也无法帮助你解决其他问题。抱歉。 - dmckee --- ex-moderator kitten

这些被称为流算法，顺便说一下 [http://scholar.google.com/scholar?&q=streaming+algorithm]。其中一些很容易，一些可以在任意ε内进行逼近（高概率），而有些则难以逼近。 - ShreevatsaR

关于您对库的请求，我有一个Julia包可以解决这些问题：https://github.com/joshday/OnlineStats.jl。OnlineStats实现了用于均值、方差、偏度/峰度、直方图（从中可以估计分位数）和几种分位数近似算法的在线算法。 - joshday

Answer 1

我使用这些增量/递归的均值和中位数估计算法，它们都使用恒定的存储空间：

mean += eta * (sample - mean)
median += eta * sgn(sample - median)

这里的eta是一个小的学习率参数(例如0.001)，sgn()是符号函数，返回{-1, 0, 1}中的一个。(如果数据是非静态的，并且想要跟踪随时间的变化，请使用常量eta; 否则，对于静态源，可以使用类似于eta=1/n的东西作为均值估计器，其中n是迄今为止看到的样本数...不幸的是，这对于中位数估计器似乎不起作用。)

这种增量均值估计器似乎在各个领域都被使用，例如在无监督神经网络学习规则中，但中位数版本似乎要少得多，尽管它具有优点(对异常值具有鲁棒性)。看起来中位数版本可以在许多应用程序中替代均值估计器。

我希望能看到一个类似形式的增量模估计器...

更新(2011-09-19)

我刚刚修改了增量中位数估计器以估算任意分位数。一般而言，分位函数告诉您将数据分成两个部分的值：p和1-p。以下增量估计此值:

quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)

变量p的取值应该在[0,1]之间。这基本上将sgn()函数输出的对称结果{-1,0,1}偏向一侧，将数据样本分成两个不均匀大小的箱（分别是数据小于/大于分位数估计的比例p和1-p）。请注意，当p=0.5时，这就简化为中位数估计。

更新（2021-11-19）

有关此处描述的中位数估计器的更多详细信息，请参见下面评论中链接的这篇论文：Bylander＆Rosen，1997年，用于跟踪中位数的感知器样式在线算法。作者网站上提供了postscript版本。

- Tyler Streeter

---

这个中位数估计器很棒。你知道是否有类似的估计器可以用于0.25/0.75分位数吗？ - Gacek

@Gacek，是的，请参见http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=614292&tag=1 - Ian Goodfellow

@Gacek：我刚刚更新了我的答案，提供了一种增量方法来估计任何分位数，您可以将p设置为0.25、0.75或[0,1]范围内的任何值。 - Tyler Streeter

这对于平均数非常有效，但我没有看到它如何产生任何与中位数接近的结果。举个例子，考虑一个毫秒级时间戳序列：[1328083200000、981014400000、-628444800000、318240000000、949392000000]，其中的中位数为318240000000。这个公式通过+/- eta来移动先前的中位数，推荐值为0.001。对于像这样的大数字来说，这不会起作用，而对于真正的小数字来说，它可能太大了。如果没有事先知道答案，你该如何选择一个能给出正确答案的eta值呢？ - mckamey

假设数字有单位，例如毫米。那么估计中位数的 eta 必须与测量具有相同的单位，因此类似于 0.001 的一般值根本没有任何意义。一种看似更好的方法是从绝对偏差的运行估计中设置 eta：对于每个新值 sample，更新 cumadev += abs（sample-median）。然后设置 eta = 1.5*cumadev/(k*k)，其中 k 是到目前为止观察到的样本数量。 - tholy

Answer 2

偏度和峰度

关于偏度和峰度的在线算法（类似于方差），请参见同一维基页面here中高阶矩统计量的并行算法。

中位数

没有排序数据的中位数很难求。如果你知道有多少数据点，理论上只需要部分排序，例如使用selection algorithm。然而，对于数十亿个值来说，这并没有太大帮助。我建议使用频率计数，请参见下一节。

使用频率计数的中位数和众数

如果是整数，我会计算frequencies，可能会截断高于和低于某个确定不再相关的值。对于浮点数（或过多的整数），我可能会创建桶/间隔，然后使用与整数相同的方法。基于频率表，近似的模式和中位数计算就变得容易了。

正态分布随机变量

如果数据是正态分布的话，我会使用总体样本均值, 方差, 偏度, 和 峰度 作为小样本的最大似然估计量。你已经知道了计算这些估计量的（在线）算法。例如，读入几十万或几百万个数据点，直到你的估计误差足够小。只需确保从你的集合中随机选择（例如，不要通过选择前100,000个值引入偏差）。同样的方法也可用于正态情况下估计众数和中位数（对于两者，样本均值都是估计量）。 进一步评论 以上所有算法都可以并行运行（包括许多排序和选择算法，例如QuickSort和QuickSelect），如果这有助于解决问题。

我一直假设（除了正态分布部分之外），我们谈论的是样本矩、中位数和众数，而不是在已知分布的情况下给出的理论矩的估计量。

通常情况下，由于数据量足够大，只需对数据进行抽样（即仅查看子集）就可以取得很好的效果，前提是所有观测值都是同一随机变量的实现（具有相同的分布），并且该分布的矩、众数和中位数确实存在。最后一个限制条件并不无关紧要。例如，Cauchy Distribution的平均值（以及所有更高阶矩）不存在。在这种情况下，“小”子集的样本均值可能与整个样本的样本均值相差甚远。

- stephan

Answer 3

我在自己编写的名为LiveStats的Python模块中实现了P-Square算法，用于动态计算分位数和直方图而无需存储观测值。它应该可以相当有效地解决您的问题。该库支持您提到的每个统计量，除了众数。我尚未找到满意的众数估计解决方案。

- Sean

---

FYI：p-square算法在C++ Boost中，路径为<boost/accumulators/statistics/weighted_p_square_cumul_dist.hpp>。 - Neil G

Answer 4

Ryan，恐怕你没有正确计算均值和方差... 几周前在这里讨论过。在线版本（实际上称为Welford方法）的一个强项是其特别准确和稳定，详情请参见此处的讨论。其中一个优点是您不需要存储总和或平方总和...

我想不出任何关于众数和中位数的在线方法，它们似乎需要同时考虑整个列表。但很可能类似于方差和均值的方法也适用于偏度和峰度...

- Jaime

---

关于偏度和峰度的问题。是的。请参考这篇文章：https://www.johndcook.com/blog/skewness_kurtosis/ - Jesse Chisholm

Answer 5

在提问中引用的维基百科文章包含计算偏度和峰度的公式。

至于众数——我相信——没有办法在线计算。为什么？假设您输入的所有值都不同，除了最后一个重复了之前的一个值。在这种情况下，您必须记住已经看到的输入中的所有值，以便检测到最后一个值重复了先前看到的值，并使其成为最频繁出现的值。

对于中位数来说，情况几乎相同——在最后一个输入之前，如果所有输入值都不同，则不知道哪个值将成为中位数，因为它可能在当前中位数之前或之后。如果您知道输入的长度，则可以在不存储所有值的情况下找到中位数，但是仍然必须存储其中许多值（我猜大约一半），因为不良的输入序列可能会严重移动中位数在第二半部分，可能使第一半的任何值成为中位数。

（请注意，我仅参考精确计算。）

- Daniel Brückner

Answer 6

如果您有数十亿个数据点，那么您不太可能需要精确答案，而是接近答案。通常，如果您有数十亿个数据点，则生成它们的基本过程可能遵循某种统计稳定性/遍历性/混合性质。此外，您期望分布是否相对连续也可能很重要。

在这些情况下，存在算法可用于在线、低内存、估计分位数（中位数是0.5分位数的特殊情况）以及众数，如果您不需要精确答案。这是统计学的一个活跃领域。

分位数估计示例：http://www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014 模式估计示例：Bickel DR. Robust estimators of the mode and skewness of continuous data. Computational Statistics and Data Analysis. 2002;39:153-163. doi: 10.1016/S0167-9473(01)00057-3.

这些是计算统计学的活跃领域。您正在进入没有任何单一最佳精确算法的领域，而是有各种性质、假设和性能不同的统计估计量。这是实验数学。可能有数百到数千篇关于此主题的论文。

最后一个问题是，您是否真的需要倾斜度和峰度本身，或者更有可能需要一些其他参数来更可靠地描述概率分布（假设您有概率分布！）。您期望高斯分布吗？

您有清理/预处理数据以使其大多数情况下呈现高斯分布的方法吗？ （例如，金融交易金额在取对数后通常有点高斯）。您期望有限标准差吗？您期望有肥尾巴吗？您关心的数量是在尾部还是在主体中？

- Matt Kennel

Answer 7

大家一直在说在线模式是不可能完成的，但这并不是真的。这里有一篇文章描述了一个算法来解决这个问题，该算法是由耶鲁大学的Michael E. Fischer和Steven L. Salzberg于1982年发明的。从文章中可以得知：

大多数查找算法使用其中一个寄存器作为流中单个项目的临时存储；该项目是当前的大多数元素候选项。第二个寄存器是初始化为0的计数器。对于流的每个元素，我们要求算法执行以下例程。如果计数器读取0，则将当前流元素安装为新的大多数候选项（替换可能已经在寄存器中的任何其他元素）。然后，如果当前元素与大多数候选项匹配，则增加计数器；否则，减少计数器。在循环的这一点上，到目前为止看到的流具有大多数元素，则该元素位于候选寄存器中，并且计数器保持大于0的值。如果没有大多数元素怎么办？在不进行第二次数据遍历的情况下 - 这在流环境中是不可能的 - 该算法不能始终在这种情况下给出明确的答案。它只承诺正确识别大多数元素（如果存在）。

它还可以扩展以使用更多内存查找前N个元素，但这应该解决模式问题。

- hackartist

---

这是一个有趣的算法，但除非我漏掉了什么，虽然所有的众数都将成为主要值，但并非所有的主要值都是众数。 - jkebinger

链接已经失效，所以我很高兴描述被包括在内。但是，正如所描述的那样，只有当大多数候选人的第二次出现与第一次出现相邻时，计数器才会递增。这意味着数据已排序，而在线（流式）数据情况下不能保证。对于随机排序的数据，这不太可能找到任何模式。 - Jesse Chisholm

Answer 8

最终，如果您没有分布的先验参数知识，我认为您必须存储所有值。

话虽如此，除非您处理某种病态情况，否则中位数（Rousseuw和Bassett 1990）可能足以满足您的需求。

非常简单，它涉及计算各个中位数批次的中位数。

- Edward Grace

Answer 9

中位数

这里可以找到两个最近的百分位数近似算法及其Python实现：

t-Digests

• https://arxiv.org/abs/1902.04023
• https://github.com/CamDavidsonPilon/tdigest

DDSketch

• https://arxiv.org/abs/1908.10693
• https://github.com/DataDog/sketches-py

这两种算法都是用于数据分桶的。由于T-Digest在尾部使用较小的分区，因此在极端值处精度更高（并且在中位数附近较弱）。DDSketch还提供相对误差保证。

- rosa b.

Answer 10

好的，伙计，试试这些：

对于C++：

double skew(double* v, unsigned long n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow((v[i] - mu)/sigma, 3);
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

double kurt(double* v, double n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow( ((v[i] - mu[i]) / sigma) , 4) - 3;
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

你说你已经可以计算样本方差（svar）和平均值（avg），那么你可以将它们指向你的函数来完成。

此外，看一下皮尔逊近似公式。在这么大的数据集上，它会非常相似。 3（平均数-中位数）/标准差 你的中位数是最大值-最小值/2

对于浮点数，众数没有意义。通常会将它们放入一个显著大小的箱子中（例如1/100 *（最大值-最小值））。

- peter