KS_test.d

import std.stdio;
import std.math;
import std.algorithm;

//正規分布の下側累積確率
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;
  }
}

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();

}