多次元ガウス分布の平均と共分散の事後分布

本記事では共分散行列 ${\bf\Sigma}$ の代わりにその逆行列である精度行列 ${\boldsymbol\Lambda}={\bf \Sigma}^{-1}$ を用います。

多次元ガウス分布

多次元ガウス分布は以下のように表されるのでした。

$\begin{eqnarray} {\mathcal N}({\bf x}|{\boldsymbol\mu},{\bf\Sigma})=\frac{1}{\sqrt{(2\pi)^D|{\bf\Sigma}|}}\exp\left(-\frac{1}{2}({\bf x}-{\boldsymbol\mu})^\top{\bf\Sigma}^{-1}({\bf x}-{\boldsymbol\mu})\right)\tag{1} \end{eqnarray}$

${\bf x}\in\mathbb{R}^D,\ {\boldsymbol\mu}\in\mathbb{R}^D,\ {\bf\Sigma}\in\mathbb{R}^{D\times D}$ です。
(1)を精度行列 ${\boldsymbol\Lambda}={\bf \Sigma}^{-1}$ を用いて表すと以下のようになります。

$\begin{eqnarray} {\mathcal N}({\bf x}|{\boldsymbol\mu},{\bf\Lambda}^{-1})=\frac{\sqrt{|{\bf\Lambda}|}}{\sqrt{(2\pi)^D}}\exp\left(-\frac{1}{2}({\bf x}-{\boldsymbol\mu})^\top{\bf\Lambda}({\bf x}-{\boldsymbol\mu})\right)\tag{2} \end{eqnarray}$

(2)の式変形では $|{\bf\Sigma}|=\frac{1}{|{\bf\Lambda}|}$ を使用しました。
対数を取ってみます。

$\begin{eqnarray} \ln{\mathcal N}({\bf x}|{\boldsymbol\mu},{\bf\Lambda}^{-1})=-\frac{1}{2}\left(({\bf x}-{\boldsymbol\mu})^\top{\bf\Lambda}({\bf x}-{\boldsymbol\mu})-\ln|{\bf\Lambda}|+D\ln2\pi\right)\tag{3} \end{eqnarray}$

ガウス-ウィシャート分布

ガウス-ウィシャート分布は以下のように表されるのでした。

$\begin{eqnarray} p({\boldsymbol\mu},{\bf\Lambda})&=&{\rm NW}({\boldsymbol\mu},{\bf\Lambda}|{\bf m},\beta,\nu,{\bf W})\\ &=&\mathcal {N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})\mathcal{W}({\bf\Lambda|\nu,{\bf W}})\tag{4} \end{eqnarray}$

${\bf m}\in\mathbb{R}^D,\beta\in\mathbb{R}^+,\nu\in\mathbb{R}, \nu>D-1,{\bf W}\in\mathbb{R}^{D\times D}, {\bf W}>0$ です。

ウィシャート分布

ウィシャート分布は以下のように表されるのでした。

$\begin{eqnarray} {\mathcal W}({\bf\Lambda}|{\nu},{\bf W})=C_{\mathcal W}(\nu,{\bf W})|{\bf\Lambda}|^{\frac{\nu-D-1}{2}}\exp\left(-\frac{1}{2}{\rm Tr}({\bf W}^{-1}{\bf\Lambda})\right)\tag{5} \end{eqnarray}$

${\bf\Lambda}\in\mathbb{R}^{D\times D}, {\bf\Lambda}>0, \nu\in\mathbb{R}, \nu>D-1, {\bf W}\in\mathbb{R}^{D\times D}, {\bf W}>0$ です。
※ ${\bf\Lambda}>0,{\bf W}>0$ は行列 ${\bf\Lambda}$ 、行列 ${\bf W}$ が正定値であることを表しています。

$C_{\mathcal W}(\nu, {\bf W})$ はウィシャート分布の正規化項です。

$\begin{eqnarray} C_{\mathcal W}(\nu, {\bf W})=|{\bf W}|^{-\nu/2}2^{-\nu D/2}\pi^{-D(D-1)/4}\prod_{d=1}^D\Gamma(\frac{\nu+1-d}{2})^{-1}\tag{6}\\ \end{eqnarray}$

対数を取ると以下のようになります。

$\begin{eqnarray} \ln {\mathcal W}({\bf\Lambda}|{\nu},{\bf W})=\frac{\nu-D-1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}{\rm Tr}({\bf W}^{-1}{\bf\Lambda})+\ln C_{\mathcal W}(\nu,{\bf W})\tag{7}\\ \end{eqnarray}$

事後分布

多次元ガウス分布に従う $N$ 個のD次元の連続値データ ${\bf X}=\{{\bf x}_1,\ldots,{\bf x}_N\}$ が与えられた時の平均 ${\boldsymbol\mu}\in\mathbb{R}^D$ と精度行列 ${\bf\Lambda}\in{\mathbb{R}}^{D\times D}$ の事後分布を求めてみます。

f:id:olj611:20210224070232p:plain

${\boldsymbol\mu},{\bf\Lambda}$ の事前分布には多次元ガウス分布の平均と精度の共役事前分布であるガウス-ウィシャート分布を用います。

${\bf m}\in\mathbb{R}^D,\beta\in\mathbb{R}^+,\nu\in\mathbb{R}, \nu>D-1,{\bf W}\in\mathbb{R}^{D\times D}, {\bf W}>0$ は固定されたハイパーパラメータです。

また、 $p({\boldsymbol\mu},{\bf\Lambda})$ は確率の乗法定理を使って以下のように書けます。

$\begin{eqnarray} p({\boldsymbol\mu},{\bf\Lambda})&=&p({\boldsymbol\mu}|{\bf\Lambda})p({\bf\Lambda})\tag{9} \end{eqnarray}$

(8)、(9) より

$\begin{eqnarray} &&p({\boldsymbol\mu}|{\bf\Lambda})=\mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})\tag{10}\\ &&p({\bf\Lambda})=\mathcal{W}({\bf\Lambda}|\nu,{\bf W})\tag{11}\\ \end{eqnarray}$

です。

事後分布 $p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})$ を確率の乗法定理を用いて分解します。

$\begin{eqnarray} p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})&=&p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})p({\bf\Lambda}|{\bf X})\tag{12} \end{eqnarray}$

(12)より $p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})$ と $p({\bf\Lambda}|{\bf X})$ を別々に求めることにします。

まず、 $p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})$ を求めます。ベイズの定理を用いると以下の式になります。

$\begin{eqnarray} p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})&=&\frac{p({\bf X}|{\boldsymbol\mu},{\bf\Lambda})p({\boldsymbol\mu}|{\bf\Lambda})}{p({\bf X}|{\bf\Lambda})}\\ &\propto&p({\bf X}|{\boldsymbol\mu},{\bf\Lambda})p({\boldsymbol\mu}|{\bf\Lambda})\\ &=&\left(\prod_{n=1}^N p({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda})\right)p({\boldsymbol\mu}|{\bf\Lambda})\\ &=&\left(\prod_{n=1}^N \mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})\right)\mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})\tag{13}\\ \end{eqnarray}$

(13)に対数を取って、 ${\boldsymbol\mu}$ についてまとめます。

$\begin{eqnarray} \ln p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})&=&\sum_{n=1}^N \ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})+\ln\mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})+{\rm const.}\\ &=&-\frac{1}{2}\sum_{n=1}^N({\bf x}_n-{\boldsymbol\mu})^\top{\bf\Lambda}({\bf x}_n-{\boldsymbol\mu})-\frac{1}{2}\beta({\boldsymbol\mu}-{\bf m})^\top{\bf\Lambda}({\boldsymbol\mu}-{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\sum_{n=1}^N({\bf x}_n^\top{\bf\Lambda}{\bf x}_n-2{\boldsymbol\mu}^\top{\bf\Lambda}{\bf x}_n+{\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu})-\frac{1}{2}\beta({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu^\top}{\bf\Lambda}{\bf m}+{\bf m}^\top{\bf\Lambda}{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\sum_{n=1}^N(-2{\boldsymbol\mu}^\top{\bf\Lambda}{\bf x}_n+{\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu})-\frac{1}{2}\beta({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu^\top}{\bf\Lambda}{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\left((N+\beta){\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}^\top{\bf\Lambda}\left(\sum_{n=1}^N{\bf x}_n+\beta{\bf m}\right)\right)+{\rm const.}\tag{14}\\ \end{eqnarray}$

(14)より $p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})$ が多次元ガウス分布になっていることが分かります。
$p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})$ の平均と精度を $\hat{\bf m},\hat{\beta}{\bf\Lambda}$ とおき、 ${\boldsymbol\mu}$ について整理します。

$\begin{eqnarray} \ln p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})&=&\ln\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1})\\ &=&-\frac{1}{2}\hat{\beta}({\boldsymbol\mu}-\hat{\bf m})^\top{\bf\Lambda}({\boldsymbol\mu}-\hat{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\hat{\beta}({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}{\bf\Lambda}\hat{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}(\hat{\beta}{\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}{\bf\Lambda}\hat{\beta}\hat{\bf m})+{\rm const.}\tag{15}\\ \end{eqnarray}$

(14)と(15)の係数を比較します。

$\begin{eqnarray} &&\hat{\beta}=N+\beta\tag{16}\\ &&\hat{\bf m}=\frac{1}{\beta}\left(\sum_{n=1}^N{\bf x}_n+\beta{\bf m}\right) \tag{17}\\ \end{eqnarray}$

次に、 $p({\bf\Lambda}|{\bf X})$ を求めます。
その前に $p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})$ にベイズの定理を適用します。

$\begin{eqnarray} p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})&=&\frac{p({\bf X}|{\boldsymbol\mu},{\bf\Lambda})p({\boldsymbol\mu},{\bf\Lambda})}{p({\bf X})}\tag{18} \end{eqnarray}$

(12)より

$\begin{eqnarray} p({\bf\Lambda}|{\bf X})&=&\frac{p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})}{p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})}\tag{19} \end{eqnarray}$

となります。(19)に(18)を代入します。

$\begin{eqnarray} p({\bf\Lambda}|{\bf X})&=&\frac{p({\bf X}|{\boldsymbol\mu},{\bf\Lambda})p({\boldsymbol\mu},{\bf\Lambda})}{p({\bf X})p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})}\\ &\propto&\frac{p({\bf X}|{\boldsymbol\mu},{\bf\Lambda})p({\boldsymbol\mu},{\bf\Lambda})}{p({\boldsymbol\mu}|{\bf\Lambda},{\bf X})}\\ &=&\frac{\left(\prod_{n=1}^N \mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})\right)\mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})\mathcal{W}({\bf\Lambda}|\nu,{\bf W})}{\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1})}\tag{20}\\ \end{eqnarray}$

対数を取ります。

$\begin{eqnarray} \ln p({\bf\Lambda}|{\bf X})&=&\sum_{n=1}^N\ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})+\ln \mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1}) + \ln\mathcal{W}({\bf\Lambda}|\nu,{\bf W})-\ln\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1})+{\rm const.}\tag{21}\\ \end{eqnarray}$

(21)を計算していくのですが、そのまま計算すると式が見にくくなってしまうので各項毎に ${\bf\Lambda}$ についてまとめつつ計算します。
$\sum_{n=1}^N\ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})$ について計算します。

$\begin{eqnarray} \sum_{n=1}^N\ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu},{\bf\Lambda}^{-1})&=&\frac{N}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\sum_{n=1}^N({\bf x}_n-{\boldsymbol\mu})^\top{\bf\Lambda}({\bf x}_n-{\boldsymbol\mu})+{\rm const.}\\ &=&\frac{N}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\left(\sum_{n=1}^N{\bf x}_n^\top{\bf\Lambda}{\bf x}_n-2{\boldsymbol\mu}^\top{\bf\Lambda}\sum_{n=1}^N{\bf x}_n+N{\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}\right)+{\rm const.}\\ &=&\frac{N}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\left({\rm Tr}\left(\sum_{n=1}^N{\bf x}_n{\bf x}_n^\top{\bf\Lambda}\right)-2{\boldsymbol\mu}^\top{\bf\Lambda}\sum_{n=1}^N{\bf x}_n+N{\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}\right)+{\rm const.}\tag{22}\\ \end{eqnarray}$

$\ln \mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1})$ について計算します。

$\begin{eqnarray} \ln \mathcal{N}({\boldsymbol\mu}|{\bf m},(\beta{\bf\Lambda})^{-1}) &=&\frac{1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\beta({\boldsymbol\mu}-{\bf m})^\top{\bf\Lambda}({\boldsymbol\mu}-{\bf m})+{\rm const.}\\ &=&\frac{1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\beta({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}^\top{\bf\Lambda}{\bf m}+{\bf m}^\top{\bf\Lambda}{\bf m})+{\rm const.}\\ &=&\frac{1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}\beta\left({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}^\top{\bf\Lambda}{\bf m}+{\rm{Tr}}({\bf m}{\bf m}^\top{\bf\Lambda})\right)+{\rm const.}\tag{23}\\ \end{eqnarray}$

$\ln\mathcal{W}({\bf\Lambda}|\nu,{\bf W})$ について計算します。

$\begin{eqnarray} \ln\mathcal{W}({\bf\Lambda}|\nu,{\bf W})=\frac{\nu-D-1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}{\rm Tr}({\bf W}^{-1}{\bf\Lambda})+{\rm const.}\tag{24}\\ \end{eqnarray}$

$-\ln\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1})$ について計算します。

$\begin{eqnarray} -\ln\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1})&=&-\frac{1}{2}\ln|{\bf\Lambda}|+\frac{1}{2}\hat{\beta}({\boldsymbol\mu}-\hat{\bf m})^\top{\bf\Lambda}({\boldsymbol\mu}-\hat{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\ln|{\bf\Lambda}|+\frac{1}{2}\hat{\beta}({\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}^\top{\bf\Lambda}\hat{\bf m}+\hat{\bf m}^\top{\bf\Lambda}\hat{\bf m})+{\rm const.}\\ &=&-\frac{1}{2}\ln\lambda+\frac{1}{2}\left((\beta+N){\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2\hat{\beta}{\boldsymbol\mu}^\top{\bf\Lambda}\left(\frac{\sum_{n=1}^N{\bf x}_n+\beta {\bf m}}{\hat{\beta}}\right)+\hat{\beta}\hat{\bf m}^\top{\bf\Lambda}\hat{\bf m}\right)+{\rm const.}\\ &=&-\frac{1}{2}\ln\lambda+\frac{1}{2}\left((\beta+N){\boldsymbol\mu}^\top{\bf\Lambda}{\boldsymbol\mu}-2{\boldsymbol\mu}^\top{\bf\Lambda}\left(\sum_{n=1}^N{\bf x}_n+\beta {\bf m}\right)+\hat{\beta}{\rm Tr}(\hat{\bf m}\hat{\bf m}^\top{\bf\Lambda})\right)+{\rm const.}\tag{25}\\ \end{eqnarray}$

(22)、(23)、(25)の式変形にはトレースの公式を用いました。

$\begin{eqnarray} &&{\bf x}^\top{\bf A}{\bf x}={\rm Tr}({\bf A}{\bf x}{\bf x}^\top)\tag{26}\\ &&{\rm Tr}({\bf A})={\rm Tr}({\bf A}^\top)\tag{27}\\ &&{\rm Tr}({\bf AB})={\rm Tr}({\bf A})+{\rm Tr}({\bf B})\tag{28}\\ \end{eqnarray}$

(22)、(23)、(24)、(25)を(21)に代入します。

$\begin{eqnarray} \ln p({\bf\Lambda}|{\bf X})=\frac{N+\nu-D-1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}{\rm Tr}\left(\left(\sum_{n=1}^N{\bf x}_n{\bf x}_n^\top+\beta{\bf m}{\bf m}^\top-\hat{\beta}\hat{\bf m}\hat{\bf m}^\top+{\bf W}^{-1}\right){\bf\Lambda}\right)+{\rm const.}\tag{29}\\ \end{eqnarray}$

(29)より $p({\bf\Lambda}|{\bf X})$ がウィシャート分布になっていることが分かります。
$p({\bf\Lambda}|{\bf X})$ のパラメータを $\hat{\nu},\hat{\bf W}$ とおき、 ${\bf\Lambda}$ について整理します。

$\begin{eqnarray} \ln p({\bf\Lambda}|{\bf X})&=&\ln {\mathcal W}({\bf\Lambda}|\hat{\nu},\hat{\bf W})\\ &=&\frac{\hat{\nu}-D-1}{2}\ln|{\bf\Lambda}|-\frac{1}{2}{\rm Tr}(\hat{\bf W}^{-1}{\bf\Lambda})+{\rm const.}\tag{30}\\ \end{eqnarray}$

(29)、(30)を係数比較します。

$\begin{eqnarray} &&\hat{\nu}=N+\nu\tag{31}\\ &&\hat{\bf W}^{-1}=\sum_{n=1}^N{\bf x}_n{\bf x}_n^\top+\beta{\bf m}{\bf m}^\top-\hat{\beta}\hat{\bf m}\hat{\bf m}^\top+{\bf W}^{-1}\tag{32}\\ \end{eqnarray}$

(12)、(16)、(17)、(31)、(32)より、 $p(\mu,\lambda|{\bf X})$ がガウス-ウィシャート分布になることが分かりました。

$\begin{eqnarray} &&p({\boldsymbol\mu},{\bf\Lambda}|{\bf X})=\mathcal{N}({\boldsymbol\mu}|\hat{\bf m},(\hat{\beta}{\bf\Lambda})^{-1}){\mathcal W}({\bf\Lambda}|\hat{\nu},\hat{\bf W})\tag{33}\\ &&\hat{\beta}=N+\beta\tag{34}\\ &&\hat{\bf m}=\frac{1}{\beta}\left(\sum_{n=1}^N{\bf x}_n+\beta{\bf m}\right) \tag{35}\\ &&\hat{\nu}=N+\nu\tag{36}\\ &&\hat{\bf W}^{-1}=\sum_{n=1}^N{\bf x}_n{\bf x}_n^\top+\beta{\bf m}{\bf m}^\top-\hat{\beta}\hat{\bf m}\hat{\bf m}^\top+{\bf W}^{-1}\tag{37}\\ \end{eqnarray}$

偉人の名言

f:id:olj611:20210224171647p:plain
天才？そんなものは決してない。ただ勉強です。方法です。不断に計画しているということです。
ロダン

参考文献

ベイズ推論による機械学習入門

機械学習基礎理論独習

誤りがあればご指摘いただけると幸いです。数式が整うまで少し時間かかります。リンクフリーです。