データが逐次追加されていく際に、追加されるたびにその時点での「分散」や「標準偏差」を計算したい場合がある。その時点での全てのデータから毎度計算しても良いが、やはり計算量が馬鹿らしい。そこで欲しくなるのが、これらの量をオンライン（ストリーム処理）で計算できるアルゴリズムだ。
安心してください、ありますよ。「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\)である。
この定義通りに分散を計算する場合、以下の２ステップを辿る。
- データ全体の平均\(\bar{x}\)を計算。
- 各データ\(x_i\)と平均\(\bar{x}\)の差分の二乗を計算。
■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投げて頂ければ嬉しいです。
■参考
