標準偏差(ExcelのSTDEV関数)の算出をC言語で実装する

標準偏差の算出を、チープなシステムにC言語で実装した際の覚書です。

標準偏差の算出

まず、分散(var)を求めます。
分散をデータ数-1で割ったルートが標準偏差(stdev)となります。

var = Σ[i=1→n]((data[i] - mean)^2)
stdev = √(var / (n - 1))

var: 分散
mean: 平均値
stdev: 標準偏差

実装

1. 対象のデータをすべてメモリに取得出来る場合

標準偏差を算出する関数は、データとデータ数を引数に取ります。

stdev.c

 #include <math.h> // sqrt()

 double stdev(int *data, int n) {
   int sum = 0;
   double mean;
   double var = 0.0;
   int i;

   for (i = 0; i < n; i++) {
     sum += data[i];
   }
   mean = (double)sum / n;

   for (i = 0; i < n; i++) {
     var += (data[i] - mean) * (data[i] - mean);
   }
   return sqrt(var / (n - 1));
 }

使い方は、

// データ数は20
int data[] = {
    0, 0, 0, -1, -1, 2, -1, -1, 2, 0,
    -2, 0, -2, 1, 0, 1, -2, 0, 0, 2
};
double sd = stdev(data, 20); // => sd: 1.252366

2. データ数が多くて、すべてをメモリに保持できない場合

私の実装システムでは、データ数が多く、データを保持することができませんでした。
分配法則 (a - b)^2 = a^2 - 2ab + b^2 でvarを展開すると

var = Σ[i=1→n]((data[i] - mean)^2)
    = Σ[i=1→n](data[i]^2 - 2(data[i] * mean) + mean^2)
    = Σ[i=1→n](data[i]^2) - (2 * Σ[i=1→n](data[i]) * mean) + n * mean^2

となり、データの総和とデータを2乗したものの合計から、標準偏差を算出できます。
これなら、データを保持するのと比べて、必要なメモリが少なくて済みます。

stdev.c

#include <math.h> // sqrt()

double stdev(int sum, int sqsum, int n) {
  double mean;
  double var;

  mean = (double)sum / n;

  var = sqsum - (2 * sum * mean) + (mean * mean * n);
  return sqrt(var / (n - 1));
}

使い方は、

// データ数は20
int data[] = {
    0, 0, 0, -1, -1, 2, -1, -1, 2, 0,
    -2, 0, -2, 1, 0, 1, -2, 0, 0, 2
};
int sum = 0;    // データの合計
int sqsum = 0;  // データを2乗したものの合計

for (i = 0; i < 20; i++) {
  sum += data[i];
  sqsum += data[i] * data[i];
}

// 渡すときに、データは保持しておく必要はなし
double sd = stdev(sum, sqsum, 20); // => sd: 1.252366

メモ

• varをnで割ると、ExcelのSTDEVP関数となる
  sqrt(var / n); // STDEVP
• 対象が標本データ(サンプリング)の場合はSTDEV、母集団そのもの(全数)の場合はSTDEVP