Dufor, J.M., Farhat, A. (2002) Exact Nonparametric Two-Sample Homogeneity Tests. In: Huber-Carol C., Balakrishnan N., Nikulin M.S., Mesbah M. (eds) Goodness-of-Fit Tests and Model Validity. Chap.33.
仕事の都合でがんばって読んだ奴。2標本Kolmogov-Smirnov検定が、離散変数のときにどうなるのかを知りたくて、なにか書いてあるかなあと思ったのである。
1. イントロダクション
2つの無作為標本\(X_1, \ldots, X_n\)と\(Y_1, \ldots, Y_m\)があるとしよう。累積分布関数(cdf)を $$ F(x) = P[X_i \leq x], \ \ i = 1, \ldots, n $$ $$ G(x) = P[Y_j \leq x], \ \ j = 1, \ldots, m $$としよう。cdfにこれ以上の制約はかけない。それは連続かもしれないし離散かもしれない。
さて、\(H_0: F=G, H_1:F \neq G\)の検定をしたい。どうするか。
多くの人が思いつくのは、2標本Kolmogov-Smirnov(KS)検定、ないしCramer-von Mises(CM)であろう。ほかにも、pdfのカーネル推定量の間のL1距離なりL2距離なりに基づく数え上げ検定とか、2標本の平均差に基づく検定とかが考えられる。
最後に挙げた平均差の奴を別にすれば、検定統計量の正確分布が手に入るのは標本サイズが小さいときだけで、たいていは漸近分布ベースの表を使う羽目になり、検定量が低くなったりする。
そこで本論文では、(1)cdfを問わずに使えて(2)有限標本において正確な、検定手続きを提案する。
なお、KS統計量やCM統計量だって分布フリーじゃんと思ったあなた、分布フリーなのは観察が連続分布にiidに従うときだけなんです。カーネルタイプの密度推定量に基づく統計量は、観察が連続分布にiidに従っていたとしてもなお分布フリーでない。この難点は数え上げを使うと解決できるけど、それは連続分布のときだけだ。というわけで、(1)(2)の両方を満たす方法はこれまでなかった。
2. 検定統計量
先行する検定統計量について。
- KS検定。$$ KS = sup_x |F_n(x) – G_m(x)| $$ ただし\(F_n(x), G_m(x)\)は経験分布関数(edf)。\(H_0: F=G\)でかつ\(F\)が連続分布のときには分布フリーになる。離散のときはそうでない。
- 2標本CM検定。$$ CM = \frac{mn}{(m+n)^2} \left( \sum_{i=1}^n [F_n(X_i) – G_m(X_i)]^2 + \sum_{j=1}^m [F_n(Y_j) – G_m(Y_j)]^2 \right) $$ [えーっと… すべての観察値について、それを\(X\)のedfに放り込んだときと\(Y\)のedfに放り込んだときとの差の二乗を求めて、平均して\(mn/(m+n)\)を掛けるわけね]
これも、\(F\)が連続だったら\(H_0\)の下で分布フリー。 - カーネル・ベースのPDF推定量の間の距離。まず\(F\)のpdf推定量を $$ f_n(x) = \frac{C_X}{n} \sum_{i=1}^n K[C_X(x-X_i)] $$ $$ C_X = n^{1/5} / (2S_X) $$ とする。\(K(x)\) は \(|x| \leq 1\)のときに\(1/2\), そうでないときに\(0\)となる関数。\(S_X\)はふつうのSD推定量である。同様に\(G\)のpdf推定量を\(g_m(x)\)とする。距離としては $$ \hat{L}_1 = \sum_{i=1}^n |f_n(X_i) – g_m(X_i)| + \sum_{j=1}^m |f_n(Y_j) – g_m(Y_j)|$$ $$ \hat{L}_2 = \left( \sum_{i=1}^n (f_n(X_i) – g_m(X_i))^2 + \sum_{j=1}^m (f_n(Y_j) – g_m(Y_j))^2 \right)^{1/2} $$ $$ \hat{L}_\infty = sup_x |f_n-g_m| $$ がある[3本目のやつは最大値の式に書き換えられるんだけどめんどくさいのでメモ省略]。\(F\)と\(G\)が離散だったら、pdfの代わりに確率質量関数 \(p(x) = P[X = x], q(x) = P[Y=x]\)を使う。これは\(F\)が連続だろうが離散だろうが分布フリーでない。
- 標本平均の差 \(\hat{\theta}_1 = \bar{X}-\bar{Y}\)。さらに、平均が違うかどうかだけではなくて、分散が違うか、尖度が違うか、歪度が違うかという検定についても考える。それぞれの検定統計量を\(\hat{\theta}_2, \hat{\theta}_3, \hat{\theta}_4\)とする[面倒くさいので数式省略]。
3. 正確ランダム化数え上げ検定
[ここから苦手分野に突入するので、ほぼ逐語訳…]
Dwass(1957)の手続き[ \( \hat{\theta}_1\) のこと] を別にすれば、前節で述べた検定はすべて、不完全にtabulateされた帰無分布を含む、ないし、有限標本において分布フリーでない。その結果、後者は恣意的に大きなサイズの歪みに繋がる。
既知のサイズの有限標本において分布フリーな検定を得るという観点からは、まず次の点に注目する。統計量 \( KS, CM, L_1, L_2, L_\infty, \hat{t}, \hat{\theta}_1, \hat{\theta_2}, \hat{\theta}_3, \hat{\theta}_4 \) に基づく、(いかなる標本サイズにおいても)真に分布フリーな検定を手に入れるには、\(m+n\)個の観察\(X_1, \ldots, X_n, Y_i, \ldots, Y_m\)を、すべての可能な方法で(等確率に)数え上げて得られる分布を考えれば良い。こうした数え上げは、未知の分布\(F\)がどうであれ、帰無仮説\(H_0\)の下で等しい確率を持つ。従って、その数え上げ分布(すなわち、その観察の順序統計量の下での条件分布)から得られた正確臨界値を使って\(H_0\)を棄却する検定は、have the same level conditionally (on the ordered statistics) as well as unconditinally である。[←学力不足でよく意味がつかめないけど、まあとにかく、permutationで得た分布で正確検定できますよってことなんだろうな]
\(T\)がpivotalな検定統計量を表すとして(すなわち、その分布が帰無仮説のもとで未知のパラメータに依存しないとして)、以下の方法でMC検定を進めることができる。
観察された標本から求めた検定統計量を\(T_0\)とする。\(T_0\)が大きな値を持つときに帰無仮説が棄却される場合、サイズ\(\alpha\)の棄却域は\(G(T_0) \leq \alpha\)と表現できる。ただし、\(G(x) = P[T \geq x|H_0]\)はp値関数である。
まず、指定された帰無分布\(F_0\)から\(N\)個の独立標本 \(X_{i1}, \ldots, X_{in}, Y_{i1}, \ldots, Y_{im}, 1\leq i \leq N\)を生成する。ここから、\(N\)個の独立な実現値 \( T_i = T(X_{i1}, \ldots, X_{in}, Y_{i1}, \ldots, Y_{im}), 1\leq i \leq N\) を生成する。ここから経験p値関数 $$ \hat{p}_N(x) = \frac{N \hat{G}_n(x)+1}{N+1} $$ を得る。ただし $$ \hat{G}_N(x) = \frac{1}{N} \sum_{i=1}^N \mathbf{1}_{[0, \infty)} (T_i-x) $$ $$ \mathbf{1}_A{x} = \begin{cases} 1, x \in A \\ 0, x \notin A \end{cases} $$ [写経ばかりでもなんなので意味を考えると、\(\hat{G}_N(x)\)は検定統計量の実現値\(T_1, \ldots, T_N\)のうち\(T_i -x\)が正である割合を表している。おおざっぱに言うとこれがp値になるわけね。なぜこれそのものをp値とせずにちょっと修正しているのか、よく分からないけれど]
MC臨界域は $$ \hat{p}_N(T_0) \leq \alpha $$ と定義できる。この\( \hat{p}_N(T_0) \) は \(G(T_0)\)の推定値と解釈しても良い。\(T\)が連続分布ならば下式が成り立つことを示せる: $$ P[\hat{p}_N(T_0) \leq \alpha | H_0] = \frac{I(\alpha(N+1))}{N+1}, \ \ 0 \leq \alpha \leq 1$$ ただし\(I[x]\)は\(x\)を越えない最大の整数である。従って、\(\alpha(N+1)\)が整数になるように\(N\)を選んでおけば、臨界域は\( \hat{p}_N(T_0) \leq \alpha \) は臨界域 \(G(T_0) \leq \alpha\)と同じサイズになる。こうして得られるMC検定は、\(N\)に関わらず理論的に正確である。
[学力不足で話の流れがつかみづらいんだけど… いま手元の2群の標本から検定統計量\(T_0\)を求めたとして、そのモンテカルロ検定を考えようという訳ね。そのために、2群の標本をあわせたCDFから同サイズの標本を抽出して検定統計量を得るというのを\(N\)回繰り返す。で、\(N\)個のうち \(T_0\)より値が大きい奴が占める割合(をちょっと修正した奴)をp値とし、それが\(\alpha\)を下回っていたら2群が異なると判断する、ということであろう]
[ブートストラップ検定とのちがいについての説明。省略]
上で定義した\( \hat{p}_N(x) \)は、連続分布からの統計量に基づく検定において成り立つ。この場合[←??? 文脈からいえば離散分布の場合のことを指していると思うのだが…]、タイが生じる確率はゼロではない。しかしMC検定の手法は、以下のランダム化したタイ・ブレーク手続きを使えば、離散分布にも適用できる。
一様分布する\( N+1 \)個の変量\( U_0, U_1, \ldots, U_N \)を、\( T_i \)とは独立に抽出する。ペア\( (T_i, U_i) \)を、以下のlexicographicな順序で並べる。$$ (T_i, U_i) < (T_j, U_j) \Leftrightarrow [T_i < T_j \ or \ (T_i = T_j \ and \ U_i < U_j) ] $$ [ひらたくいうと、基本的には\(T\)の順に並べます、もしタイだったら\(U\)の順にします、ってことね]
次に、連続分布の場合と同様に$$ \tilde{p}_N(x) = \frac{N \tilde{G}_n(x)+1}{N+1} $$を求める。ただし、$$ \tilde{G}_N(x) = 1 – \frac{1}{N} \sum_{i=1}^N \mathbf{1}_{[0, \infty)} (T_i-x) + \frac{1}{N} \sum_{i=1}^N \mathbf{1}_{[0]} (T_i-x) \mathbf{1}_{[0,\infty)} (U_i – U_0) $$ [こんどは、\(N\)個の\(T_i\)のうち (a)\(T_0\)より小さい奴が占める割合, (b)\(T_0\)と等しくてかつ\(U_0\)より\(U_i\)が大きい奴が占める割合、を求めて、1 – (a) + (b) を求める、ってことね]
結果として得られる臨界域 \( \tilde{p}_N (T_0) \leq \alpha \) は、\( \alpha(N+1) \)が整数であれば領域\(G(T_0) \leq \alpha\)と同じレベルになる。もっと正確に言うと[…略]。
もし、帰無仮説によって「無作為標本は交換可能な変数からなっている」ということが保証されており、かつ検定統計量の値が大きいときに帰無仮説が棄却されるのならば、その仮説のMC検定は次の5段階で実行される。
- 観察された標本から検定統計量の値を求める。これを\(T_0\)としよう。
- 標本のすべてのpermulationsから、\(N\)個のpermutationsを復元無しで無作為抽出する。
- それぞれのpermutateされた標本から検定統計量を再計算する。これを\(T_1, \ldots, T_N\)としよう。
- \(\{T_0, T_1, \ldots, T_n\}\)における\(T_0\)のランクを\(R_0\)とする(タイがある場合は無作為化手法を使う)。帰無仮説のMC検定のp値は\(1-R_0/(N+1)\)となる。
- 選択したレベルに従って判断する。
4. シミュレーション研究 [パス]
5. 結論 [略]
… 苦労して読んだけど、2標本Kolmogorov-Smirnov検定量って離散変数のときにどうなっちゃうのか、結局よくわからないまま終わってしまった… がっくり… 読む文献を間違えたかも…