滚动方差算法

Question

我正在尝试寻找一种高效、数值稳定的算法来计算滚动方差（例如，一个20周期的滚动窗口内的方差）。我知道Welford算法可以有效地计算数字流的运行方差（只需要一次遍历），但不确定是否可以将其适应于滚动窗口。我还希望解决方案避免John D. Cook在这篇文章顶部讨论的精度问题。任何语言的解决方案都可以。

- Abiel

---

+1 因为提到了 Welford 算法；我知道它在 Knuth 的书里，但从未知道它的原始来源。 - Jason S

你好，你最终做了什么？你采用了Chan的算法吗？顺便说一下，“卡汉求和”在使用“朴素”方法（跟踪值及其平方的总和）时，不应该能够克服数值不稳定性吗？ - Arthur

另一个选择是指数加权移动方差，它将产生与简单移动平均不同的值，但不需要循环缓冲区，因此更节省内存。 - impopularGuy

Answer 1

我一直在处理相同的问题。

计算平均数可以采用迭代的方法，但需要在循环缓冲区中保留完整的数值历史记录。

next_index = (index + 1) % window_size;    // oldest x value is at next_index, wrapping if necessary.

new_mean = mean + (x_new - xs[next_index])/window_size;

我改编了Welford算法，并且对我测试过的所有值都有效。

varSum = var_sum + (x_new - mean) * (x_new - new_mean) - (xs[next_index] - mean) * (xs[next_index] - new_mean);

xs[next_index] = x_new;
index = next_index;

要获得当前方差，请将varSum除以窗口大小：variance = varSum / window_size;

- DanS

---

在更新varSum时，使用varSum += (x_new + x_old - mean - new_mean) * (x_new - x_old)可能会更加稳定，其中x_old = xs[next_index]，因为您从两个相减的项目中去除了一个有可能很大的mean * new_mean项。 除此之外，在这里这是最正确的答案，非常遗憾它没有得到更多的认可。 - Jaime

为了澄清Jaime的回答，他进行了一些代数运算，采用了DanS的varSum方程并分配了乘法。一些项会被取消，但你还需要使用x_new * x_old - x_new * x_old的技巧来得出他的结果。 - Ryan J McCall

非常晚的评论：为什么你要除以 window_size 而不是 window_size-1。换句话说：为什么你不使用贝塞尔校正。我注意到 John D. Cook 在他的运行方差代码中包括了贝塞尔校正。 - hansfn

你能不能通过 variance += (x_new + x_old - mean - new_mean) * (x_new - x_old) / window_size 来完全移除 varSum？ - Guiorgy

Answer 2

我也遇到了这个问题。 有一些关于计算累积方差的很好的帖子，例如John Cooke的 准确计算运行方差 帖子和Digital explorations的帖子，Python代码用于计算样本和总体方差、协方差和相关系数。只是找不到适用于滚动窗口的任何内容。

Subluminal Messages的Running Standard Deviations 帖子对于使滚动窗口公式工作至关重要。Jim取值与Welford方法使用均值的平方差的总和不同，而是取值的平方差的幂和。 公式如下：

 

今天PSA = 昨天PSA + (((x今天 * x今天)-x昨天)) / n

    

 
• x =时间序列中的值
 
• n =您迄今分析的值的数量。
 

但是，要将Power Sum Average公式转换为窗口化的版本，您需要将公式调整为以下内容：

 

今天PSA = 昨天PSA + (((x今天 * x今天) - (x昨天 * x昨天) / n

    

 
• x =时间序列中的值
 
• n =用于滚动窗口的期间。
 

您还需要滚动简单移动平均公式：

 

今天SMA = 昨天SMA + ((x今天-x今天-n) / n

    

 
• x =时间序列中的值
 
• n =用于您的滚动窗口的期间。
 

从那里，您可以计算Rolling Population Variance：

 

Population Var今天=（PSA今天*n-n*SMA今天*SMA今天）/ n

或者使用滚动样本方差:

今天的样本方差=（今天的PSA * n - n * 今天的SMA * 今天的SMA）/（n-1）

我在几年前的博客文章中涵盖了这个主题，并提供了示例Python代码，Running Variance。

希望这可以帮到你。

请注意：我为这个答案提供了所有博客文章和LaTeX数学公式的链接（图片）。但由于我的声望太低（＜10），我只能提供2个超链接和绝对没有图片。对此我深感抱歉，希望这不会影响内容。

- Mike Taylor

---

在这个公式中：Population Var today = (PSA today * n - n * SMA today * SMA today) / n，为什么不去掉 n？Population Var today = (PSA today - SMA today * SMA today)。 - astef

由于在公式中对样本进行了平方处理，这个算法展示了正是 OP 试图避免的数值不准确性。 - marton78

好的，我会尽力进行翻译。这句话的意思是：“是的，这不是一种数值稳定的方法。最接近正确答案的是@DanS下面的回答。” - Jaime

感谢您的解释，这是一个C＃实现 https://gist.github.com/mattdot/d459b1cb15480fefd953841a1ac70be8 - Matt Dotson

Answer 3

如果你更喜欢代码而不是文字（基于DanS的文章）：http://calcandstuff.blogspot.se/2014/02/rolling-variance-calculation.html

本文介绍了滚动方差计算的方法，通过使用代码进行计算。请参考以上链接获取更多信息。

public IEnumerable RollingSampleVariance(IEnumerable data, int sampleSize)
{
    double mean = 0;
    double accVar = 0;

    int n = 0;
    var queue = new Queue(sampleSize);

    foreach(var observation in data)
    {
        queue.Enqueue(observation);
        if (n < sampleSize)
        {
            // Calculating first variance
            n++;
            double delta = observation - mean;
            mean += delta / n;
            accVar += delta * (observation - mean);
        }
        else
        {
            // Adjusting variance
            double then = queue.Dequeue();
            double prevMean = mean;
            mean += (observation - then) / sampleSize;
            accVar += (observation - prevMean) * (observation - mean) - (then - prevMean) * (then - mean);
        }

        if (n == sampleSize)
            yield return accVar / (sampleSize - 1);
    }
}

- Joachim

Answer 4

实际上，据我所知，Welford算法可以轻松地适应于计算加权方差。通过将权重设置为-1，您应该能够有效地取消掉元素。我还没有检查过数学是否允许负权重，但初步看应该可以！

我使用ELKI进行了一个小实验：

void testSlidingWindowVariance() {
MeanVariance mv = new MeanVariance(); // ELKI implementation of weighted Welford!
MeanVariance mc = new MeanVariance(); // Control.

Random r = new Random();
double[] data = new double[1000];
for (int i = 0; i < data.length; i++) {
  data[i] = r.nextDouble();
}

// Pre-roll:
for (int i = 0; i < 10; i++) {
  mv.put(data[i]);
}
// Compare to window approach
for (int i = 10; i < data.length; i++) {
  mv.put(data[i-10], -1.); // Remove
  mv.put(data[i]);
  mc.reset(); // Reset statistics
  for (int j = i - 9; j <= i; j++) {
    mc.put(data[j]);
  }
  assertEquals("Variance does not agree.", mv.getSampleVariance(),
    mc.getSampleVariance(), 1e-14);
}
}

我相对于精确的两遍算法获得了约14位数字的精度；这大约是双倍精度所能期望的。请注意，由于额外的除法，Welford确实会带来一些计算成本——它需要比精确的两遍算法多大约两倍的时间。如果您的窗口大小很小，实际上重新计算平均值，然后在第二次通过方差每次进行可能更加明智。

我已将此实验作为单元测试添加到ELKI中，您可以在此处查看完整源代码：http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki/math/TestSlidingVariance.java 它还与精确的两遍方差进行比较。

但是，在偏斜数据集上，行为可能会有所不同。这个数据集显然是均匀分布的；但我也尝试过排序数组，并且它也有效。 更新：我们发表了一篇详细介绍不同加权方案的（协）方差的论文：

Schubert, Erich和Michael Gertz. "数值稳定的(协)方差并行计算。" 第30届科学和统计数据库管理国际会议论文集。ACM，2018年。(获得SSDBM最佳论文奖)。此外，本文还讨论了如何使用加权方法进行计算的并行化，例如AVX、GPU或集群。

- Erich Schubert

---

将ELKI MeanVariance.java类移植到JS，添加了值缓冲区，并使用-1的权重来删除值。我发现结果精度取决于通过累加器运行多少个值。在通过1-10M个值之后，我看到了约12位数字的精度。（即“足够好”）感谢您提供使用-1权重的提示！ - broofa

如果您需要比这更高的精度，您可能需要使用Kahan求和或Shewchuk算法。它们使用额外的浮点数来存储丢失的数字，因此可以提供更高的精度。但是实现会变得更加混乱和缓慢。有关详细信息，请参见我在帖子中添加的参考文献。 - Erich Schubert

Answer 5

我知道这个问题很久了，但如果还有其他人感兴趣，这里是Python代码。它受到johndcook博客文章、@Joachim的、@DanS的代码和@Jaime评论的启发。下面的代码仍然会在小数据窗口大小时产生一些不精确性。祝使用愉快。

from __future__ import division
import collections
import math


class RunningStats:
    def __init__(self, WIN_SIZE=20):
        self.n = 0
        self.mean = 0
        self.run_var = 0
        self.WIN_SIZE = WIN_SIZE

        self.windows = collections.deque(maxlen=WIN_SIZE)

    def clear(self):
        self.n = 0
        self.windows.clear()

    def push(self, x):

        self.windows.append(x)

        if self.n <= self.WIN_SIZE:
            # Calculating first variance
            self.n += 1
            delta = x - self.mean
            self.mean += delta / self.n
            self.run_var += delta * (x - self.mean)
        else:
            # Adjusting variance
            x_removed = self.windows.popleft()
            old_m = self.mean
            self.mean += (x - x_removed) / self.WIN_SIZE
            self.run_var += (x + x_removed - old_m - self.mean) * (x - x_removed)

    def get_mean(self):
        return self.mean if self.n else 0.0

    def get_var(self):
        return self.run_var / (self.WIN_SIZE - 1) if self.n > 1 else 0.0

    def get_std(self):
        return math.sqrt(self.get_var())

    def get_all(self):
        return list(self.windows)

    def __str__(self):
        return "Current window values: {}".format(list(self.windows))

- ewerlopes

---

感谢您提供的Python综合想法。我不喜欢在进入else块时窗口大小变为WIN_SIZE - 1。因此，如果在调用push时WIN_SIZE为10并且我们进行了附加，则由于使用了deque构造函数选项，它仍然为10，然后在else块中popleft会进一步减小大小到9。所以也许可以使用maxlen=WIN_SIZE + 1？或者不使用maxlen选项。此外，可以放弃n变量并使用len(self.windows)。 - Ryan J McCall

在 get_var 方法中，分母应该是 self.n 或者 len(self.windows)。 - Ryan J McCall

Answer 6

这里有一种分治方法，可以实现时间复杂度为 O(log k) 的更新，其中 k 是样本数。由于与成对求和和快速傅里叶变换（FFT）相同的原因，它应该是相对稳定的，但它有点复杂，常数并不好。

假设我们有一个长度为 m、均值为 E(A)、方差为 V(A) 的序列 A，以及一个长度为 n、均值为 E(B)、方差为 V(B) 的序列 B。让 C 成为 A 和 B 的串联。我们有：

p = m / (m + n)
q = n / (m + n)
E(C) = p * E(A) + q * E(B)
V(C) = p * (V(A) + (E(A) + E(C)) * (E(A) - E(C))) + q * (V(B) + (E(B) + E(C)) * (E(B) - E(C)))

现在，将元素放入红黑树中，每个节点都装饰有以该节点为根的子树的平均值和方差。向右插入；向左删除。（由于我们只访问端点，因此伸展树可能是O(1) 摊销的，但我猜测摊销对于您的应用程序来说是一个问题。）如果k在编译时已知，您可能可以像FFTW一样展开内部循环。

- userOVER9000

---

（注意：除非k非常大，否则可以计算q = 1 - p。） - userOVER9000

好的，这基本上是Chan等人在维基百科上描述的并行算法。这就是我没有向下滚动得到的结果... - userOVER9000

你能否稍微详细地解释一下如何将这个算法应用于移动窗口的方差计算？我对Chan等人的方法有些了解，但是认为它是一种单次通过方法，用于计算整个样本的单个方差，并且具有问题可以分成并行运行的部分的额外优势。 - Abiel

Chan等人提供了一种计算部分串联的统计数据的方法，给定部分的统计数据。高层次的想法是维护一个部分集合（实际上只是它们的统计数据），使得任何窗口都是O(log k)个部分的串联。一种方法是使用平衡二叉树，但正如Rex所指出的那样，这是过度的，我们可以仅维护大小为2的幂次方的对齐部分的统计数据（例如[0,1)，[1,2)，[0,2)，[2,3)，[3,4)，[2,4)，[0,4)等）。 - userOVER9000

Answer 7

这只是对DanS提供的出色答案的小补充。以下方程式用于从窗口中删除最旧的样本并更新均值和方差。这很有用，例如，如果您想在输入数据流的右侧附近采取较小的窗口（即仅删除最旧的窗口样本而不添加新样本）。

window_size -= 1; % decrease window size by 1 sample
new_mean = prev_mean + (prev_mean - x_old) / window_size
varSum = varSum - (prev_mean - x_old) * (new_mean - x_old)

这里，x_old是您希望移除的窗口中最旧的样本。

- vibe

Answer 8

我猜跟踪你的20个样本、计算1..20中X的平方和、计算1..20中X的和，然后在每次迭代中连续重新计算这两个总和不够高效？实际上可以在不每次加起来、平方等所有样本的情况下重新计算新的方差。

具体操作：

Sum(X^2 from 2..21) = Sum(X^2 from 1..20) - X_1^2 + X_21^2
Sum(X from 2..21) = Sum(X from 1..20) - X_1 + X_21

- John

---

我认为这个解决方案容易受到我在原帖中提到的稳定性问题的影响（http://www.johndcook.com/standard_deviation.html）。特别是当输入值很大且它们之间的差异很小时，结果实际上可能是负数。由于我无法控制输入，因此我更愿意避免采用这种方法。 - Abiel

哦，我明白了。你能说一些关于输入的事情吗？预期使用方式？这是一个可以通过增加更多位数来解决的问题吗（64位浮点数、任意精度算术等）？如果你在有效数字方面超越输入，舍入误差会消失，不是吗？ - John

同意 - 这会有稳定性问题。想象一下1000个接近1,000,000.0的样本，然后是20个接近零的样本。 - Jason S

@Jason S：滚动方差就是它的本质。在从100万到~零的过渡中可能会发生很多事情，但这就是问题的本质。此外，在变化发生时，前980个100万值中的第一个已经不重要了。我的评论建议如果你的计算中有足够的有效数字，那么这些都不重要。 - John

输入可以是任何东西。值的数量级肯定可以达到万亿级别，虽然原始数据只有几个小数点的精度，但用户可以在计算方差之前转换他们的数据（例如除以任何标量）。 - Abiel

Answer 9

这是另一种O(log k)的解决方案：先找到原始序列的平方，然后对成对、四个一组等进行求和。（您需要一些缓冲区来有效地找到所有这些。）然后将您需要的这些值相加以得出答案。例如：

|||||||||||||||||||||||||  // Squares
| | | | | | | | | | | | |  // Sum of squares for pairs
|   |   |   |   |   |   |  // Pairs of pairs
|       |       |       |  // (etc.)
|               |
   ^------------------^    // Want these 20, which you can get with
        |       |          // one...
    |   |       |   |      // two, three...
                    | |    // four...
   ||                      // five stored values.

现在你可以使用标准的E(x^2)-E(x)^2公式，然后就完成了。（如果您需要处理少量数字的稳定性较好，这样做可能不太可靠；这仅适用于滚动误差的累积导致问题的情况。）

话虽如此，在大多数体系结构上，对20个平方数求和是非常快速的。如果您要进行更多操作——比如说，几百次——更高效的方法显然更好。但我不确定这里是否采用蛮力法才是正确的方式。

- Rex Kerr

---

不要使用标准的E(x^2)-E(x)^2公式，它完全不稳定。使用更好的算法进行调整。 - userOVER9000

@userOVER9000 - 你为什么担心超过20个项目的稳定性？在数百万条目中累积的累积误差是一个问题（特别是在制作滚动窗口时），但这不是这里的问题。 - Rex Kerr

我对此感到担忧，因为这是一个问题。去阅读维基百科的文章，如果你仍然不信，请尝试计算20个iid样本的N（1，1e-10）的方差。 - userOVER9000

我实际上还没有看到任何具有合理单位和原点的现实数据集出现这个问题，但如果这是 OP 想要的，那也没什么问题... - Rex Kerr

Answer 10

我希望我是错的，但我认为这不能“快速”完成。话虽如此，计算的大部分是跟踪窗口内的期望值，这可以轻松完成。

我留下一个问题：你确定你需要一个窗口函数吗？除非你正在使用非常大的窗口，否则最好只使用一个众所周知的预定义算法。

- Andrew White