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();
}
実装は以上となります。全体の実装は以下のサイトに置いてあります。
