flint>flint blog>2013年> 8月>28日>相関係数のインクリメンタル計算

相関係数のインクリメンタル計算

先月と今月は、ソフトウェア開発 (この記事で触れた案件の続き) のお仕事を在宅で。 今回作成したプログラムは、最初に先方から説明を受けた時点では処理の内容も単純で、それほど難度の高くないものだと思っていましたが、実は意外なところに罠が潜んでいました。

サンプル (標本) が逐次追加されていく集合の間の相関係数を計算し、その結果をリアルタイム表示する、というのが本件で開発するソフトウェアに要求される主要な機能だったのですが、これがなかなかのクセモノ。 扱うサンプルの数は数千~数万個のオーダなのですが、その集合に対して、リアルタイムで4,096通りの組み合わせで相関係数を算出する必要があります。 処理対象となるサンプルの数が増えるにつれ、計算処理に要する時間は急速に増加。 表示の更新間隔はどんどん間延びしていき、傍から見ているぶんには、フリーズしている状態と区別が付きません。 (一応、計算処理はスレッド化しているので、任意のタイミングで中断させることはできますが、根本的な解決にはなっていない。)

そこで、前回の反復 (イテレーション) における計算結果を記憶しておき、これに次の反復で新しく入ってきたサンプル値による影響を加えていくことで計算量を削減する、という手法を導入することにしました。 そのために、ノートに数ページに渡って数式の変形過程を書き連ねてるハメになった (主に計算ミスが多かったため) わけですが、その苦労の甲斐あって、色々と応用が利きそうなアルゴリズムを習得することができたので、記事にまとめておこうと思った次第です。

平均は簡単

相関係数を求めるには、まず処理の対象となる集合の平均 (相加平均) を計算する必要があります。 ご存知の通り、集合 Xn = { x0, x1, ..., xn-1 } の平均は次の式で与えられます。

平均 (n)

これに新しい要素 xn を加えた集合 Xn+1 の平均は、先に求めた Xn の平均を用いて次のように求められます。

平均 (n + 1)

以前の平均値に n/(n + 1)、新しい要素値に 1/(n + 1) の重みを付けて平均を計算すれば新しい平均値を得ることができる、という寸法。 直感的にも理解しやすい式ですね。

分散を求める

相関係数の計算式を見てみると、全体の分母に当たる部分が対象となる二つの集合 Xn, Yn それぞれの標準偏差の積になっていることに気が付くでしょう。 標準偏差を以前の反復の結果から直接計算するのは大変なので、その代わりとして、標準偏差の平方 (自乗) にあたる分散を求めます。 分散の定義は以下の通り。

分散 (n)

次の反復における分散の値は、これまでに計算された値を用いて次のように求めることができます。

分散 (n + 1)

平均の場合と同様に、前の反復における分散に n/(n + 1)、新しく追加されたサンプルと新旧ふたつの平均の差の積に 1/(n + 1) の重みを付けて、平均を計算してやればよいわけです。 個人的には、とてもエレガントな式だと感じられるのですが、如何でしょうか。

共分散を求める

今度は相関係数の全体の分子に当たる部分に着目してみると、対象となる集合 Xn と、Yn共分散になっていることが分かります。

共分散 (n)

実は前節で示した分散の計算式は、共分散のそれから求めたもの。(分散は共分散の特殊な形 - 等しい二つの集合に対する共分散であるため。) 細かい説明は省きますが、前の反復で求めた共分散の値から次の共分散の値を求める式は以下の様になります。

共分散 (n + 1)

二つの式のどちらを使っても構いません。 一応、式の中の相違部分が等値であることの確認もしておきました。

「答合わせ」をしよう

これで、平均, 標準偏差 (実際にはその平方である分散) および共分散をインクリメンタルに計算できるようになったので、これらを組み合わせて最終的な目標である相関係数のインクリメンタルな計算を行うことも可能になりました。 テストを行っても、10,000サンプル程度では普通に計算した結果とのズレは検出されません。 (検証には倍精度浮動小数点数 (double 型) を使用。) しかし、まだ消えない不安。 一般的に良く知られた、より簡潔な別の計算方法があるかも知れない。 だとしたら、車輪の再発明をした上に、先発の手法よりも劣化した計算方法をドヤ顔で採用した自分はとんだおマヌケさん、ということになってしまいます。

そんなときにはウェブ検索。 日本語ではそれらしい情報は見当たりません。 そこでキーワードを英語 ("incremental calculation", "variance" など) に切り替えてみると、まず論文が出てきました。 著者はケンブリッジ大学の方のようです。

Incremental calculation of weighted mean and variance
http://nfs-uxsup.csx.cam.ac.uk/~fanf2/hermes/doc/antiforgery/stats.pdf

式の形は異っていますが、読み替えと変形によって本質的に同じアルゴリズムであることが確認できました。 ですが、この論文には共分散の計算アルゴリズムについての言及がありません。 そこでもう少し検索を続けたたところ......なんと、Wikipedia にズバリそのもののページがあるではありませんか。

Algorithms for calculating variance - Wikipedia, the free encyclopedia
http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

共分散を計算するアルゴリズムについても、バッチリ記載されています。 いいんだ......自分で計算してみることに意義があったんだい。(涙)

余談ですが、英語版の Wikipedia は、数学系の記事がサブカル系の記事ばかり載っている日本語版と比較してとても充実しています。 以下の記事などは、眺めているだけでワクワクしますよね。 そして、あっという間に過ぎる時間。 中毒性が高いので要注意です。

List_of_fractals_by_Hausdorff_dimension - Wikipedia, the free encyclopedia
http://en.wikipedia.org/wiki/List_of_fractals_by_Hausdorff_dimension
成田 (添え字が 0 ~ n - 1 なのはプログラマの習性。)
このエントリーをはてなブックマークに追加

コメント

投稿者
URI
メールアドレス
表題
本文