前回記事では行列ガンマ分布とその規格化定数の計算方法を紹介しました.
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*}
\mathcal{W}(\Lambda|\nu, W)
&=C(\nu,W)|\Lambda|^\frac{\nu-d-1}{2}\exp\{-\mathrm{tr}(W^{-1}\Lambda)\},\\
C(\nu,W)
&= 2^{-\frac{\nu d}{2}}|W|^{-\frac{\nu}{2}}\pi^{-\frac{d(d-1)}{4}}
\prod_{k=1}^d\Gamma\left(\dfrac{\nu+1-k}{2}\right)
\end{align*}
が確率密度であることがわかります.
この分布を自由度 のウィシャート分布と呼びましょう.
多次元ガウス分布を用いた数理統計ではしばしば表れるもので, この分布の便利さの一端はこの記事を通して理解できます.
平均・分散が未知のときの多次元ガウス分布に対するベイズ推論をしよう
平均 , 精度行列
がともに未知の
次元ガウス分布を考えましょう.
すなわち, 観測モデルを
\begin{align*}
p(X|\nu,\Lambda)
=N(x|\mu,\Lambda^{-1})
\end{align*}
と定めましょう. これは
\begin{align*}
f(x)&=(2\pi)^{-\frac{d}{2}}, \\
\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\langle\mu,\Lambda\mu\rangle ,\\
g_1(x) &= \dfrac12, \quad g_2(x)=\dfrac12 x\otimes x,\quad
g_3(x)=x,\quad g_4(x)=1
\end{align*}
なる指数型分布族です.
指数型分布族については
wakabame.hatenablog.com
で解説を行いました. ここで導出した公式を用いてベイズ推論をしましょう. 1次元のガウス分布については
【ベイズ推論】一次元ガウス分布の平均・分散についてのベイズ推論の公式 - わかばめにっき
【ベイズ推論】一次元ガウス分布の平均・分散についてのベイズ推論の意味 - わかばめにっき
で行いましたが, 多次元の場合は行列で表している分やや複雑に感じられますが, 前回得た公式を活用すれば以下のように計算を進めることができます.
まず, 事前分布はパラメータ により,
\begin{align*}
p(\mu,\Lambda|v )
&= z(v)\exp\left\{v_1\ln{|\Lambda|} -\langle V_2,\Lambda\rangle
+\langle v_2,\Lambda\mu\rangle -\dfrac{v_4}{2}\langle \mu,\Lambda\mu \right\},\\
z(v)^{-1}
&= \int_{\mathbb{R}^d} \int_{\Lambda\in\mathrm{Sym}_+(d)} \exp\left\{v_1\ln{|\Lambda|} -\langle V_2,\Lambda\rangle
+\langle v_2,\Lambda\mu\rangle -\dfrac{v_4}{2}\langle \mu,\Lambda\mu \rangle\right\} \, d\lambda d\mu \\
&=\int_{\Lambda\in\mathrm{Sym}_+(d)}|\Lambda|^{v_1} e^{-\langle V_2,\Lambda\rangle}
\int_{\mathbb{R}^d}\exp\left\{\langle v_3.\Lambda\mu\rangle -\dfrac{v_4}{2}\langle \mu,\Lambda\mu\rangle \right\} \,d\mu \,d\lambda\\
&=\sqrt{\dfrac{(2\pi)^d}{v_4}} \int_{\Lambda\in\mathrm{Sym}_+(d)}
|\Lambda|^{v_1-\frac12} \exp \left\{-(V_2-\frac{v_3\otimes v_3}{2v_4}),\Lambda\right\} \,d\Lambda\\
&=\sqrt{\dfrac{(2\pi)^d}{v_4}}\left|V_2-\frac{v_3\otimes v_3}{2v_4}\right|^{-(v_1+\frac{d}{2})}
\pi^\frac{d(d-1)}{4}\prod_{k=1}^d \Gamma(v_1+\frac12+\frac{d-k}{2})
\end{align*}
となります. 最後の式変形で前回記事中で証明した公式を用いました.
事前分布は次のように書き換えられます.
\begin{align*}
p(\mu,\Lambda|v )
&= \textrm{MGam}\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)\\
&= \mathcal{W}\left(\Lambda|d+2v_1, 2\left(V_2-\frac{v_3\otimes v_3}{2v_4}\right)^{-1}\right)
N\left(\mu| \frac{v_3}{v_4}, (v_4\lambda)^{-1}\right)
\end{align*}
であることから, 共役事前分布はハイパーパラメータを元に精度行列を求め, その精度行列を用いて平均を求めるグラフィカルモデルと位置づけることができます.
さらに, パラメータを直観的に理解できるようにするために, 事前分布が となるようなハイパーパラメータ(
)置き換えます;
\begin{equation}
\left\{\begin{split}
&\nu=d+2v_1 \\
&W=\frac12 \left(V_2-\frac{v_3\otimes v_3}{2v_4}\right)^{-1} \\
&m=\frac{v_3}{v_4}\\
&(\beta\lambda)^{-1} = (v_4\lambda)^{-1}
\end{split}\right.
\Rightarrow
\left\{\begin{split}
&v_1= \dfrac{\nu-d}{2}\\
&V_2= \frac12 W^{-1}+\frac{\beta}{2}m\otimes m\\
&v_3 = m\beta\\
&v_4= \beta
\end{split}\right.
\end{equation}
事前分布を用いた予測分布は,
\begin{align*}\displaystyle
p(X_\ast|v)
&= \dfrac{f(X_\ast)z(v)}{z(v+g(X_\ast))}\\
&= \dfrac{\dfrac{1}{(2\pi)^\frac{d}{2}}
\sqrt{\dfrac{v_4}{(2\pi)^d}}\left|V_2-\frac{v_3\otimes v_3}{2v_4}\right|^{v_1+\frac{d}{2}}
\pi^{-\frac{d(d-1)}{4}}\displaystyle\prod_{k=1}^d \Gamma(v_1+\frac12+\frac{d-k}{2})^{-1}}
{\sqrt{\dfrac{v_4+1}{(2\pi)^d}}\left|V_2+\frac12 X_\ast\otimes X_\ast
-\frac{(v_3+X_\ast)\otimes (v_3+X_\ast)}{2(v_4+1)}\right|^{v_1+\frac{d+1}{2}}
\pi^{-\frac{d(d-1)}{4}}\displaystyle\prod_{k=1}^d \Gamma(v_1+1+\frac{d-k}{2})^{-1}}\\
&= \dfrac{1}{(2\pi)^\frac{d}{2}}\sqrt{\dfrac{v_4}{v_4+1}}\dfrac{\Gamma(v_1+\frac{d+1}{2})}{\Gamma(v_1+\frac12)}
\left|V_2 -\dfrac{v_3\otimes v_3}{2v_4}\right|^{v_1+\frac{d}{2}}\\
&\quad \times \left| V_2+\frac12 X_\ast\otimes X_\ast
-\dfrac{(v_3+X_\ast)\otimes(v_3+X_\ast)}{2(v_4+1)}\right|^{-(v_1+\frac{d+1}{2})}\\
&= \dfrac{1}{(2\pi)^\frac{d}{2}}\sqrt{\dfrac{\beta}{\beta+1}}\dfrac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu-d+1}{2})}
\left|\dfrac12 W^{-1} +\frac{\beta}{2}m\otimes m -\frac{\beta}{2}m\otimes m \right|^{\frac{\nu}{2}}\\
&\quad \times \left| \dfrac12 W^{-1} +\frac12 X_\ast\otimes X_\ast+\frac{\beta}{2}m\otimes m
-\frac{\beta m\otimes \beta m-2X_\ast\otimes \beta m +X_\ast\otimes X_\ast}{2(\beta+1)}\right|^{-\frac{\nu+1}{2}}
\end{align*}
です. ここで,
\begin{align*}
&\left| \dfrac12 W^{-1} +\frac12 X_\ast\otimes X_\ast+\frac{\beta}{2}m\otimes m
-\frac{\beta m\otimes \beta m-2X_\ast\otimes \beta m X_\ast\otimes X_\ast}{2(\beta+1)}\right|\\
&\quad =
\left| \frac12 W^{-1}+\frac{\beta }{2(\beta+1)}(X_\ast-m)\otimes (X_\ast-m)\right|
\end{align*}
と計算をしておくと,
\begin{align*}
p(X_\ast|v)
&= \sqrt{\dfrac{\beta}{(2\pi)^\frac{d}{2}(\beta+1)}}
\dfrac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu-d+1}{2})}
\dfrac{\left| \frac12 W^{-1}\right|^{\frac{\nu}{2}}}{\left| \frac12 W^{-1}+\frac{\beta }{2(\beta+1)}(X_\ast-m)\otimes (X_\ast-m)\right|^{\frac{\nu+1}{2}}}\\
&= \sqrt{\dfrac{\beta|W|}{\pi^d(\beta+1)}}\dfrac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu-d+1}{2})}
\left| 1+\frac{\beta }{\beta+1}
\langle X_\ast-m, W(X_\ast-m)\rangle\right|^{-\frac{\nu+1}{2}}
\end{align*}
となり, 事前分布に対する予測分布が得られました.
さらに, 事後分布 による予測分布
を求めましょう.
の代わりに
を代入すればよいので,
新たにパラメータ を次のように置きます;
\begin{equation}
\left\{\begin{split}
\hat{\nu}&=\nu+N \\
\hat{W}&=\dfrac12\left(V_2+\frac12\sum X_n\otimes X_n -\dfrac{(v_3+\sum X_n)\otimes(v_3+\sum X_n)}{2(v_4+N)}\right) \\
& = \left( W^{-1} +\sum X_n\otimes X_n +\beta m\otimes m -\hat{\beta} \hat{m}\otimes \hat{m} \right)^{-1}\\
\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)
&=\sqrt{\dfrac{(\beta+N)|\hat{W}|}{\pi^d(\beta+N+1)}}\dfrac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu-d+1}{2})}
\left| 1+\frac{\beta+N }{\beta+N+1}
\langle X_\ast-\hat{m}, \hat{W}(X_\ast-\hat{m})\rangle\right|^{-\frac{\nu+N+1}{2}}
\end{align*}
です.
したがって, 事後分布による予測分布はスチューデントの -分布に従い,
を十分大きくとる(つまり, 観測データを多くしていく)と, 事後分布は真の平均・真の分散を持つような正規分布に漸近することがわかりました.
