分散を計算する
一般に、手元にあるデータ（標本） $ x_1, x_2, x_3 ...x_n $の$ N $個データがあり、その平均を $ \mu $とした時、その分散 $ { s }^{ 2 } $は、以下のように得られます。
{s}^{ 2 }=\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ { \left( \mu -x_{ i } \right) }^{ 2 } }
これをそのままプログラミングで平凡に計算しようとした場合、例えば以下のようになります。
#include <stdio.h>
#include <random>
int main()
{
std::vector<double> data;
std::mt19937 engine;
std::uniform_real_distribution<> random(-1.0, 1.0);
for (int i = 0; i < 100000; ++i) {
data.push_back(4.0 + random(engine) + random(engine) + random(engine) + random(engine));
}
double sum = 0.0;
for (int i = 0; i < data.size(); ++i) {
sum += data[i];
}
double average = sum / data.size();
double variance_sum = 0.0;
for (int i = 0; i < data.size(); ++i) {
variance_sum += (average - data[i]) * (average - data[i]);
}
double variance = variance_sum / data.size();
printf("avg = %.4f\n", average);
printf("variance = %.4f\n", variance);
printf("sd = %.4f\n", sqrt(variance));
return 0;
}
インクリメンタル？
これでめでたしめでたし、となれば別にこれで良いです。
でも場合によっては、データが計算時にすべて揃っていない場合や、どんどんデータが後から増える場合、dataを２回走査することを嫌いたい場合、などなどこの方法が最適では無い場合もあります。
来たデータからどんどん分散を計算できて、なおかつ新しいデータが入ってきても問題なく対応できるような方法は無いでしょうか？
分散 = 2乗の平均 - 平均の2乗
そこで分散が、２乗の平均から平均の２乗を引いたものに等しくなるという性質を使って、式を変形してみましょう。
参考
http://www.geisya.or.jp/~mwm48961/kou3/prob_variance1.htm
{ s }^{ 2 }=\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ { \left( \mu -x_{ i } \right) }^{ 2 } } \\ =\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ { x_{ i } }^{ 2 } } -{ \left( \frac { 1 }{ N } \sum _{ i=1 }^{ N }{ x_{ i } } \right) }^{ 2 }\\ =\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ { x_{ i } }^{ 2 } } -\frac { 1 }{ N } \frac { 1 }{ N } { \left( \sum _{ i=1 }^{ N }{ x_{ i } } \right) }^{ 2 }\\ =\frac { 1 }{ N } \left\{ \sum _{ i=1 }^{ N }{ { x_{ i } }^{ 2 } } -\frac { 1 }{ N } { \left( \sum _{ i=1 }^{ N }{ x_{ i } } \right) }^{ 2 } \right\}
試しにプログラムに落としてみます。
#include <stdio.h>
#include <random>
int main()
{
std::vector<double> data;
std::mt19937 engine;
std::uniform_real_distribution<> random(-1.0, 1.0);
for (int i = 0; i < 100000; ++i) {
data.push_back(4.0 + random(engine) + random(engine) + random(engine) + random(engine));
}
int n = 0;
double sum_of_squared = 0.0;
double sum = 0.0;
for (double value : data) {
sum_of_squared += value * value;
sum += value;
n++;
}
double variance = (sum_of_squared - sum * sum / n) / n;
double average = sum / n;
printf("avg = %.4f\n", average);
printf("variance = %.4f\n", variance);
printf("sd = %.4f\n", sqrt(variance));
return 0;
}
これももちろん最初の計算と結果が一致します。（ただし当然ではありますが、コンピュータによる浮動小数点数演算は誤差が生じるので、完全には一致しないことは注意が必要です。）
さて、最初の方法と比較してみると、大きな違いがあることに気づきます。
プログラムに落とし込む時にループが１つで良いことと、$\mu$が分解されたことです。
これを例えばクラスで表現すると、より本質が浮き彫りになります。
#include <stdio.h>
#include <random>
class IncrementalStatatics {
public:
void addSample(double value) {
_sum_of_squared += value * value;
_sum += value;
_n++;
}
double variance() const {
return (_sum_of_squared - _sum * _sum / _n) / _n;
}
double avarage() const {
return _sum / _n;
}
private:
double _sum_of_squared = 0.0;
double _sum = 0.0;
int _n = 0;
};
int main()
{
std::vector<double> data;
std::mt19937 engine;
std::uniform_real_distribution<> random(-1.0, 1.0);
for (int i = 0; i < 100000; ++i) {
data.push_back(4.0 + random(engine) + random(engine) + random(engine) + random(engine));
}
IncrementalStatatics statatics;
for (double value : data) {
statatics.addSample(value);
}
printf("avg = %.4f\n", statatics.avarage());
printf("variance = %.4f\n", statatics.variance());
printf("sd = %.4f\n", sqrt(statatics.variance()));
return 0;
}
これでたとえループの途中であっても、いつでもその時点での統計量を取り出すことが可能になりました。
母集団の分散について
手元にある分散を考えた時は、上記で問題ないのですが、母集団の分散を推定するという立場からすると、標本の分散の期待値というのは母集団の分散と若干の差があります。
なので、それをunbiasedにするためには、
{ s }^{ 2 }=\frac { 1 }{ N-1 } \left\{ \sum _{ i=1 }^{ N }{ { x_{ i } }^{ 2 } } -\frac { 1 }{ N } { \left( \sum _{ i=1 }^{ N }{ x_{ i } } \right) }^{ 2 } \right\}
を使うと良いそうです。
ただし、今回は手元にあるデータ自身の統計をとることにフォーカスするため、先の式を用いています。
出典 - Bessel's correction
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
補正加算
新しい方法では、2乗の平均を使用するので、最初のナイーブな方法よりも一時変数の数値が大きくなりがちです。
一般に大きな数字と小さな数字を加算するときには誤差が大きくなります。
そこでカハンの加算アルゴリズムという誤差を減少させる工夫を入れておくと、少し安心です。wikipediaからコピペしてきてみましょう。
うまくoperatorなど実装すれば、コードの変更点を最小にすることができます。
#include <stdio.h>
#include <random>
class Kahan {
public:
Kahan() {}
Kahan(double value) {}
void add(double x) {
double y = x - _c;
double t = _sum + y;
_c = (t - _sum) - y;
_sum = t;
}
void operator=(double x) {
_sum = x;
_c = 0.0;
}
void operator+=(double x) {
add(x);
}
operator double() const {
return _sum;
}
private:
double _sum = 0.0;
double _c = 0.0;
};
class IncrementalStatatics {
public:
void addSample(double value) {
_sum_of_squared += value * value;
_sum += value;
_n++;
}
double variance() const {
return (_sum_of_squared - _sum * _sum / _n) / _n;
}
double avarage() const {
return _sum / _n;
}
private:
Kahan _sum_of_squared = 0.0;
Kahan _sum = 0.0;
int _n = 0;
};
int main()
{
std::vector<double> data;
std::mt19937 engine;
std::uniform_real_distribution<> random(-1.0, 1.0);
for (int i = 0; i < 100000; ++i) {
data.push_back(4.0 + random(engine) + random(engine) + random(engine) + random(engine));
}
IncrementalStatatics statatics;
for (int i = 0; i < data.size(); ++i) {
statatics.addSample(data[i]);
}
printf("avg = %.4f\n", statatics.avarage());
printf("variance = %.4f\n", statatics.variance());
printf("sd = %.4f\n", sqrt(statatics.variance()));
return 0;
}
さらに誤差を抑えるアルゴリズム
有識者の方より、もっと良いアルゴリズムあるよ！という情報を頂きまして、そちらも紹介してみます。
こちらはもういきなりコードにしました。
class OnlineVariance {
public:
void addSample(double x) {
_n++;
double delta = x - _mean;
_mean += delta / _n;
double delta2 = x - _mean;
_M2 += delta * delta2;
}
double variance() const {
// return _M2 / (_n - 1);
return _M2 / _n;
}
double avarage() const {
return _mean;
}
private:
Kahan _M2 = 0.0;
Kahan _mean = 0.0;
int _n = 0;
};
試してみると、確かに他の結果とも一致しました。
ざっくり見てみると、平均値を計算するときに、合計を蓄積しないように上手に回避しています。分散についても同様の考えのようです。
なんで合計を保持しなくて済むかといえば、例えば現在の平均値を $m_1$、一つサンプルが増えたときの平均値を$m_2$とすると、
\begin{align}
m_{ 1 }&=\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ x_{ i } } \\ m_{ 2 }&=\frac { 1 }{ N+1 } \sum _{ i=1 }^{ N+1 }{ x_{ i } } \\ &=\left( \frac { N }{ N+1 } \right) \frac { 1 }{ N } (\sum _{ i=1 }^{ N }{ x_{ i } } +x_{ N+1 })\\ &=\left( \frac { N }{ N+1 } \right) \frac { 1 }{ N } \sum _{ i=1 }^{ N }{ x_{ i } } +\left( \frac { N }{ N+1 } \right) \frac { 1 }{ N } x_{ N+1 }\\ &=\left( \frac { N }{ N+1 } \right) m_{ 1 }+\frac { 1 }{ N+1 } x_{ N+1 }\\ &=\left( \frac { N+1-1 }{ N+1 } \right) m_{ 1 }+\frac { 1 }{ N+1 } x_{ N+1 }\\ &=\left( \frac { N+1 }{ N+1 } \right) m_{ 1 }+\left( \frac { -1 }{ N+1 } \right) m_{ 1 }+\frac { 1 }{ N+1 } x_{ N+1 }\\ &=m_{ 1 }+\left( \frac { -1 }{ N+1 } \right) m_{ 1 }+\frac { 1 }{ N+1 } x_{ N+1 }\\ &=m_{ 1 }+\left( \frac { 1 }{ N+1 } \right) \left( x_{ N+1 }-m_{ 1 } \right)
\end{align}
なるほど確かにこれで先のプログラムと一致します。これは面白いですね。
分散もなんでこれで良いのかは・・・ちょっと統計を勉強することにします。
これによって数字がとても小さく抑えることができ、誤差は大変少なくなることでしょう。
出典
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
マージする
マルチスレッド等の事情で、いくつか個別に計算した値を最後にマージしたいケースもあります。先のIncrementalStataticsの場合は自明ですが、上記OnlineVariance も少し工夫すればマージ可能にすることができそうです。その場合の対応を見ていきます。
まずnは個数なので、単純に加算すれば大丈夫です。
次に変数 M2 ですが、これはvariance()にて最後にM2がnで除算されているだけなので、こちらも単純に足すだけで良さそうです。
最後にmeanですが、これはちょっとだけ工夫が必要です。
直観的には個数から足したあとの割合を考慮して足せば行けそうですが、一応数式で考えてみましょう。
マージする前の値をそれぞれ、$ m_{ a } $ $ m_{ b } $とすると、先の計算と少し似たかんじで展開すると、以下のように落ち着きます。
\begin{align}
m_{ a }&=\frac { 1 }{ N } \sum _{ i=1 }^{ N }{ x_{ i } } \\ m_{ b }&=\frac { 1 }{ M } \sum _{ i=1 }^{ M }{ y_{ i } } \\ m_{ 2 }&=\frac { 1 }{ N+M } \sum _{ i=1 }^{ N+M }{ x_{ i } } \\ &=\frac { 1 }{ N+M } \sum _{ i=1 }^{ N }{ x_{ i } } +\frac { 1 }{ N+M } \sum _{ i=1 }^{ M }{ y_{ i } } \\ &=\frac { N }{ N+M } \frac { 1 }{ N } \sum _{ i=1 }^{ N }{ x_{ i } } +\frac { M }{ N+M } \frac { 1 }{ M } \sum _{ i=1 }^{ M }{ y_{ i } } \\ &=\frac { N }{ N+M } m_{ a }+\frac { M }{ N+M } m_{ b }
\end{align}
となると、さらにマージメンバ関数を追加すると以下のようになります。
class OnlineVariance {
public:
void addSample(double x) {
_n++;
double delta = x - _mean;
_mean += delta / _n;
double delta2 = x - _mean;
_M2 += delta * delta2;
}
double variance() const {
// return _M2 / (_n - 1);
return _M2 / _n;
}
double avarage() const {
return _mean;
}
OnlineVariance merge(const OnlineVariance &rhs) const {
OnlineVariance r;
double ma = _mean;
double mb = rhs._mean;
double N = _n;
double M = rhs._n;
double N_M = N + M;
double a = N / N_M;
double b = M / N_M;
r._mean = a * ma + b * mb;
r._M2 = _M2 + rhs._M2;
r._n = N_M;
return r;
}
private:
Kahan _mean = 0.0;
Kahan _M2 = 0.0;
int _n = 0;
};
これでかなり利便性が上がったのではないでしょうか。
まとめ
ちょっとした工夫ですが、面白い計算方法だと思います。
またコードとしてもシンプルに分散を求められるようになったと思います。
みなさんもガンガン分散を計算してみてはいかがでしょうか
サンプルコードまとめ
https://gist.github.com/Ushio/0203917f5d9e599ede5bdac579afd00e