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の累積確率分布 と、標本Yの累積確率分布を求める。
2. 上記の2つの累積確率分布の差の絶対値の最大値であるKS統計量Dを求める。
3. Dから検定統計量、相関係数pを算出する。
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(); }
実装は以上となります。全体の実装は以下のサイトに置いてあります。