在循环缓冲区上计算标准差

Question

我需要计算存储在循环缓冲区中的值的标准差。最终算法将在资源受限设备上运行，因此我希望它尽可能轻量级。天真的方法是每次推入新值时重新评估整个缓冲区的标准差，但这样会非常慢。理想情况下，我希望有一种算法，可以在推入新值时动态更新当前标准差的值。 维基百科报告了一些快速计算技术，但它们只能用于流：在我的情况下，当推入新值时，应该计算标准差，就好像已弹出的最后一个值从未存在过。

tl;dr：如何以最小的计算量计算循环缓冲区的标准差？

- F. Cali

Answer 1

标准差可以表示为：

 stddev = sqrt(1/N * SUM(x[i]^2) - 1/N^2 * SUM(x[i])^2)

对于未修正的样本标准差，可以写成：

对于修正后的样本标准差，可以写成：

 stddev = sqrt(1/(N-1) * SUM(x[i]^2) - 1/(N^2-N) * SUM(x[i])^2)

如果你维护两个累加器，一个用于SUM(x[i]^2)，另一个用于SUM(x[i])，那么每次更新这些累加器都是很简单的（减去最老的值的影响，加上最新值的影响）。N当然是你缓冲区的长度。

当然，如果你在浮点数中进行此操作，则可能会遇到舍入误差（一般情况下，(x + y) - x != y）。我认为没有任何简单的解决方法。

- Oliver Charlesworth

---

有一些解决方案。可以使用Kahan求和算法来减轻问题。这将消除由于删除一个元素而产生的舍入误差。还有其他需要更多存储空间的方法。 - Alexandre C.

@Alexandre：我指的是真正意义上的修复。Kahan求和可以减少问题，但不能消除它。如果你知道零误差修复方法，我会非常感兴趣阅读相关内容。 - Oliver Charlesworth

如果卡汉算法不够用，还有舒克算法可供选择。但要注意空间需求，因为资源受到限制。卡汉算法每个总和只需要一个额外的变量。 - Alexandre C.

这里有一个（但是O(n^2)时间和n空间）的算法！这个问题讨论了Shewchuk的算法。其思想是存储所有的部分和（Kahan只存储最后一个部分和）。但是它需要额外的n空间。在我看来，Kahan是一个很好的折衷方案，我在这里全力推荐它。 - Alexandre C.

无论如何，http://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/ 是一篇很棒的文章。我可能对n^2运行时间有所错误，这值得检查。不管怎样，这是一个有趣的问题。 - Alexandre C.

Answer 2

保留三个数字的集合：值的数量、值的总和以及值的平方总和。为方便起见，我们将它们称为k、sn和sn2。

如果像我认为你在暗示的那样，一个新值总是替换环形队列中的旧值，那么也许数量是恒定的。或者可能存在不满的队列。无论哪种方式：

每次向队列添加值时，将计数加1、将该值加到总和中，并将该值的平方加到平方总和中。也就是说，如果添加一个新值“n”，则k=k+1，sn=sn+n，而sn2=sn2+n^2。

每次从队列中删除值时，将计数减1、从总和中减去该值，并从平方总和中减去该值的平方。也就是说，如果删除一个值“n”，则k=k-1，sn=sn-n，而sn2=sn2-n^2。

这样，您可以在没有重新计算总和的情况下轻松地在每次更改后重新计算标准偏差。

请注意，这意味着您必须能够在实际删除之前“捕获”一个值。

注：我怀疑这就是我的口袋计算器是如何做到的，因为它有转储sum(n)和sum(n²)的功能，我可以“删除”我从未添加过的值，比如我可以说将2和4添加到集合中，然后删除3，并且它会说好的没问题。因此，我认为它没有保留列表：它必须只是保留计数和总和。

- Jay

Answer 3

最简单的方法是反转 Knuth 的公式（来自关于方差计算的维基百科文章：Wikipedia article：

M_k-1 = M_k - (x_k - M_k)/(k-1)

S_k-1 = S_k - (x_k - M_k-1)(x_k - M_k)

但是，请注意，浮点误差将在整个过程中累积！这意味着您的平均数和方差数字将倾向于漂移，而此方法不能作为准确的在线方法使用。

一种稳定的在线方法，每个样本需要 O(log N) 操作（其中 N 是队列中元素的数量）是使用统计项的二进制合并树，使用维基百科文章中的“并行合并”公式，如下所示（以 C++ 形式）：

struct Statistic {
  int k;
  Element M;
  Element S;

  Statistic(Element x)
    : k(1)
    , M(x)
    , S(0)
  {}
  Statistic(Statistic a, Statistic b)
    : k(a.k + b.k)
    , M(a.M*a.k + b.M*b.k)/float(k)
    , S(a.S + b.S + (a.M-b.M)*(a.M-b.M)*(a.k*b.k/float(k)))
  {}
};

为了实现一个稳定的O(log N)在线算法，需要维护一个平衡的二叉树来存储上述统计信息，其中叶子节点代表单个元素；根节点将提供所需的在线统计信息。当您以旋转缓冲区的方式更新元素时，每次从叶子节点到根节点的更新传播将需要O(log N)次操作。

- comingstorm

---

不错！还可以看看这里的类似想法：http://blog.sigfpe.com/2010/11/statistical-fingertrees.html - Alexandre C.