小智 23
我也遇到过这个问题.在计算运行累积方差方面有一些很好的帖子,例如John Cooke的准确计算运行方差帖子和数字探索的帖子,用于计算样本和人口方差的Python代码,协方差和相关系数.只是找不到任何适合滚动窗口的东西.
Subluminal Messages发布的运行标准偏差对于使滚动窗口公式起作用至关重要.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 - 今天是 - n)/ n
- x =您的时间序列中的值
- n =用于滚动窗口的周期.
从那里你可以计算滚动人口差异:
今日人口变量=(PSA今天*n - n*SMA今天*SMA今天)/ n
或滚动样本差异:
今天的样本Var =(PSA今天*n - n*SMA今天*SMA今天)/(n - 1)
几年前我在博客文章" 运行方差"中介绍了这个主题以及示例Python代码.
希望这可以帮助.
请注意:我提供了Latex(图像)中所有博客文章和数学公式的链接.但是,由于我的声誉很低(<10); 我只限于2个超链接,绝对没有图像.为此表示歉意.希望这不会带走内容.
Dan*_*anS 20
我一直在处理同样的问题.
均值很容易迭代计算,但您需要将值的完整历史记录保存在循环缓冲区中.
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的帖子):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);
}
}
实际上，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 ，它还与精确的两遍方差进行比较。
但是，在倾斜的数据集上，行为可能会有所不同。该数据集显然是均匀分布的；但是我也尝试了一个排序数组，它起作用了。
更新：我们发表了一篇论文，详细介绍了（协方差）的不同加权方案：
Schubert，Erich和Michael Gertz。“ （协）方差的数值稳定并行计算。 ”第30届科学与统计数据库管理国际会议论文集。ACM，2018年。（获得SSDBM最佳论文奖。）
这也讨论了如何使用加权来并行化计算，例如在AVX，GPU或群集上。
小智 5
这是一个有O(log k)时间更新的分而治之的方法,其中k是样本数.它应该是相对稳定的,因为成对求和和FFT是稳定的,但它有点复杂,常数不是很大.
假设我们有一个序列A长度的m均值E(A)和方差V(A),以及序列B长度的n均值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风格.
我知道这个问题很老，但如果其他人对此感兴趣，请遵循 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))
我期待着在这一点上被证明是错误的，但我不认为这可以“很快”完成。也就是说，计算的很大一部分是跟踪窗口上的 EV，这很容易完成。
我将留下一个问题：您确定需要窗口函数吗？除非您使用非常大的窗口，否则最好使用众所周知的预定义算法。