ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。

ベイズ統計でよく出てくる「正規化定数を一旦無視して計算」について

この記事はベイズ統計の本を読んでるとたまに出会うであろう「正規化定数を一旦無視して計算」するやつについて, 「なんで無視できるんだろう」と思った人(いるのか?)向けの投稿です.

最初の例:分割表

ベイズの定理の説明でよく出てくる2×2の分割表を, 私の好みに基づき次のように書いた.

曝露(X) 疾病(Y) 数(Z)
0 0 z_{00}
0 1 z_{01}
1 0 z_{10}
1 1 z_{11}

この表では曝露や疾病がない場合を0, ある場合を1とコード化している.

また, ここでの「曝露」と「疾病」は例としてイメージしやすそうな言葉を入れただけで以降で特に疫学的な話題を扱うわけではない.

全体の数を  N= z_{00}+z_{01}+z_{10}+z_{11} と置く.

疾病あり(Y=1)の集団から無作為に一人とってきたとき、その人が曝露あり(X=1)だった確率を考えてみる.

 \displaystyle P(X=1|Y=1)=\frac{z_{11}}{z_{01}+z_{11}}.

右辺の分母は疾病ありの人全員の数で, 分子はそのうち曝露ありの人の数である. 分子と分母をN で割って, 次のように書いてもいい.

 \displaystyle P(X=1|Y=1)=\frac{z_{11}/N}{(z_{10}+z_{11})/N} \\
\displaystyle =\frac{P(X=1, Y=1)}{P(Y=1)}.

右辺は条件付き確率の定義そのものである.

もう少し変形してみるのでお付き合いください. 右辺の分子に  1=(z_{10}+z_{11})/(z_{10}+z_{11})をかけて,

 \displaystyle P(X=1|Y=1)\\
\displaystyle =\left( \frac{z_{11}}{z_{10}+z_{11}} \cdot \frac{z_{10}+z_{11}}{N} \right) \big/ \left(z_{10}+z_{11})/N\right) \\
\displaystyle = \frac{P(Y=1|X=1)P(X=1)}{P(Y=1)}.

この等式は一般の確率変数の場合,

 \displaystyle P(X=x|Y=y) = \frac{P(Y=y|X=x)P(X=x)}{P(Y=y)}

にはベイズの定理と呼ばれる.

ここであらためて最初に書いた式( P(X=1|Y=1)=z_{11}/(z_{01}+z_{11}))を見ると, 条件付き確率に関する計算は比の関係だけ考えて, あとから全体(z_{01}+z_{11}N)で割って0から1の範囲に収まるようにしてやると, だいぶ楽になるように感じると思う.

そこからさらに進んで, ベイズの定理は条件付き確率の定義からほぼ明らかな計算に名前をつけたに過ぎないとまで感じて貰えれば幸甚である.

指数型分布族の共役事前分布

数学が絡む授業の(嫌な)あるあるとして, 最初のほう割合を計算してみようみたいな簡単な例からはじまったのに, 途中から急に積分とか出てきて難しくなる, というのがある気がする. 残念なことに, ここでもそれをやる.

さっきはすべての分布がわかっているものとして計算したが, ここでは未知の潜在変数があってそれを推定したいというような動機がある.

さて, 指数型分布族とは確率(密度)関数が0でない値を取る範囲において, ベクトル値関数(ベクトルを返すような関数) f(\theta) g(x) を用いて,

\displaystyle p(x|\theta) = \exp( \langle f(\theta) , g(x)\rangle)

と表せるような分布のことである. ここでは f(\theta) g(x)内積 \langle f(\theta) , g(x)\rangle と表記した.

これに対して, もう一つの確率分布,

\displaystyle \pi(\theta | \phi) = \frac{1}{Z(\phi)}\exp( \langle \phi , f(\theta)\rangle)

を考える.

ベイズ統計では考えたい対象の変数すべてに確率分布を設定する.

ここでは  x をデータとして観測される変数,  \theta を推定したい未知の潜在変数とし, p(x|\theta) \pi(\theta | \phi) をそれぞれ x\theta の分布として設定する.

 x を所与としたときの  \theta の分布(これを \hat pi(\theta | x)と置くことにし, 事後分布と呼ぶ)は, 条件付き確率から次のように書ける.

\displaystyle \hat \pi(\theta | x) \propto \pi(\theta | \phi) p(x|\theta).

この  \propto は左辺が右辺に比例(適切な定数を掛けたら一致)することを意味する.

右辺を右辺の「全体で割って0から1の範囲に収まるようにする」ことを目標に, 定義域全区間に渡る積分 \hat Z(\phi) と置く),

 \displaystyle \hat Z(\phi)= \int \pi(\theta | \phi) p(x|\theta) \, d \theta

を考える. ここで p(x|\theta) \pi(\theta | \phi) の定義を思い出すと,

 \displaystyle \hat Z(\phi)= \frac{1}{Z(\phi)} \int \exp( \langle f(\theta) , g(x)\rangle) \cdot \exp( \langle \phi , f(\theta)\rangle) \, d \theta \\ \displaystyle = \frac{1}{Z(\phi)}\int \exp( \langle \phi +g(x) , f(\theta)\rangle) \, d \theta \\
\displaystyle = \frac{1}{Z(\phi)}Z(\phi+g(x)),

ここから,

\displaystyle \hat \pi(\theta | x) = \frac{1}{\hat Z(\phi)} \pi(\theta | \phi) p(x|\theta) \\
\displaystyle = \frac{1}{Z(\phi+g(x))} \langle \phi+g(x), f(\theta)\rangle \\
\displaystyle = \pi (\theta| \phi + g(x))

という結果を得る.

この結果から指数型分布族には, パラメータ \phi (事前分布のパラメータは特にハイパーパラメータと呼ばれる)に統計量 g(x) (この統計量は十分統計量と呼ばれる)を足すだけで事後分布が求められるような事前分布が存在することと, その事前分布が(上と同じ記号で)\pi(\theta | \phi) であることがわかる. この事前分布を共役な事前分布と呼ぶ.

正規化定数を一旦無視して計算する例:ベルヌーイ分布

共役な事前分布の簡単な例として, ベルヌーイ分布に対するベータ分布を考えてみる.

ベルヌーイ分布の確率関数は x が0または1のときに,

\displaystyle p(x|\theta) = \theta^x (1-\theta)^{1-x}

それ以外の場合に0を取る.

ベータ分布の確率密度関数は,

\displaystyle \pi(\theta| a, b ) = \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1}

である. ここで B(a,b) はベータ関数とした.

指数型分布族一般の共役事前分布を考えたときの先の表記に合わせる場合は,

\displaystyle \phi=(a-1, a+b-2)',

\displaystyle f(\theta) = \left(\log (\theta/(1-\theta)), \log(1-\theta)\right)',

\displaystyle g(x)=(x,1)'

とすれば良い.

一方で, 「正規化定数を一旦無視して計算」すると,

\displaystyle \hat \pi(\theta | x) \propto \int p(x|\theta)  \pi(\theta| a, b ) \, d \theta \\
= \theta^{a+x} (1-\theta)^{b+1-x}

であるが, 右辺を全区間にわたって積分すると

\displaystyle \int p(x|\theta) \, \pi(\theta| a, b ) \, d \theta \\
\displaystyle =\int  \theta^x (1-\theta)^{1-x} \cdot \theta^{a} (1-\theta)^{b}  \, d \theta \\
\displaystyle = \int  \theta^{a+x} (1-\theta)^{b+1-x} \\
\displaystyle = B(a+x, b+1-x)

である. こうしてやってもやはり事後分布は, \pi (\theta| a+x, b+1-x) とハイパーパラメータを置き換えるだけで得られることがわかる.