74

ローリング分散 (たとえば、20 期間のローリング ウィンドウにわたる分散) を計算するための効率的で数値的に安定したアルゴリズムを見つけようとしています。一連の数値の連続分散を効率的に計算するWelford アルゴリズムは知っていますが (必要なパスは 1 つだけです)、ローリング ウィンドウに適用できるかどうかはわかりません。また、この記事の冒頭でJohn D. Cookが説明した精度の問題を回避するためのソリューションも希望しています。どの言語でのソリューションでも問題ありません。

4

11 に答える 11

28

私もこの問題に遭遇しました。実行中の累積分散の計算に関する素晴らしい投稿がいくつかあります。たとえば、John Cooke のAccurately Computing running Variationの投稿や、デジタル探索、サンプルと母集団の分散を計算するための Python コード、共分散と相関係数の投稿などです。ローリングウィンドウに適応したものを見つけることができませんでした.

Subluminal Messages によるRunning Standard Deviationsの投稿は、ローリング ウィンドウ式を機能させる上で重要でした。Jim は、平均値の差の 2 乗の和を使用するという Welford のアプローチに対して、値の差の 2 乗の累乗和を取ります。式は次のとおりです。

今日の 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 = ローリング ウィンドウに使用される期間。

そこから、Rolling Population Variance を計算できます。

今日の人口変数 = (今日の PSA * n - n * 今日の SMA * 今日の SMA) / n

またはローリング サンプル分散:

今日のサンプル変数 = (今日の 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 に答える
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);
}
}

正確な2パスアルゴリズムと比較して、約14桁の精度が得られます。これは、ダブルスから期待できる程度です。ウェルフォードは、余分な分割のために計算コストがかかることに注意してください。正確な2パスアルゴリズムの約2倍の時間がかかります。ウィンドウサイズが小さい場合は、実際に平均を再計算してから、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、またはクラスタなどで、重み付けを使用して計算を並列化する方法についても説明します。

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

ここで、要素を赤黒木に詰め込みます。各ノードは、そのノードをルートとするサブツリーの平均と分散で装飾されます。右側に挿入します。左側を削除します。(端にしかアクセスしていないので、スプレイ ツリーは 償却される可能性がありますが、アプリケーションにとっては償却が問題になると思います) 。O(1)k

于 2011-02-28T21:53:47.307 に答える
1

別のO(log k)解決策は次のとおりです。元のシーケンスの2乗を見つけ、次にペアを合計し、次に4倍などを見つけます（これらすべてを効率的に見つけるには、少しのバッファーが必要です）。次に、必要な値を合計します。あなたの答えを得るために。例えば：

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

これが間違っていることが証明されることを楽しみにしていますが、これが「すぐに」できるとは思いません。とはいえ、計算の大部分は、簡単に実行できるウィンドウ上の EV を追跡することです。

私は質問を残します: 本当にウィンドウ関数が必要ですか? 非常に大きなウィンドウで作業している場合を除き、よく知られている事前定義されたアルゴリズムを使用することをお勧めします。

于 2011-02-28T20:57:42.793 に答える
1

値が 20 個しかない場合は、ここで公開されているメソッドを適応させるのは簡単です(ただし、高速とは言いませんでした)。

RunningStatこれらのクラスの 20 個の配列を簡単に取得できます。

ストリームの最初の 20 要素は多少特殊ですが、これが完了すると、はるかに単純になります。

  • 新しい要素が到着すると、現在のインスタンスをクリアし、その要素を 20 個のインスタンスすべてに追加し、新しい「フル」インスタンスRunningStatを識別する「カウンター」(モジュロ 20) をインクリメントします。RunningStat
  • いつでも、現在の「完全な」インスタンスを参照して、実行中のバリアントを取得できます。

このアプローチは実際にはスケーラブルではないことに気付くでしょう...

また、保持している数値には冗長性があることにも注意してください (RunningStat完全なクラスを使用する場合)。明らかな改善は、20 のラストMkをSk直接保持することです。

この特定のアルゴリズムを使用したより良い式は思いつきません。残念ながら、その再帰的な定式化は私たちの手を縛っています。

于 2011-03-01T09:01:37.143 に答える