分散を計算するときの桁落ち対策

分散 = 二乗平均 - 期待値の二乗

という式を使って分散を計算すると、右辺の2つの項の値が非常に近い場合、桁落ちが発生する。

たとえば以下のようなデータを考える。

x₁ = 10¹⁰
x₂ = 10¹⁰ + 1

x₃ = 10¹⁰ + 2

x₄ = 10¹⁰ + 3

x₅ = 10¹⁰ + 4

この場合、二乗平均も期待値の二乗も10²⁰程度である。

両者の差は2だが、倍精度浮動小数点の有効桁数は15桁程度なので、この差の情報は失われてしまう。

分散は平行移動しても変わらないので、平均が0になるように中心化を行う。

x₁ = -2
x₂ = -1

x₃ = 0

x₄ = 1

x₅ = 2

すると、二乗平均は2、期待値の二乗は0となり、差の情報は失われない。

以下のプログラムを実行してみると、桁落ちの様子を確認できる。

#include <iostream>
#include <vector>
#include <iomanip>

using namespace std;

// データを中心化する
vector<double> centerize(const vector<double> &v) {
  int n = v.size();
  double avg = 0.0;
  for (auto &x : v)
    avg += x;
  avg /= n;
  vector<double> w(n);
  for (int i = 0; i < n; i++)
    w[i] = v[i] - avg;
  return w;
}

// 分散を計算する
double calc(const vector<double> &v) {
  int n = v.size();
  double avg_sq = 0.0;
  double avg = 0.0;
  for (auto &x : v) {
    avg += x;
    avg_sq += x * x;
  }
  avg_sq /= n;
  avg /= n;

  return avg_sq - avg * avg;
}

// データ生成
vector<double> gen() {
  vector<double> v;
  for (int i = 0; i < 5; i++)
    v.push_back(1e10 + i);
  return v;
}

int main(int argc, char *argv[])
{
  vector<double> v = gen();
  cout << fixed << setprecision(15);
  cout << calc(v) << endl;
  cout << calc(centerize(v)) << endl;
  return 0;
}

以下実行結果。

$ ./main       
0.000000000000000
2.000000000000000