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

Posted 2023-02-16

【中文标题】用于估计统计中位数、众数、偏度、峰度的“在线”（迭代器）算法？

【英文标题】："On-line" (iterator) algorithms for estimating statistical median, mode, skewness, kurtosis?

【发布时间】：2010-11-06 17:25:54

【问题描述】：

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

我想计算基本统计数据：

平均值：算术平均值 方差：与均值的平方偏差的平均值 标准差：方差的平方根 中位数：将较大一半数字与较小一半数字分开的值 mode：在集合中找到的最频繁的值 偏度：tl;博士 峰度：tl；博士

计算任何这些的基本公式是小学算术，我知道它们。还有许多实现它们的统计库。

我的问题是我正在处理的集合中有大量（数十亿）值：在 Python 中工作，我不能只用数十亿个元素制作一个列表或散列。即使我用 C 写了这个，十亿元素的数组也不太实用。

数据未排序。它是由其他过程随机、即时生成的。每个集合的大小是高度可变的，并且大小不会提前知道。

我已经想出了如何很好地处理均值和方差，以任意顺序遍历集合中的每个值。 （实际上，在我的例子中，我按照它们生成的顺序来处理它们。）这是我正在使用的算法，礼貌http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm：

初始化三个变量：count、sum 和 sum_of_squares 对于每个值： 递增计数。 将该值相加。 将值的平方添加到 sum_of_squares。 将总和除以计数，存储为变量均值。 将 sum_of_squares 除以计数，存储为变量 mean_of_squares。 均方，存储为 square_of_mean。 从 mean_of_squares 中减去 square_of_mean，存储为方差。 输出均值和方差。

这种“在线”算法有弱点（例如，准确性问题，因为 sum_of_squares 迅速增长到大于整数范围或浮点精度），但它基本上可以满足我的需要，而无需存储每个集合中的每个值。

但我不知道是否存在用于估计附加统计数据（中位数、众数、偏度、峰度）的类似技术。只要处理 N 个值所需的内存大大小于 O(N)，我就可以接受有偏差的估计器，甚至可以在一定程度上损害准确性的方法。

如果该库具有“在线”计算这些操作中的一项或多项的功能，将我指向现有的统计库也会有所帮助。

【问题讨论】：

传入的数据会不会排序，会提前知道输入的个数吗？

*** 上有用的现有链接：***.com/questions/895929/…

是整数数据还是浮点数据？你有最大值还是最小值？

dmckee：我实际上是在使用 Welford 方法来计算标准差。但我在该链接中没有看到任何关于众数、中位数、峰度或偏度的信息……我错过了什么吗？

stephan：有些数据集是整数，有些是浮点数。总体分布非常接近正态分布（高斯分布），因此我们可以建立置信区间，但没有硬性范围边界（在某些情况下 x > 0 除外）。

【参考方案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）

我刚刚修改了增量中位数估计器来估计任意分位数。通常，quantile function 告诉您将数据分成两个部分的值：p 和 1-p。下面逐步估计这个值：

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

值 p 应在 [0,1] 范围内。这实质上将 sgn() 函数的对称输出 -1,0,1 向一侧倾斜，将数据样本分成两个大小不等的 bin（分数 p 和 1-p数据分别小于/大于分位数估计）。请注意，对于 p=0.5，这会简化为中值估计量。

更新（2021-11-19）

有关此处描述的中值估计器的更多详细信息，我想重点介绍以下 cmets 中链接的这篇论文：Bylander & Rosen, 1997, A Perceptron-Like Online Algorithm for Tracking the Median。这是作者网站上的postscript version。

【讨论】：

这个中值估计器很棒。你知道 0.25/0.75 分位数是否有类似的估计器？

@Gacek，当然：将输入流拆分为 Lohalf median，并在每一半上使用 running-median。

@Gacek：我刚刚用增量方法更新了我的答案来估计任何分位数，您可以将 p 设置为 0.25、0.75 或 [0,1 内的 any 值].

这对均值很有用，但我没有看到它如何产生任何接近中位数的东西。以毫秒时间戳序列为例：[1328083200000, 981014400000, -628444800000, 318240000000, 949392000000]，其中位数为318240000000。此等式将先前的中位数移动 +/- eta，其中推荐值为 0.001。对于像这样的大数字，这不会有任何作用，对于非常小的数字来说，它可能太大了。你会如何选择一个eta，它实际上给了你正确的答案而不知道先验的答案？

想象数字有单位，例如毫米。然后很明显 eta （用于中位数的估计）必须与测量具有相同的单位，因此像 0.001 这样的通用值根本没有任何意义。一个看似更好的方法是从绝对偏差的运行估计中设置 eta：对于每个新值 sample，更新 cumadev += abs(sample-median)。然后设置eta = 1.5*cumadev/(k*k)，其中k是目前看到的样本数。

【参考方案2】：

偏度和峰度

对于偏度和峰度的在线算法（沿着方差线），请参阅同一 wiki 页面here 更高矩统计的并行算法。

中位数

没有排序的数据，中位数很难。如果你知道你有多少数据点，理论上你只需要部分排序，例如通过使用selection algorithm。但是，这对数十亿的价值并没有太大帮助。我建议使用频率计数，请参阅下一节。

具有频率计数的中位数和众数

如果是整数，我会数 frequencies，可能会在我确信它不再相关的某个值之外切断最高和最低值。对于浮点数（或太多整数），我可能会创建桶/间隔，然后使用与整数相同的方法。 （近似）众数和中位数计算比基于频率表更容易。

正态分布随机变量

如果它是正态分布的，我将使用人口样本 mean、variance、skewness 和 kurtosis 作为一小部分子集的最大似然估计量。计算这些的（在线）算法，你现在已经知道了。例如。读取数十万或数百万个数据点，直到您的估计误差变得足够小。只需确保您从集合中随机选择（例如，您不会通过选择前 100'000 个值来引入偏差）。同样的方法也可用于估计正常情况下的众数和中位数（因为样本均值都是估计量）。

更多课程

如果有帮助，上述所有算法都可以并行运行（包括许多排序和选择算法，例如 QuickSort 和 QuickSelect）。

我一直假设（除了关于正态分布的部分）我们谈论的是样本矩、中位数和众数，而不是给定已知分布的理论矩的估计量。

一般来说，考虑到数据量，对数据进行采样（即仅查看子集）应该是相当成功的，只要所有观察结果都是相同随机变量的实现（具有相同的分布）和矩, 众数和中位数对于这个分布实际上是存在的。最后一个警告并非无害。例如，Cauchy Distribution 的均值（以及所有更高的矩）不存在。在这种情况下，“小”子集的样本均值可能与整个样本的样本均值相差很大。

【参考方案3】：

我在一个名为LiveStats 的简洁Python 模块中实现了P-Square Algorithm for Dynamic Calculation of Quantiles and Histograms without Storing Observations。它应该可以非常有效地解决您的问题。该库支持您提到的所有统计信息，但模式除外。我还没有找到一个令人满意的模式估计解决方案。

【讨论】：

仅供参考：p-square 算法在 C++ 中提升：<boost/accumulators/statistics/weighted_p_square_cumul_dist.hpp>.

【参考方案4】：

Ryan，恐怕你没有正确计算均值和方差...这是几周前出现的 here。在线版本（实际上以 Welford 方法的名称命名）的优点之一是它特别准确和稳定，请参阅讨论 here。优点之一是您不需要存储总和或总平方和......

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

【讨论】：

回复：skewness and kurtosis 是的。见这篇文章：johndcook.com/blog/skewness_kurtosis

【参考方案5】：

问题中引用的***文章包含在线计算偏度和峰度的公式。

对于模式 - 我相信 - 没有办法在线执行此操作。为什么？假设您输入的所有值都不同，除了与前一个重复的最后一个值。在这种情况下，您必须记住在输入中已经看到的所有值，以检测最后一个值是否与之前看到的值重复并使其成为最常见的值。

对于中位数几乎相同 - 直到最后一个输入，如果所有输入值都不同，您不知道哪个值将成为中位数，因为它可能在当前中位数之前或之后。如果您知道输入的长度，您可以在不将所有值存储在内存中的情况下找到中值，但您仍然需要存储其中的许多（我猜大约是一半），因为错误的输入序列可能会在后半部分可能使前半部分的任何值成为中位数。

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

【参考方案6】：

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

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

分位数估计示例：http://www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014

模式估计示例：Bickel DR。连续数据的模式和偏度的稳健估计。计算统计和数据分析。 2002；39：153–163。 doi: 10.1016/S0167-9473(01)00057-3.

这些是计算统计的活跃领域。您正在进入的领域没有任何单一的最佳精确算法，而是它们的多样性（实际上是统计估计器），它们具有不同的属性、假设和性能。是实验数学。可能有成百上千篇关于这个主题的论文。

最后一个问题是您是否真的需要偏度和峰度本身，或者更可能需要一些其他参数，这些参数在表征概率分布时可能更可靠（假设您有一个概率分布！）。你期待高斯吗？

您是否有办法清理/预处理数据以使其大部分为高斯分布？ （例如，金融交易金额在取对数后通常呈高斯分布）。你期望有限的标准偏差吗？你期待肥尾巴吗？您关心的数量是尾部还是散装？

【参考方案7】：

每个人都在说您不能以在线方式进行该模式，但事实并非如此。这是一个article，描述了一个算法来解决这个问题，该算法由耶鲁大学的 Michael E. Fischer 和 Steven L. Salzberg 于 1982 年发明。来自文章：

多数查找算法使用它的一个寄存器来临时 存储流中的单个项目；这个项目是当前的 多数元素的候选人。第二个寄存器是一个计数器 初始化为0。对于流的每个元素，我们询问算法 执行以下例程。如果计数器读数为 0，请安装 当前流元素作为新的多数候选者（取代任何 可能已经在寄存器中的其他元素）。那么，如果 当前元素匹配多数候选者，增加计数器； 否则，递减计数器。在循环的这一点上，如果 到目前为止看到的流的一部分具有多数元素，该元素是 在候选寄存器中，并且计数器的值大于 0. 如果没有多数派怎么办？无需对数据进行第二次传递——这在流环境中是不可能的—— 算法不能总是给出明确的答案 环境。它只是承诺正确识别大多数 元素，如果有的话。

它也可以扩展以找到具有更多内存的前 N ​​但这应该可以解决模式。

【讨论】：

这是一个有趣的算法，但除非我遗漏了什么，虽然所有多数值都是众数，但并非所有众数都是众数。

链接已失效，所以我很高兴包含说明。但是，如上所述，仅当多数候选第二次出现与第一次出现相邻时，计数器才会增加。这意味着对数据进行了排序。在在线（流式传输）数据案例中不能保证这一点。对于随机排序的数据，这不太可能找到任何模式。

【参考方案8】：

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

也就是说，除非您正在处理某种病态情况，否则补救措施（Rousseuw 和 Bassett 1990）可能就足以满足您的目的。

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

【参考方案9】：

中值和众数不能仅使用可用的常量空间在线计算。但是，由于中位数和众数无论如何都比“定量”更具“描述性”，因此您可以估计它们，例如通过对数据集进行采样。

如果数据从长远来看是正态分布的，那么你可以只用你的平均值来估计中位数。

您还可以使用以下技术估计中值：为数据流中的每 1,000,000 个条目建立一个中值估计 M[i]，以便 M[0] 是前一百万个条目的中值，M[ 1] 第二个一百万条目的中值等。然后使用 M[0]...M[k] 的中值作为中值估计量。这当然会节省空间，您可以通过“调整”参数 1,000,000 来控制要使用多少空间。这也可以递归地推广。

【参考方案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) 您将这些指向您的功能。

另外，看看 Pearson 的近似值。在如此大的数据集上，它会非常相似。 3（平均值 - 中位数）/标准差 您的中位数为 max - min/2

对于浮动模式没有意义。人们通常会将它们放在一个非常大的容器中（例如 1/100 * (max - min)）。

【参考方案11】：

Pebay 等人解决了这个问题：

https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf

【参考方案12】：

中位数

可以在此处找到最近的两个百分位近似算法及其 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 还提供相对错误保证。

【参考方案13】：

我倾向于使用可以自适应的存储桶。桶大小应该是您需要的精度。然后随着每个数据点的进入，您将一个添加到相关存储桶的计数中。 这些应该为您提供中值和峰度的简单近似值，方法是将每个桶计算为其计数加权的值。

一个问题可能是数十亿次运算后浮点分辨率的损失，即加一不会再改变值！为了解决这个问题，如果最大存储桶大小超过某个限制，您可以从所有计数中取出一个很大的数字。

【参考方案14】：

for j in range (1,M):
    y=np.zeros(M) # build the vector y
    y[0]=y0

    #generate the white noise
    eps=npr.randn(M-1)*np.sqrt(var)

    #increment the y vector
    for k in range(1,T):
        y[k]=corr*y[k-1]+eps[k-1]

    yy[j]=y

list.append(y)

【讨论】：

可以使用一些解释来更好地将其与原始问题联系起来。