롤링 분산 알고리즘

롤링 분산 알고리즘

---

롤링 분산 (예 : 20주기 롤링 윈도우에 대한 분산)을 계산하기 위해 효율적이고 수치 적으로 안정적인 알고리즘을 찾으려고합니다. 숫자 스트림에 대한 실행 분산을 효율적으로 계산 하는 Welford 알고리즘 을 알고 있지만 (한 번의 패스 만 필요함) 이것이 롤링 창에 적용될 수 있는지 확실하지 않습니다. 또한 John D. Cook 이이 기사 의 상단에서 논의한 정확도 문제를 피할 수있는 솔루션을 원합니다 . 모든 언어의 솔루션이 좋습니다.

---

이 문제도 겪었습니다. John Cooke의 Accurately compute running variance post 및 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 = 지금까지 분석 한 값의 수.

또한 Rolling Simple Moving Average 공식이 필요합니다.

오늘 SMA = 어제 SMA + ((오늘 x-x 오늘-n) / n

• x = 시계열 값
• n = 롤링 윈도우에 사용 된 기간.

여기에서 Rolling Population Variance를 계산할 수 있습니다.

오늘의 인구 변동 = (오늘 PSA * n-n * 오늘 SMA * 오늘 SMA) / n

또는 롤링 표본 분산 :

오늘 샘플 변수 = (오늘 PSA * n-n * 오늘 SMA * 오늘 SMA) / (n-1)

몇 년 전, Running Variance 라는 블로그 게시물에서 샘플 Python 코드와 함께이 주제를 다루었습니다 .

도움이 되었기를 바랍니다.

참고 :이 답변에 대한 모든 블로그 게시물 및 Latex (이미지)의 수학 공식에 대한 링크를 제공했습니다. 하지만 저의 평판이 낮기 때문에 (<10); 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;

---

단어보다 코드를 선호하는 경우 (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);
    }
}

---

Actually Welfords algorithm can AFAICT easily be adapted to compute weighted Variance. And by setting weights to -1, you should be able to effectively cancel out elements. I havn't checked the math whether it allows negative weights though, but at a first look it should!

I did perform a small experiment using 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);
}
}

정확한 2- 패스 알고리즘에 비해 약 14 자리의 정밀도를 얻습니다. 이것은 복식에서 기대할 수있는 정도입니다. Welford 는 추가 분할로 인해 계산 비용이 발생합니다. 정확한 2- 패스 알고리즘보다 약 두 배가 걸립니다. 창 크기가 작 으면 실제로 평균을 다시 계산 한 다음 매번 분산을 전달하는 것이 훨씬 더 합리적 일 수 있습니다 .

이 실험을 ELKI에 단위 테스트로 추가했습니다. 여기에서 전체 소스를 볼 수 있습니다. http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki /math/TestSlidingVariance.java 또한 정확한 2 단계 분산과 비교합니다.

그러나 치우친 데이터 세트에서는 동작이 다를 수 있습니다. 이 데이터 세트는 분명히 균일하게 분포되어 있습니다. 그러나 나는 또한 정렬 된 배열을 시도했으며 작동했습니다.

업데이트 : (공) 분산에 대한 다양한 가중치 체계에 대한 세부 정보가 포함 된 논문을 게시했습니다.

슈베르트, 에리히, 마이클 거츠. " (공) 분산의 수치 적으로 안정된 병렬 계산. "과학 및 통계 데이터베이스 관리에 관한 제 30 차 국제 컨퍼런스의 진행. ACM, 2018. (SSDBM 최우수 논문상 수상)

또한 가중치를 사용하여 계산을 병렬화하는 방법 (예 : AVX, GPU 또는 클러스터)에 대해서도 설명합니다.

---

다음 은 샘플 수인 O(log k)시간 업데이트 가있는 분할 및 정복 접근 방식입니다 k. pairwise summation과 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)))

Now, stuff the elements in a red-black tree, where each node is decorated with mean and variance of the subtree rooted at that node. Insert on the right; delete on the left. (Since we're only accessing the ends, a splay tree might be O(1) amortized, but I'm guessing amortized is a problem for your application.) If k is known at compile-time, you could probably unroll the inner loop FFTW-style.

---

I know this question is old, but in case someone else is interested here follows the python code. It is inspired by johndcook blog post, @Joachim's, @DanS's code and @Jaime comments. The code below still gives small imprecisions for small data windows sizes. Enjoy.

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

---

I look forward to be proven wrong on this but I don't think this can be done "quickly." That said, a large part of the calculation is keeping track of the EV over the window which can be done easily.

I'll leave with the question: are you sure you need a windowed function? Unless you are working with very large windows it is probably better to just use a well known predefined algorithm.

---

I guess keeping track of your 20 samples, Sum(X^2 from 1..20), and Sum(X from 1..20) and then successively recomputing the two sums at each iteration isn't efficient enough? It's possible to recompute the new variance without adding up, squaring, etc., all of the samples each time.

As in:

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

---

Here's another O(log k) solution: find squares the original sequence, then sum pairs, then quadruples, etc.. (You'll need a bit of a buffer to be able to find all of these efficiently.) Then add up those values that you need to to get your answer. For example:

|||||||||||||||||||||||||  // 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.

Now you use your standard E(x^2)-E(x)^2 formula and you're done. (Not if you need good stability for small sets of numbers; this was assuming that it was only accumulation of rolling error that was causing issues.)

That said, summing 20 squared numbers is very fast these days on most architectures. If you were doing more--say, a couple hundred--a more efficient method would clearly be better. But I'm not sure that brute force isn't the way to go here.

---

For only 20 values, it's trivial to adapt the method exposed here (I didn't say fast, though).

You can simply pick up an array of 20 of these RunningStat classes.

The first 20 elements of the stream are somewhat special, however once this is done, it's much more simple:

• when a new element arrives, clear the current RunningStat instance, add the element to all 20 instances, and increment the "counter" (modulo 20) which identifies the new "full" RunningStat instance
• at any given moment, you can consult the current "full" instance to get your running variant.

You will obviously note that this approach isn't really scalable...

You can also note that there is some redudancy in the numbers we keep (if you go with the RunningStat full class). An obvious improvement would be to keep the 20 lasts Mk and Sk directly.

I cannot think of a better formula using this particular algorithm, I am afraid that its recursive formulation somewhat ties our hands.

---

This is just a minor addition to the excellent answer provided by DanS. The following equations are for removing the oldest sample from the window and updating the mean and variance. This is useful, for example, if you want to take smaller windows near the right edge of your input data stream (i.e. just remove the oldest window sample without adding a new sample).

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)

Here, x_old is the oldest sample in the window you wish to remove.

참고URL : https://stackoverflow.com/questions/5147378/rolling-variance-algorithm