滚动方差算法

Posted 2023-03-12

【中文标题】滚动方差算法

【英文标题】：Rolling variance algorithm

【发布时间】：2011-07-06 01:33:56

【问题描述】：

我正在尝试找到一种高效、数值稳定的算法来计算滚动方差（例如，20 周期滚动窗口上的方差）。我知道Welford algorithm 可以有效地计算数字流的运行方差（它只需要一次通过），但不确定这是否可以适用于滚动窗口。我还想要解决方案来避免 John D. Cook 在this article 顶部讨论的准确性问题。任何语言的解决方案都可以。

【问题讨论】：

+1 用于提及 Welford 算法；我知道它在 Knuth，但我不知道原始来源

你好，你最后做了什么？你适应了陈的算法吗？顺便说一句，在使用“幼稚”方法（跟踪值的总和及其平方）时，kahan sum 不应该能够克服数值不稳定性吗？

【参考方案1】：

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

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 = 用于滚动窗口的时间段。

您可以从那里计算滚动人口方差：

今日人口变异 =（今日 PSA * n - n * 今日 SMA * 今日 SMA）/n

或滚动样本方差：

今天的样本 Var = (今天的 PSA * n - n * 今天的 SMA * 今天的 SMA) / (n - 1)

几年前，我在一篇博文Running Variance 中介绍了这个主题以及示例 Python 代码。

希望这会有所帮助。

请注意：我提供了所有博客文章和数学公式的链接 在乳胶（图像）中为这个答案。但是，由于我的声誉低（

【讨论】：

在这个公式中：Population Var today = (PSA today * n - n * SMA today * SMA today) / n - 为什么不删除n？ Population Var today = (PSA today - SMA today * SMA today).

由于公式中的样本平方，该算法表现出 OP 试图避免的非常数值的不准确性。

是的，这不是一种数值稳定的方法。最接近正确答案的是下面的@DanS。

感谢解释，这是一个C#实现gist.github.com/mattdot/d459b1cb15480fefd953841a1ac70be8

【参考方案2】：

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

平均值很容易迭代计算，但您需要将值的完整历史记录保存在循环缓冲区中。

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;

【讨论】：

执行varSum += (x_new + x_old - mean - new_mean) * (x_new - x_old)（其中x_old = xs[next_index]）可能会稍微稳定一些，因为您从减去更新varSum 的两项中删除了一个可能很大的mean * new_mean 和。除此之外，这是这里最正确的答案，可惜没有得到更多的爱。

为了澄清 Jaime 的答案，他用 DanS 的 varSum 方程进行了一些代数运算并分配了乘法。一些条款取消，但您还必须执行添加 x_new * x_old - x_new * x_old 的技巧才能得出他的结果

很晚的评论：你为什么要跳到window_size 而不是window_size-1。换句话说：你为什么不使用贝塞尔的校正。我注意到 John D. Cook 确实在他的运行方差代码中包含了 Bessel 的校正。

你不能只用variance += (x_new + x_old - mean - new_mean) * (x_new - x_old) / window_size 完全删除varSum 吗？

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

【参考方案4】：

实际上，Welfords 算法可以很容易地适应 AFAICT 以计算 加权 方差。 通过将权重设置为 -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 它还与精确的两遍方差进行比较。

但是，在倾斜的数据集上，行为可能会有所不同。这个数据集显然是均匀分布的；但我也尝试了一个排序数组并且它有效。

更新：我们发表了一篇论文，详细介绍了（共）方差的不同加权方案：

舒伯特、埃里希和迈克尔·格茨。 “（协）方差的数值稳定并行计算。”第 30 届科学与统计数据库管理国际会议论文集。 ACM，2018。（获得 SSDBM 最佳论文奖。）

本文还讨论了如何使用加权来并行化计算，例如，使用 AVX、GPU 或在集群上。

【讨论】：

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

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

【参考方案5】：

这是一种分而治之的方法，它具有O(log k)-time 更新，其中k 是样本数。它应该是相对稳定的，原因与成对求和和 FFT 稳定的原因相同，但它有点复杂，常数也不是很大。

假设我们有一个长度为m的序列A，其均值E(A)和方差V(A)，以及一个长度n的序列B，均值E(B)和方差V(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 样式。

【讨论】：

（注意：计算 q = 1 - p 没问题，除非 k 非常大。）

好的，这基本上是 Chan 等人在 Wikipedia 上描述的并行算法。这就是我不向下滚动所得到的......

您能否更详细地解释一下如何将此算法应用于移动窗口上的方差？我对 Chan 等人的方法稍微熟悉，但认为它是一种一次性方法，用于计算整个样本的单个方差，其额外优势是可以将问题分解为并行运行的部分。

Chan 等人给出了一种方法，可以在给定部件统计信息的情况下计算部件连接的统计信息。高级的想法是维护一个部分的集合（实际上只是它们的统计数据），这样任何窗口都是 O(log k) 部分的串联。一种方法是使用平衡二叉树，但正如 Rex 指出的那样，这太过分了，我们可以只维护大小为 2 的幂的对齐部分的统计信息（例如，[0, 1）、[1, 2）、[0 , 2), [2, 3), [3, 4), [2, 4), [0, 4), 等等)

【参考方案6】：

我知道这个问题很老，但如果其他人对此感兴趣，请遵循 python 代码。它的灵感来自johndcook 博客文章、@Joachim 的、@DanS 的代码和@Jaime cmets。下面的代码仍然对小数据窗口大小给出了小的不精确性。享受吧。

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))

【讨论】：

感谢python中的idea综合。我不喜欢在输入 else 块的情况下窗口的大小如何变为WIN_SIZE - 1。因此，如果 WIN_SIZE 在调用 push 并且我们追加时为 10，由于使用了 deque 构造函数选项，它仍然是 10，然后在 else 块中 popleft 将大小进一步减小到 9。所以也许 maxlen=WIN_SIZE + 1 ?或者不使用maxlen 选项。此外，可以删除n 变量并使用len(self.windows)。

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

【参考方案7】：

我期待在这方面被证明是错误的，但我认为这不能“快速”完成。也就是说，计算的很大一部分是在窗口上跟踪 EV，这很容易完成。

我会带着问题离开：你确定你需要一个窗口函数吗？除非您使用非常大的窗口，否则最好只使用众所周知的预定义算法。

【参考方案8】：

我猜想跟踪您的 20 个样本 Sum(X^2 from 1..20) 和 Sum(X from 1..20)，然后在每次迭代中连续重新计算这两个总和是否不够高效？可以重新计算新的方差，而无需每次将所有样本相加、平方等。

如：

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

【讨论】：

我相信这个解决方案容易受到我原始帖子 (johndcook.com/standard_deviation.html) 中链接中提到的稳定性问题的影响。特别是，当输入值和大且它们的差值小于时，结果实际上可能是负数。我无法控制输入，所以我宁愿避免这种方法。

哦，我明白了。关于输入你有什么可以说的吗？有可能的使用？您可以在（64 位浮点数、任意精度算术等）上抛出更多位，这是一个问题吗？如果您以有效数字胜过输入，舍入错误就会消失，不是吗？

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

@Jason S：滚动方差就是这样。从 100 万到零的过渡可能会发生很多事情，但这就是野兽的本性。那，当无论如何发生变化时，1000 ~ 100 万个值中的前 980 个都不存在。我的评论表明，如果您的计算中有足够的有效数字，这些都无关紧要。

输入可以是任何东西。价值量级当然可以达到数万亿，虽然原始数据只能精确到小数点后几位，但用户可以在计算方差之前转换他们的数据（例如除以任何标量）。

【参考方案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 个平方数相加非常快。如果你做得更多——比如说几百个——更有效的方法显然会更好。但我不确定蛮力不是这里的方法。

【讨论】：

"使用您的标准 E(x^2)-E(x)^2 公式" 不，不要；它甚至都不是很稳定。采用一种更好的算法。

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

我很担心，因为这是个问题。去阅读***的文章，如果你仍然不相信，试着计算 N(1, 1e-10) 的 20 个独立同分布样本的方差。

对于任何具有合理单位和来源的现实数据集，我还没有看到这实际上是一个问题，但足够公平，如果这就是 OP 想要的......

【参考方案10】：

对于只有 20 个值，调整暴露的方法 here 是微不足道的（不过我没说快）。

您可以简单地选择包含 20 个 RunningStat 类的数组。

流的前 20 个元素有些特殊，但是一旦完成，就简单多了：

当一个新元素到达时，清除当前的RunningStat 实例，将该元素添加到所有 20 个实例中，并递增标识新的“完整”RunningStat 实例的“计数器”（模 20） 在任何给定时刻，您都可以查阅当前的“完整”实例以获取正在运行的变体。

您显然会注意到，这种方法并不是真正可扩展的......

您还可以注意到，我们保留的数字有些冗余（如果您选择 RunningStat 全班）。一个明显的改进是直接保留Mk 和Sk 的20 个last。

我想不出使用这种特定算法的更好公式，恐怕它的递归公式有点束缚我们的手。

【参考方案11】：

这只是对 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 是您希望删除的窗口中最旧的样本。