pooneilの脳科学論文コメント

■ 適合度の尤度比検定ってこれで良いの?

(追記：latexでPDFファイルを作成しました。そちらの方が数式が読みやすいかと思います。こちらから：loglikelihood.pdf)

宿題持ってくるんじゃねー、とか言われたりしないのでブログは良い。
以下の問題を解こうとして、ネットを調べてたんですが、よくわからないので自作してみました。これでいいんでしょうか、先生(<-だれよ)：

こういう問題です：
Observed dataがN個あってそれのヒストグラムを作りました。
このヒストグラムをfittingしてやるために二種類のモデルによるヒストグラムを作りました。
モデルH0とモデルH1があって、H1はH0のモデルにパラメータをひとつ付加したもの。それぞれのモデルからデータをgenerateしてやって、ヒストグラムを作る。
モデルH0によるfittingの適合度とモデルH1によるfittingの適合度とを比べて有意差があるか知りたい。
そこで、ふたつの適合度の尤度比検定をしてやろうというわけです。
それぞれのfittingの尤度関数さえ計算できれば、モデル間で差がないとする帰無仮説において「-2*loglikelihoodの差」が自由度1(H0とH1での自由度の差)のカイ二乗分布に従うことを使って検定できます。
(なお、全部ヒストグラムではなくて確率密度関数にしても成り立つはず。あと、要はモデル選択なので、AICとかBICでやるかという話もあるけど、とりあえず簡単のため尤度比検定で。)

Obserbed dataはN個：

[data(1),data(2),data(3),....data(N)]

ヒストグラムのbinはn個：

[x(1),x(2),x(3),....x(n)]

Observed data、
モデルH0でのbest-fitted data、
モデルH1でのbest-fitted data、
のそれぞれのヒストグラムでの各binの頻度：

[ NO(x(1)), NO(x(2)),..., NO(x(n))]
[NE0(x(1)),NE0(x(2)),...,NE0(x(n))]
[NE1(x(1)),NE1(x(2)),...,NE1(x(n))]

それぞれのヒストグラムのデータの総数は

∑( NO(x(i)))=N
∑(NE0(x(i)))=N
∑(NE1(x(i)))=N

です。(Fitting時に同じになるようにnormalizeしてある。)

まずそれぞれのfittingの尤度関数を作ってみる。
まずObserved dataのある1個data(j)がヒストグラムのbin x(i)に落ちるとすると、
data(j)がモデルH0の分布のbin x(i)に落ちる確率は

NE0(x(i))/∑(NE0(x(i)))=NE0(x(i))/N

と書けます。
だから、尤度関数L(H0)は全observed dataの数だけこれを掛け合わせたものです。
Observed dataの[data(1),data(2),data(3),....data(N)]
がそれぞれ落ちるbinがたとえば
x(m),x(n),x(q),x(r)
とかだとすると、

L(H0)=(NE0(x(m))/N)*(NE0(x(n))/N)*(NE0(x(q))/N)*...*(NE0(x(r))/N)
=NE0(x(m))*NE0(x(n))*NE0(x(q))*...*NE0(x(r)/(N^N)

といったかんじになります。
これをlog-likelihoodに変換すると、

LL(H0)=ln(NE0(x(m)))+ln(NE0(x(n)))+ln(NE0(x(q)))+...+ln(NE0(x(r))-N*ln(N)

です。
全observed dataはx(i)のbinのどこかに落ちるから、
たとえば、x(1)のbinに落ちるobserved dataはNO(x(1))個ある、というふうに整理すると、

LL(H0)=NO(x(1))*ln(NE0(x(1)))+NO(x(2))*ln(NE0(x(2)))+...+NO(x(n))*ln(NE0(x(n)))-N*ln(N)

となります。これを∑を使ってbinごとに(i=1:n)足し合わせるように表記すると

LL(H0)=∑(NO(x(i))*ln(NE0(x(i))))-N*ln(N)

と書けます。
同様にH1のモデルのときは

LL(H1)=∑(NO(x(i))*ln(NE1(x(i))))-N*ln(N)

となります。
あとは2LLを作るだけ。二番目の項が消えます。

2LL = -2*(LL(H1)-LL(H0))
=2*∑(NO(x(i))*ln(NE0(x(i))))-2*∑(NO(x(i))*ln(NE1(x(i))))

∑はどちらも共通のbin x(i)で足し合わせているから合体できます。

2LL = 2*∑(NO(x(i))*ln(NE0(x(i))/NE1x(i)))

ということで計算できました。
これで合ってるのか答えがネットを探していても見つからないのだけれど、G-testやKL diveregenceの値と似ているから、間違った方向には行ってないでしょう。ということで検算のためにG-testの式に近づけてみる：
上記の2LL式を変形してやると、

2LL= 2*∑(NO(x(i))* ln(NE0(x(i))/NE1(x(i))*NO(x(i))/NO(x(i))))
   = 2*∑(NO(x(i))*(ln(NE0(x(i))/ NO(x(i)))-ln(NE1(x(i))/ NO(x(i)))))
   =-2*∑(NO(x(i))*(ln( NO(x(i))/ NO(x(i)))-ln( NO(x(i))/NE1(x(i))))
   = 2*∑(NO(x(i))*(ln( NO(x(i))/NE1(x(i))))-2*∑(NO(x(i))*(ln(NO(x(i))/NE0(x(i))))
   = G(H1)-G(H0)

となって、二つのG検定値の差となっています。なんてこったorz
ちなみにモデルH0でobserved dataをfittingしたときのG検定のG-statisticsは、

G（H0)=2*∑(NO(x(i))*(ln(NO(x(i))/NE0(x(i))))

となります。

というわけでたぶんこんな長々と計算しなくても、二つのfittingをして、G-statisticsを計算して、その差をカイ二乗検定すれば良かったというオチ、のようです。じつははじめに計算をしたときは、G statisticsではなくてカイ二乗の方を使ったのだけれど、カイ二乗の差をまたカイ二乗検定、ってなんかへんではないか?と思って上のような計算をしてみた次第。あと、G-test自体がobserved dataとfitted dataとのあいだでの尤度比検定なのは知っていたのですが、二つのモデルの比較で単純にそれらをさし引いて良いかどうかがわからなかったわけです。以上の計算からすると、最尤推定の考えからしてもさし引いて良い、ということになりそうですが。

ま、いつも通り自分で疑問出して自分で納得してるってかんじなのですが、まだ納得しているわけではないんです。というのも、計算される値が大きすぎるし、ググっても、関係してくる項目が見つからない。まだなんか間違えている気がします。

また、G検定はいわばKL divergenceの離散バージョンと式の上では同じですから(ところでこのことをwebで探しても明確に書かれているものを見つけられない)、情報幾何で使われるイメージ化の方法が使えます。ここでしている-2LLのカイ二乗検定というのは、[Observed data -> model H0でのbest fitted data]の距離と[Observed data -> model H1でのbest fitted data])の距離との差を取っているということになります。たとえば「神経回路網とEMアルゴリズム」とかにあるようなとり扱いをすると、observed dataの集団を多様体の中のある一点として、H0のモデルに属した曲面とH1のモデルに属した曲面とがあって、observed dataの点からそれぞれの曲面に垂線を下ろした当たったところがそれぞれのモデルでのbest-fitted dataの集団の位置で、その垂線(正確にはm測地線)の距離がKL-diveregenceです。

さて、そうするとわたしがわからないのは、そういったまっすぐでない空間での二つのKL-divergenceをたんにさっ引いてよいのかっていう問題になるのかも知れません。ていうかそこまでいくと、2LLの原理を調べろ、ってことでネイマンーピアソンまで戻って勉強しないといけない。こういうことやってるとJSTORにある昔の統計学の論文とか読み出してどんどんはまることになるので、このへんまでにしておきます。