滚动方差algorithm

我试图find一个有效的，数值稳定的algorithm来计算滚动方差（例如，在20周期滚动窗口上的方差）。 我知道Welfordalgorithm可以有效地计算一串数字的运行方差（它只需要一次通过），但我不确定这是否可以适应滚动窗口。 我也想要避免在本文顶部讨论的准确性问题的解决scheme。 任何语言的解决scheme都很好。

Answer 1

我也遇到了这个问题。 在计算运行累积方差（如John Cooke的精确计算运行方差post和数字探索的post， 用于计算样本和总体方差的Python代码，协方差和相关系数）方面，有一些很棒的post。 就是找不到适合滚动窗口的东西。

“Subluminal Messages”中的“ 运行标准偏差”对于使滚动窗口公式起作用至关重要。 吉姆采用数值的平方差的功率和与利用平均值的平方差的总和的韦尔福德的方法。 公式如下：

今日PSA = PSA（昨天）+（（（今天×今天）×昨天））/ n

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

但是，要将Power Sum Average公式转换为窗口types，您需要将公式调整为以下值：

今日PSA =昨日PSA +（（（x今日* x今日） – （x昨日* x昨日）/ n

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

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

今天的SMA =昨天的SMA +（（x今天 – x今天 – n）/ n

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

从那里你可以计算滚动人口差异：

今日人口Var =（今天的PSA * n – 今天的SMA *今天的SMA）/ n

或滚动样本差异：

今天的样本Var =（今天的PSA *今天的* SMA *今天的SMA）/（n – 1）

我已经在几年前的一篇博客文章中介绍了这个主题以及示例Python代码， 运行变化 。

希望这可以帮助。

请注意：我为这个答案提供了所有乳胶博客文章和math公式的链接（图片）。 但是，由于我的低声誉（<10）; 我仅限于2个超链接，绝对没有图像。 为此事道歉。 希望这不会带走内容。

Answer 2

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

平均值迭代计算很简单，但是您需要将值的完整历史logging保存在循环缓冲区中。

 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的algorithm，它适用于我所testing过的所有值。

 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;

Answer 3

如果你喜欢代码的话（主要基于DanS的post）： 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); } } 

Answer 4

这是一个分而治之的方法，它具有O(log k)时间更新，其中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))) 

现在，将这些元素填充到红黑树中，每个节点都装有以该节点为根的子树的均值和方差。 插入右侧; 删除左边。 （因为我们只是访问结束，一个splay树可能是O(1)分期付款，但我猜分期付款是您的应用程序的一个问题。）如果在编译时知道k ，你可能可以展开内部循环FFTW风格。

Answer 5

实际上Welfordsalgorithm可以很容易地将AFAICT用于计算加权方差。 通过将权重设置为-1，您应该能够有效地抵消元素。 我没有检查math是否允许负面的权重，但在第一眼看来应该！

我用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); } } 

与精确的双通道algorithm相比，我可以获得大约14位的精度。 这与双打的预期差不多。 请注意，由于额外的划分，Welford 确实会带来一些计算成本 – 这大约是确切的两遍algorithm的两倍。 如果窗口的大小很小，实际上重新计算平均值，然后在每次传递方差时可能更为明智。

我已经将这个实验作为unit testing添加到ELKI中，您可以在这里看到完整的源代码： http : //elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki /math/TestSlidingVariance.java，它也比较准确的两遍方差。

但是，对于偏斜的数据集，其行为可能会有所不同。 这个数据集明显是均匀分布的; 但我也尝试了一个sorting的数组，它的工作。

Answer 6

我期待在这方面被certificate是错误的，但我不认为这可以“迅速”完成。 也就是说，计算的很大一部分是跟踪可以轻松完成的窗口上的电动车。

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

Answer 7

我想跟踪你的20个样本Sum（X ^ 2从1..20）和Sum（X从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 

Answer 8

这里有另一个O(log k)解：find原始序列的平方，然后求和，然后四倍等。（你需要一点缓冲才能find所有这些有效的。）然后加起来那些你需要得到你的答案的价值。 例如：

Answer 9

Answer 10

Answer 11

Answer 12

Answer 13

Answer 14

Answer 15

Answer 16

| // 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个平方的数字是非常快的。 如果你做得更多 – 比方说几百个 – 更有效率的方法显然会更好。 但我不确定蛮力是不是要走到这里的路。

Answer 17

对于只有20个值，适应这里暴露的方法是微不足道的（虽然我没有说快）。

您可以简单地选取这20个RunningStat类的数组。

stream的前20个元素有些特别，但是一旦完成，就更简单了：

• 当一个新的元素到达时，清除当前的RunningStat实例，将这个元素添加到所有的20个实例中，然后递增“counter”（模20），标识新的“完整” RunningStat实例
• 在任何时候，你可以参考当前的“完整”实例来获得你的运行变体。

显然你会注意到这种方法并不是真正可扩展的。

你也可以注意到，我们保留的数字有一些减less（如果你使用RunningStat全class）。 一个明显的改进就是直接保留20个Mk和Sk 。

我想不出使用这个特定algorithm的更好的公式，恐怕它的recursion公式有点牵扯我们的手。

Answer 18

我知道这个问题是旧的，但如果有人在这里感兴趣，这里是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))