はじめに
あまり格好のよいことではないのだけど、内省のために書いた文章がそこそこの分量に達したし、誰かにとって有益になるかもしれないと思ったので自分語りをしてみます。軽い話ではないので、そういうのが好みではない人は読み飛ばしてください。
元々は自分向けの文章なので、表現が雑な箇所もあります。また、敬体は本文中では用いません。ご了承ください。
あまり格好のよいことではないのだけど、内省のために書いた文章がそこそこの分量に達したし、誰かにとって有益になるかもしれないと思ったので自分語りをしてみます。軽い話ではないので、そういうのが好みではない人は読み飛ばしてください。
元々は自分向けの文章なので、表現が雑な箇所もあります。また、敬体は本文中では用いません。ご了承ください。
前回記事では行列ガンマ分布とその規格化定数の計算方法を紹介しました.
wakabame.hatenablog.com
今回はその応用として, 多次元のガンマ分布を観測モデルとするベイズ推論を行いましょう.
次の行列ガンマ分布は, パラメータである正実数 と次の対称行列 に対して
\begin{align*}
\textrm{MGam}(\Lambda|a,B)
&= C(a,B) |\Lambda|^{a-1}\exp\{-\langle B, \Lambda\rangle\},\\
C(a,B)^{-1}
&= |B|^{-(a+\frac{d-1}{2})}\pi^\frac{d(d-1)}{4}\prod_{k=1}^d \Gamma\left(a+\frac{d-k}{2}\right)
\end{align*}
と表せたのでした.
多次元ガウス分布の精度行列のベイズ推論を行う際にあらわれる行列ガンマ分布:
\begin{align*}
\textrm{MGam}(\Lambda|a,B)
&= C(a,B) |\Lambda|^{a-1}\exp\{-\langle B, \Lambda\rangle\}
\end{align*}
の規格化定数
\begin{align*}
C(a,B)^{-1}
&= \int_{\Lambda\in \mathrm{Sym}_+(d)} |\Lambda|^{a-1}\exp\{-\langle B, \Lambda\rangle\}\,d\Lambda
\end{align*}
の計算方法について述べます.
一次元ガウス分布の精度行列のベイズ推論をする際には, 事前分布としてガンマ分布が表れました.
事前分布は正のパラメータである精度の分布であったから, 予測分布(あるいは規格化定数)を求めるために での積分を行いました.
一方, 多次元においては精度行列は一般に正定値の対称行列となり, パラメータの自由度が大きくなります.
そこで, 正定値の対称行列の値域による積分を導入し, 行列ガンマ分布を定義しましょう.
前回記事では, 一次元ガウス分布を観測モデルに持つようなベイズ推論を行うために, 事前分布はパラメータ により,
\begin{align*}
p(\lambda|v )
&= z(v)\exp\left\{v_1\ln{\lambda} -v_2\lambda +v_3\lambda\mu -\dfrac{v_4}{2}\lambda\mu^2 \right\},\\
z(v)^{-1}
&=\sqrt{\dfrac{2\pi}{v_4}}\dfrac{\Gamma(v_1+\frac12)}{v_2^{v_1+\frac12}}
\left(1 -\dfrac{v_3^2}{2v_2v_4}\right)^{-(v_1+\frac12)}
\end{align*}
置いたとき, 事前分布に対する予測分布を
\begin{align*}
p(X_\ast|v)
&= \dfrac{1}{\sqrt{2\pi}}\sqrt{\dfrac{v_4}{v_4+1}}\dfrac{\Gamma(v_1+1)}{\Gamma(v_1+\frac12)}
\left(v_2 -\dfrac{v_3^2}{2v_4}\right)^{-\frac12}\\
&\quad\times \left\{ 1+ \dfrac{v_4}{2(v_4+1)(v_2-\frac{v_3^2}{2v_4})}\left(X_\ast-\frac{v_3}{v_4}\right)^2\right\}^{-(v_1+1)}
\end{align*}
と得ました.
今回記事では, ハイパーパラメータを変数変換することにより, この意味を明確にしていきましょう.
まず, 事前分布は次のようにガンマ分布の確率密度関数 と, 正規分布の確率密度関数 との積に書き換えられます.
\begin{align*}
p(\lambda|v )
&= z(v)\sqrt{\dfrac{2\pi}{v_4}} \lambda^{v_1-\frac12}e^{-(v_2-\frac{v_3^2}{2v_4})\lambda}
\exp\left\{-\dfrac{v_4\lambda}{2}\left(\mu - \frac{v_3}{v_4}\right)^2 \right\}\\
&= \textrm{G}\left(\lambda|v_1+\frac12, v_2-\frac{v_3^2}{2v_4}\right)
N\left(\mu| \frac{v_3}{v_4}, (v_4\lambda)^{-1}\right).
\end{align*}
このことから, 共役事前分布はハイパーパラメータを元にガンマ分布に従う分散を求め,
その分散も用いて正規分布に従う平均を求めるグラフィカルモデルと位置づけることが出来ます.
それでは, 事前分布が
\begin{equation}
p(\lambda|v )=\textrm{G}\left(\lambda|a, b\right)
N\left(\mu| m, (\beta\lambda)^{-1}\right)
\end{equation}
となるようなハイパーパラメータ()に置き換えましょう;
\begin{equation}
\left\{\begin{split}
&a=v_1+\frac12 \\
&b=v_2 - \frac{v_3^2}{2v_4} \\
&m=\frac{v_3}{v_4}\\
&(\beta\lambda)^{-1} = (v_4\lambda)^{-1}
\end{split}\right.
\Rightarrow
\left\{\begin{split}
&v_1= a-\frac12\\
&v_2=b +\frac12 m^2\beta\\
&v_3 = m\beta\\
&v_4= \beta
\end{split}\right.
\end{equation}
すると, 事前分布による予測分布は,
\begin{align*}
p(X_\ast|v)
&=\dfrac{\beta}{2\pi(\beta+1)b}\dfrac{\Gamma(a+\frac12)}{\Gamma(a)}\\
&\quad\times \left\{ 1+\dfrac{\beta}{2(\beta+1)b}(X_\ast-m)^2\right\}^{-(a+\frac12)}
\end{align*}
という, 自由度 の -分布の確率密度関数に近い形をしています.
実はこれは, 自由度 の -分布に従う確率変数 に対して,
\begin{equation}
\tilde{T}=m+\dfrac{\beta a}{(\beta+1)b}T
\end{equation}
により定められる確率変数 の確率密度関数です.
さらに, 事後分布 による予測分布 を求めましょう.
の代わりに を代入すればよいので,
新たにパラメータ を次のように置きます;
\begin{equation}
\left\{\begin{split}
&\hat{a}=(v+\hat{v})_1+\frac12=a+\dfrac{N}{2} \\
&\hat{b}=(v+\hat{v})_2 - \frac{(v+\hat{v})_3^2}{2(v+\hat{v})_4}
=b+\dfrac12\left(\sum X_n^2 +m^2\beta -\hat{m}^2\hat{\beta}\right) \\
&\hat{m}=\frac{(v+\hat{v})_3}{(v+\hat{v})_4}=\dfrac{m\beta+\sum X_n}{\beta+N}\\
&\hat{\beta} = (v+\hat{v})_4=\beta+N
\end{split}\right.
\end{equation}
すると, 事後分布による予測分布は,
\begin{align*}
p(X_\ast|\bar{X}, v)
&=\dfrac{\hat{\beta}}{2\pi(\hat{\beta}+1)\hat{b}}\dfrac{\Gamma(\hat{a}+\frac12)}{\Gamma(\hat{a})}\\
&\quad\times \left\{ 1+\dfrac{\hat{\beta}}{2(\hat{\beta}+1)\hat{b}}(X_\ast-\hat{m})^2\right\}^{-(\hat{a}+\frac12)}
\end{align*}
です.
したがって, 事後分布による予測分布はスチューデントの -分布に従い, 自由度, 中心, スケールパラメータはそれぞれ
\begin{align*}
2\hat{a} &= 2a + N\\
&= N + O(1), \quad \text{as } N\rightarrow \infty,\\
\hat{m} &= \dfrac{m\beta+\sum X_n}{\beta+N} \\
&= \dfrac{1}{N}\sum_{n=1}^N X_n +o(1), \quad \text{as } N\rightarrow \infty,\\
\dfrac{\hat{\beta}}{2(\hat{\beta}+1)\hat{b}}
&=\dfrac{N}{(N+\beta)(2b+\sum X_n^2 +m^2\beta -\hat{m}^2\beta)}\\
&= \dfrac{1}{N} \dfrac{N-1}{\displaystyle \sum_{n=1}^N\left( X_n -\frac{1}{N} \sum_{k=1}^N X_k\right)^2}+o(1), \quad \text{as } N\rightarrow \infty
\end{align*}
です.
, はそれぞれ平均, 分散の一致推定量なので, 事後分布は真の平均, 真の分散を持つような正規分布に漸近することが予想できます.
正規分布についてのベイズ推論を考えましょう. 母数が つになることに気をつけましょう. 平均 , 精度 がともに未知の一次元ガウス分布を考えます.
\begin{equation}
N(x|\mu,\lambda^{-1})
= \sqrt{\dfrac{1}{2\pi}}\exp\left\{ \dfrac12\ln{\lambda} -\dfrac{\lambda}{2}x^2 + \lambda\mu x -\dfrac{1}{2}\lambda\mu^2\right\}
\end{equation}
により,
\begin{align*}
f(x)&=\sqrt{\dfrac{1}{2\pi}}, \\
\phi_1(\mu,\lambda)&= \ln{\lambda},\quad \phi_2(\mu,\lambda) =-\lambda,\quad
\phi_3(\mu,\lambda)= \lambda\mu,\quad \phi_4(\mu,\lambda) =-\frac12\lambda\mu^2,\\
g_1(x) &= \dfrac12, \quad g_2(x)=\dfrac12 x^2,\quad
g_3(x)=x,\quad g_4(x)=1
\end{align*}
なる指数型分布族です.
続きを読むベイズ推論をします. パラメータ をもつような確率変数 に対する独立同分布な族 を考えます. 事後分布を書く際にはハイパーパラメータ をしばしば省略します. すなわち とかきましょう.
実数値確率変数 が未知パラメータ をもつ(条件付き)確率分布 に従うとします(つまり, 観測モデルを と定めます).
このとき, が指数型分布族に属するとは, ある値域の次元が等しい関数 , により,
\begin{align*}
p(x| \theta)
&= f(x)\exp\langle \phi(\theta) , g(x)\rangle ,\\
f(x)
&= \left(\int \exp\langle \phi(\theta) , g(x)\rangle\, d\theta\right)^{-1}
\end{align*}
が成立することをいいます. ガウス分布やポアソン分布などの代表的な分布はこの定義に属し, 具体例については次回また書いていきます. 今回は一般論として, この性質だけからベイズ推論における共役事前分布が得られ, さらにそれを用いた予測分布の公式が得られることを示していきます.