74

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

4

11 回答 11

28

我也遇到过这个问题。在计算运行累积方差方面有一些很棒的帖子,例如 John Cooke 的准确计算运行方差帖子和来自数字探索的帖子,用于计算样本和总体方差、协方差和相关系数的 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 = 用于滚动窗口的时间段。

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

今日人口变异 = (今日 PSA * n - n * 今日 SMA * 今日 SMA) / n

或滚动样本方差:

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

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

希望这可以帮助。

请注意:我为此答案提供了 Latex(图像)中所有博客文章和数学公式的链接。但是,由于我的声誉低(< 10);我仅限于 2 个超链接,绝对没有图像。为此表示歉意。希望这不会影响内容。

于 2012-07-29T18:10:53.430 回答
27

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

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

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;

于 2011-07-12T12:33:27.267 回答
9

如果您更喜欢代码而不是文字(主要基于 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);
    }
}
于 2014-04-24T18:32:40.050 回答
6

实际上 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 或在集群上。

于 2013-01-05T13:47:57.507 回答
5

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

假设我们有一个带有均值和方差的A长度序列,以及一个带有均值和方差的长度序列。让是 和 的串联。我们有mE(A)V(A)BnE(B)V(B)CAB

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 样式。

于 2011-02-28T21:53:47.307 回答
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))
于 2017-08-29T22:26:04.140 回答
1

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

于 2011-02-28T22:26:00.743 回答
1

我猜想跟踪你的 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
于 2011-02-28T21:07:20.403 回答
1

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

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

于 2011-02-28T20:57:42.793 回答
1

对于只有 20 个值,调整这里公开的方法是微不足道的(不过我没说快)。

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

流的前 20 个元素有些特殊,但是一旦完成,它就会简单得多:

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

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

您还可以注意到,我们保留的数字有一些冗余(如果您参加RunningStat全班课程)。Mk一个明显的改进是直接保留 20 个鞋楦Sk

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

于 2011-03-01T09:01:37.143 回答
0

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

于 2018-04-09T15:38:01.660 回答