累積確率分布の「開き」の大きさ

  • Gene-set enrichment analysisの一環としてこちらでコルモゴロフ=スミルノフ検定(KS検定)の話がちらっと出た
  • KS検定は1標本と2標本の区別がある
    • 観察された標本の元の分布が、ある特定の分布ではない、と考えた方がよいかを判断する(1標本のKS検定)
    • もしくは観察された2つの標本の元の分布が同一ではない、と考えた方がよいかを判断する(2標本のKS検定)
  • KS検定について書くためには、経験分布・経験分布関数というのを知っておくのがよい
  • 経験分布・経験分布関数
    • 標本の1番小さい値より小さい範囲の「累積確率」はゼロとする。標本の1番大きい値より大きい範囲の「累積確率」は1とする。途中は、値が来るたびに、標本数nについて\frac{1}{n}だけジャンプして階段状にすることにする。
    • 標本の値をy_1,...,y_nとするとF(x)=\frac{1}{n}\sum_{i}^n \delta(i,x)(\delta(i,x)y_i \le xのときに1でy_i > xのときは0とする)とも書くが(こちら参照)、式を見ると恐ろしげだけれどやっていることは単純
    • ちなみにRなら
y <- rnorm(10)
ecdf.y <- ecdf(y)
plot(ecdf.y)

  • 2標本のKSテストをなぞってみよう
    • もちろんRならks.test()関数を使えばよいだけではあるが
    • 黒いプロットが赤いプロットより小さめ(左)に位置する
# 2標本を作る
n1 <- 20
n2 <- 30
m1 <- 0
m2 <- 1
sd1 <- 1
sd2 <- 1.5
x <- rnorm(n1,m1,sd1)
y <- rnorm(n2,m2,sd2)
# xの方がyより大きいというのが「正しい」ところ
ks.test(x,y,alternative="less")
ks.test(x,y,alternative="two.sided")
ks.test(x,y,alternative="greater")
ecdf1 <- ecdf(x)
ecdf2 <- ecdf(y)
plot(ecdf(x), xlim=range(c(x, y)))
plot(ecdf(y), add=TRUE, col=2)

    • KSテストでは、2つの経験分布関数について、標本x,yの両方の観測値において、2つの経験分布関数の値の差を調べる
    • 赤プロットが黒プロットより上にあるかどうかを知りたいとき("less")には、赤の経験分布値から黒の経験分布値を引いた値の最大値を問題にする
    • 赤プロットが黒プロットより下にあるかどうかを知りたいとき("greater")には、黒の経験分布値から赤の経験分布値を引いた値の最大値を問題にする
    • 赤・黒のどちらが上か下かは気にせず、とにかく違うのかどうかを知りたいときは、赤・黒の値の絶対値の最大値を問題にする
    • この差は次のように、各観測ポイントでの値が算出できる
      • 階段関数なので、その段差について考慮するとそうなる
n.x <- length(x)
n.y <- length(y)
w <- c(x, y)
z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))
plot(w,z)
    • 2つのecdfプロットと「差」を上下に並べるとよくわかる

    • 検定するためには、この「最大値」が「帰無仮説(2標本が同じ分布に由来する)」のときにどのような分布を取るかと引き比べる必要がある
    • これは元の分布の形によらずブラウン運動(ウィーナー過程)から導き出されるので、それに照らしてp値化する。
  • では、次の記事でこの統計量がどんな具合かを調べてみよう