読者です 読者をやめる 読者になる 読者になる

もちもちしている

おらんなの気まぐれブログ

D言語でコルモゴロフ-スミルノフ検定

この記事はD言語 Advent Calendar 2014 24日目の投稿です。

なぜこの題材を選んだのか

私は研究室内で行われたプログラミング講習にて、DDoS検知システムを作りました。その検知ロジックにコルモゴロフ-スミルノフ検定を利用したのですが、Rubyで実装したため、とても遅いのが問題でした。

そこで、ネイティヴで動き、C++よりもはるかに習得が容易だと私のまわりでHotなプログラミング言語であるD言語で実装し直してみることにしました。

コルモゴロフ-スミルノフ検定とは

コルモゴロフ-スミルノフ検定(以下、KS検定)とは、標本Xと標本Yの分布の相違を計る仮説検定のひとつです。

出典:http://sysplan.nams.kyushu-u.ac.jp/gen/edu/MarineStatistics/2007/kougi12/kougi12.pdf

計算手順

1. 標本Xの累積確率分布 S_n(x)と、標本Yの累積確率分布S_m(y)を求める。
2. 上記の2つの累積確率分布の差の絶対値の最大値であるKS統計量Dを求める。


D=max|S_n(x)-S_m(y)|

3. Dから検定統計量\chi ^2、相関係数pを算出する。


\chi ^2=4 D ^2 \frac{sum(x)sum(y)}{sum(x)+sum(y)}


p=min(1, qchi2(\chi ^2, 2) ^2)

4. 帰無仮説の採否を決める(有意水準5%の場合)

p > 0.05の場合 : 帰無仮説の採択
p <= 0.05の場合 : 帰無仮説の棄却

カイ2乗分布を求める

qchi2はカイ2乗分布の上側累積確率を求める関数です。実装は以下の通りです。

//正規分布の下側累積確率
auto p_nor(immutable double z)
in {
  assert(z >= 0.0 && z <= 3.09);
}
out(probability) {
  assert(probability >= 0.0 && probability <= 0.499);
}
body {
  double t, p;
  t = p = z * exp(-0.5 * z * z) / sqrt(2 * PI);

  for(size_t i = 3; i < 200; i+= 2) {
    immutable auto prev = p;
    t *= z * z / i;
    p += t;
    if (p == prev) { return 0.5 + p; }
  }
  return (z > 0.0) ? 1.0 : 0.0;
}

// 正規分布の上側累積確率
auto q_nor(immutable double z)
in {
  assert(z >= 0.0 && z <= 3.09);
}
out (probability) {
  assert(probability >= 0.001 && probability <= 0.5);
}
body {
  return 1.0 - p_nor(z);
}

//上側累積確率
auto q_chi2(immutable int df, immutable double chi2)
in {
  assert(df >= 1 && df <= 300);
  assert(chi2 >= 0.00003927 && chi2 <= 366.844);
}
out(probability) {
  assert(probability >= 0.0 && probability <= 1.0);
}
body {
  double s, t;
  if ((df & 1) == 1) { //自由度が奇数
    const auto chi = sqrt(chi2);
    if (df == 1) { return 2.0 * q_nor(chi); }
    s = t = chi * exp(-0.5 * chi2) / sqrt(2 * PI);
    for (size_t i = 3; i < df - 1; i += 2) {
      t *= chi2 / i;
      s += t;
    }
    return 2 * (q_nor(chi) + s);
  } else { //自由度が偶数
    s = t = exp(-0.5 * chi2);
    for (size_t i = 2; i < df - 1; i += 2) {
      t *= chi2 / i;
      s += t;
    }
    return s;
  }
}

KS検定を行う

KS検定を行うプログラムは以下の通りになります。

//累積和を求める 
auto cumsum(immutable double[] a) {
  auto b = cast(double[])a;
  foreach(i; 1 .. a.length) {
    b[i] = a[i] + a[i - 1];
  }

  return b;
}

auto ks_test(immutable double[] x, immutable double[] y)
in  {
  assert(x.length == y.length);
 }
out(results) {
  assert(results[2] >= 0.0 && results[2] <= 1.0);
}
body {
  immutable auto sum_x = x.reduce!("a + b");
  immutable auto sum_y = y.reduce!("a + b");
  auto cum_x = cumsum(x).map!(s => s / sum_x);
  auto cum_y = cumsum(y).map!(s => s / sum_y);

  double[] cum = new double[cum_x.length];

  for(size_t i = 0; i < cum.length; i++) {
    cum[i] = abs(cum_x[i] - cum_y[i]);
  }

  // 累積分布関数の差の最大値
  immutable auto d = cum.reduce!(max);
  // 検定統計量
  immutable auto chi = 4 * d * d * sum_x * sum_y / (sum_x + sum_y);
  // 相関係数(p値)
  immutable auto p = [1.0, pow(q_chi2(2, chi), 2.0)].reduce!(min);

  return [d, chi, p];
}

void main() {
  immutable auto x = [1.0,2.0,1.0,3.0,2.0,3.0,3.0,2.0,7.0,11.0,10.0,9.0,13.0,13.0,22.0,17.0,23.0,20.0,17.0,14.0,13.0,5.0,2.0,1.0,1.0,1.0];
  immutable auto y = [0.0,1.0,2.0,2.0,4.0,5.0,6.0,8.0,10.0,7.0,16.0,17.0,17.0,13.0,19.0,13.0,18.0,10.0,4.0,6.0,6.0,5.0,1.0,3.0,0.0,0.0];

  immutable auto p = ks_test(x, y)[2];

  // 帰無仮説の採否(有意水準5%)
  // p <= 0.05 (帰無仮説の棄却:二つの分布に差があると言える)
  // p > 0.05 (帰無仮説の採択:二つの分布に差があると言えない)
  (p <= 0.05).writeln();

}

実装は以上となります。全体の実装は以下のサイトに置いてあります。

KS_test.d · GitHub

感想

Rubyで実装したときは楽に実装できましたが、D言語でも同じくらい楽に実装できました。D言語最高〜〜〜✌('ω')