2020年11月15日日曜日

分散や標準偏差のオンライン計算 → Welfordアルゴリズム

データが逐次追加されていく際に、追加されるたびにその時点での「分散」や「標準偏差」を計算したい場合がある。その時点での全てのデータから毎度計算しても良いが、やはり計算量が馬鹿らしい。そこで欲しくなるのが、これらの量をオンライン(ストリーム処理)で計算できるアルゴリズムだ。

安心してください、ありますよ。「Welford アルゴリズム」というものです。

ここではそのWelfordアルゴリズムを紹介したい。

※以下、不偏分散を考えるが、標本分散の場合でも全く同じ考えが出来る。

■分散と標準偏差

データサンプルが\(x_i, \ldots, x_N\)で与えられる場合を考えられる。この時、データの不偏分散\(s^2\)の定義は

\[s^2 = \frac{\sum_{i=1}^N (x_i – \bar{x})^2}{N-1}\]

として与えられる。そして標準偏差\(s\)はその平方根\(s = \sqrt{s^2}\)だ。

ここで\(\bar{x}\)はサンプルの平均、すなわち\(\bar{x} = \frac{1}{N}\sum_{i=1}^N x_i\)である。

この定義通りに分散を計算する場合、以下の2ステップを辿る。

  1. データ全体の平均\(\bar{x}\)を計算。
  2. 各データ\(x_i\)と平均\(\bar{x}\)の差分の二乗を計算。
このような計算は明らかに無駄が多い。まずデータ全体を舐める計算を2回も繰り返す必要がある。また、データ全体に対する計算を行うためにすべてのデータを保持しないといけないので今回の主題であるオンラインの計算をするためにはこのままではダメだ。何か工夫がいる。

■Welfordアルゴリズム

そこでWelfordアルゴリズムでは、サンプルが\(N\)個の時と\(N-1\)個の時の分散の差に着目して以下の計算を行う。

\begin{align} &(N-1)s_N^2 – (N-2)s_{N-1}^2 \\ &= \sum_{i=1}^N (x_i-\bar{x}_N)^2-\sum_{i=1}^{N-1} (x_i-\bar{x}_{N-1})^2 \\ &= (x_N-\bar{x}_N)^2 + \sum_{i=1}^{N-1}\left((x_i-\bar{x}_N)^2-(x_i-\bar{x}_{N-1})^2\right) \\ &= (x_N-\bar{x}_N)^2 + \sum_{i=1}^{N-1}(x_i-\bar{x}_N + x_i-\bar{x}_{N-1})(\bar{x}_{N-1} – \bar{x}_{N}) \\ &= (x_N-\bar{x}_N)^2 + (\bar{x}_N – x_N)(\bar{x}_{N-1} – \bar{x}_{N}) \\ &= (x_N-\bar{x}_N)(x_N-\bar{x}_N – \bar{x}_{N-1} + \bar{x}_{N}) \\ &= (x_N-\bar{x}_N)(x_N – \bar{x}_{N-1}) \\ \end{align}

この結果から、下式のようにデータが\(N\)個の時の分散を\(N-1\)個の時の分散から求められることがわかる。

\[ s_N^2 = \frac{N-2}{N-1} s_{N-1}^2 + \frac{1}{N-1} (x_N-\bar{x}_N)(x_N – \bar{x}_{N-1}) \]

この結果をもとに分散をオンラインで計算するアルゴリズムに落とすと、下の疑似コードのようになる。(forで各データを逐次回している部分がそれだ)

驚くほど簡単なアルゴリズムだ。

■ライブラリ

簡単なので必要に応じて自分で実装すれば良いが、Python用のライブラリを作ってPypiに登録してみたのでよければそれを使ってみてください。改良点があればGithubリポジトリにIssue投げて頂ければ嬉しいです。


■参考

0 件のコメント:

コメントを投稿