からっぽのしょこ

読んだら書く!書いたら読む!読書読読書読書♪同じ事は二度調べ(たく)ない

3.4.1:多次元ガウス分布の学習と予測:平均が未知の場合【緑ベイズ入門のノート】

はじめに

 『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

 この記事は3.4.1項の内容です。尤度関数と事前分布を多次元ガウス分布とした場合の平均パラメータの事後分布を導出し、学習した事後分布を用いた予測分布を導出します。またその推論過程をR言語で実装します。

 省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

3.4.1 平均が未知の場合

・観測モデルの設定

 多次元ガウス分布に従うと仮定する$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$を用いて平均パラメータ$\boldsymbol{\mu}$の事後分布を求めていく。

 観測データ$\mathbf{X}$、未知の平均パラメータ$\boldsymbol{\mu}$、共分散行列$\boldsymbol{\Sigma}$、精度行列(パラメータ)を$\boldsymbol{\Lambda}$をそれぞれ

$$ \mathbf{X} = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1D} \\ x_{21} & x_{22} & \cdots & x_{2D} \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{ND} \end{bmatrix} ,\ \boldsymbol{\mu} = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \cdots \\ \mu_D \end{bmatrix} ,\ \boldsymbol{\Sigma} = \begin{bmatrix} \sigma_{11}^2 & \sigma_{12}^2 & \cdots & \sigma_{1D}^2 \\ \sigma_{21}^2 & \sigma_{22}^2 & \cdots & \sigma_{2D}^2 \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{D1}^2 & \sigma_{D2}^2 & \cdots & \sigma_{DD}^2 \end{bmatrix} ,\ \boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1} $$

としたとき、観測モデルは次のようになる。

$$ p(\mathbf{x} | \boldsymbol{\mu}) = \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1}) \tag{3.96} $$

 また$\boldsymbol{\mu}$の事前分布$p(\boldsymbol{\mu})$を、平均パラメータ$\mathbf{m}$、精度パラメータ$\boldsymbol{\Lambda}_{\boldsymbol{\mu}}$を持つ$D$次元のガウス分布

$$ p(\boldsymbol{\mu}) = \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, \boldsymbol{\Lambda}^{-1}) \tag{3.97} $$

とする。

・事後分布の導出

 観測データ$\mathbf{X}$によって学習した$\boldsymbol{\mu}$の事後分布$p(\boldsymbol{\mu} | \mathbf{X})$は、観測モデルに対してベイズの定理を用いて

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

となる。分母の$p(\mathbf{X})$は$\boldsymbol{\mu}$に影響しないため省略して、比例関係にのみ注目する。省略した部分については、最後に正規化することで対応できる。

 次にこの分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\boldsymbol{\mu}$に関して整理すると

$$ \begin{align} \ln p(\boldsymbol{\mu} | \mathbf{X}) &= \sum_{n=1}^N \ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1}) + \ln \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1}) + \mathrm{const.} \\ &= \sum_{n=1}^N - \frac{1}{2} \Bigl\{ (\mathbf{x}_n - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x}_n - \boldsymbol{\mu}) + \ln |\boldsymbol{\Lambda}^{-1}| + D \ln 2 \pi \Bigr\} \\ &\qquad - \frac{1}{2} \Bigl\{ (\boldsymbol{\mu} - \mathbf{m})^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} (\boldsymbol{\mu} - \mathbf{m}) + \ln |\boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1}| + D \ln 2 \pi \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Biggl\{ \sum_{n=1}^N \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n - \sum_{n=1}^N \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} - \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + N \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} \Biggr.\\ &\qquad \Biggl. + \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \boldsymbol{\mu} - \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} - \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \boldsymbol{\mu} + \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \Biggr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Biggl\{ - 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + N \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} + \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \boldsymbol{\mu} - 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \Biggr\} + \mathrm{const.} \\ &= - \frac{1}{2} \left\{ \boldsymbol{\mu^{\top}} (N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}}) \boldsymbol{\mu} - 2 \boldsymbol{\mu}^{\top} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \right\} + \mathrm{const.} \tag{3.98} \end{align} $$

【途中式の途中式】

  1. 具体的な式に置き換える。
  2. 式(A.1)より、括弧を展開する。また$\boldsymbol{\mu}$に影響しない項を$\mathrm{const.}$にまとめる。
  3. 次の関係から式を整理する。また$\boldsymbol{\mu}$に影響しない項を$\mathrm{const.}$にまとめる。

 前の因子について、精度パラメータを

$$ \boldsymbol{\Lambda} = \begin{bmatrix} \lambda_{11} & \cdots & \lambda_{1D} \\ \vdots & \ddots & \vdots \\ \lambda_{D1} & \cdots & \lambda_{DD} \end{bmatrix} $$

とすると

$$ \begin{aligned} \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} &= \begin{bmatrix} x_{n1} & \cdots & x_{nD} \end{bmatrix} \begin{bmatrix} \lambda_{11} & \cdots & \lambda_{1D} \\ \vdots & \ddots & \vdots \\ \lambda_{D1} & \cdots & \lambda_{DD} \end{bmatrix} \begin{bmatrix} \mu_1 \\ \vdots \\ \mu_D \end{bmatrix} \\ &= \begin{bmatrix} \sum_{d=1}^D x_{nd} \lambda_{d1} & \cdots & \sum_{d=1}^D x_{nd} \lambda_{dD} \end{bmatrix} \begin{bmatrix} \mu_1 \\ \vdots \\ \mu_D \end{bmatrix} \\ &= \sum_{d=1}^D x_{nd} \lambda_{d1} \mu_1 + \cdots + \sum_{d=1}^D x_{nd} \lambda_{dD} \mu_D \\ &= \sum_{d'=1}^D \sum_{d=1}^D x_{nd} \lambda_{dd'} \mu_{d'} \end{aligned} $$

であり、また

$$ \begin{aligned} \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_n &= \begin{bmatrix} \mu_1 & \cdots & \mu_D \end{bmatrix} \begin{bmatrix} \lambda_{11} & \cdots & \lambda_{1D} \\ \vdots & \ddots & \vdots \\ \lambda_{D1} & \cdots & \lambda_{DD} \end{bmatrix} \begin{bmatrix} x_{n1} \\ \vdots \\ x_{nD} \end{bmatrix} \\ &= \begin{bmatrix} \sum_{d=1}^D \mu_d \lambda_{d1} & \cdots & \sum_{d=1}^D \mu_d \lambda_{dD} \end{bmatrix} \begin{bmatrix} x_{n1} \\ \vdots \\ x_{nD} \end{bmatrix} \\ &= \sum_{d=1}^D \mu_d \lambda_{d1} x_{n1} + \cdots + \sum_{d=1}^D \mu_d \lambda_{dD} x_{nD} \\ &= \sum_{d'=1}^D \sum_{d=1}^D \mu_d \lambda_{dd'} x_{nd'} \end{aligned} $$

なので

$$ \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} = \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_n \tag{1} $$

である。

 この関係は次のようにも考えられる。$\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}$の項は、それぞれ$1 \times D$、$D \times D$、$D \times 1$の行列の積なので、スカラーになる。スカラーは転置しても影響しないので、式(A.2)より

$$ \begin{align} \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} &= ( \mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} )^{\top} \\ &= \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}^{\top} \mathbf{x}_n \\ &= \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_n \tag{1} \end{align} $$

となる。また精度行列$\boldsymbol{\Lambda}$は対称行列なので式(A.3)より、転置しても影響しない。

 後の因子についても同様に

$$ \begin{aligned} \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} &= ( \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} )^{\top} \\ &= \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{\top} \boldsymbol{\mu} \\ &= \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \boldsymbol{\mu} \end{aligned} $$

である。

  1. 式を整理する。


となる。

 1次元ガウス分布のときと同様に、事後分布$p(\boldsymbol{\mu} | \mathbf{X})$がガウス分布になることから、次のようにおき

$$ p(\boldsymbol{\mu} | \mathbf{X}) = \mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}) \tag{3.100} $$

この式の対数をとり、$\boldsymbol{\mu}$に関して整理すると

$$ \ln p(\boldsymbol{\mu} | \mathbf{X}) = -\frac{1}{2} \{ \boldsymbol{\mu}^{\top} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} \boldsymbol{\mu} - 2 \boldsymbol{\mu}^{\top} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} \hat{\mathbf{m}} \} + \mathrm{const.} \tag{3.101} $$

となる。

 従って式(3.99)との対応関係から、事後分布の精度パラメータは

$$ \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} = N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \tag{3.102} $$

となり、平均パラメータは

$$ \begin{align} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} \hat{\mathbf{m}} &= \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \tag{3.103} \end{align} $$

と求められる。

・予測分布

 続いて未観測のデータ$\mathbf{x}_{*}$に対する予測分布$p(\mathbf{x}_{*})$を求めていく。

 1次元ガウス分布のときと同様に、パラメータ$\boldsymbol{\mu}$を周辺化(積分計算)して予測分布

$$ p(\mathbf{x}_{*}) = \int p(\mathbf{x}_{*} | \boldsymbol{\mu}) p(\boldsymbol{\mu}) d\boldsymbol{\mu} $$

を求めることは避け、ベイズの定理を用いて

$$ p(\boldsymbol{\mu} | \mathbf{x}_{*}) = \frac{ p(\mathbf{x}_{*} | \boldsymbol{\mu}) p(\boldsymbol{\mu}) }{ p(\mathbf{x}_{*}) } $$

予測分布を求める。

 この式の両辺の対数をとり、$p(\mathbf{x}_{*})$に関して式を整理すると

$$ \begin{align} \ln p(\boldsymbol{\mu} | \mathbf{x}_{*}) &= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu}) - \ln p(\mathbf{x}_{*}) + \mathrm{const.} \\ \ln p(\mathbf{x}_{*}) &= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu}) - \ln p(\boldsymbol{\mu} | \mathbf{x}_{*}) + \mathrm{const.} \tag{3.104} \end{align} $$

で計算できることが分かる。

 2つ目の項は、データが1つ($N= 1$)の事後分布(3.100)と捉えられることから、式(3.102),(3.103)よりパラメータを

$$ \begin{align} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} &= \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \\ \mathbf{m}_{\mathbf{x}_{*}} &= \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \\ &= ( \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \tag{3.106} \end{align} $$

とおくと

$$ p(\boldsymbol{\mu} | \mathbf{x}_{*}) = \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}_{\mathbf{x}_{*}}, \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1}) \tag{3.105} $$

となる。

 この式を式(3.104)に代入して、$\mathbf{x}$に関して整理すると

$$ \begin{align} \ln p(\mathbf{x}_{*}) &= \ln \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1}) - \ln \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}_{\mathbf{x}_{*}}, \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1}) + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ (\mathbf{x}_{*} - \boldsymbol{\mu})^{\top} \boldsymbol{\Lambda} (\mathbf{x}_{*} - \boldsymbol{\mu}) + \ln |\boldsymbol{\Lambda}^{-1}| + D \ln 2 \pi \Bigr\} \\ &\qquad + \frac{1}{2} \Bigl\{ (\boldsymbol{\mu} - \mathbf{m}_{\mathbf{x}_{*}})^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} (\boldsymbol{\mu} - \mathbf{m}_{\mathbf{x}_{*}}) + \ln |\boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1}| + D \ln 2 \pi \Bigr\} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} - \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} \Bigr.\\ &\qquad \Bigl. - \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\mu} + \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \mathbf{m}_{\mathbf{x}_{*}} + \mathbf{m}_{\mathbf{x}_{*}}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\mu} - \mathbf{m}_{\mathbf{x}_{*}}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \mathbf{m}_{\mathbf{x}_{*}} \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - 2 \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} \Bigr.\\ &\qquad \Bigl. + 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \mathbf{m}_{\mathbf{x}_{*}} - \mathbf{m}_{\mathbf{x}_{*}}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \mathbf{m}_{\mathbf{x}_{*}} \Bigr\} + \mathrm{const.} \end{align} $$

【途中式の途中式】

  1. 具体的な式に置き換える。
  2. 式(A.1)より、括弧を展開する。
  3. 式(1)のときと同様に$\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu} = \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*}$、$\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \mathbf{m}_{\mathbf{x}_{*}} = \mathbf{m}_{\mathbf{x}_{*}}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\mu}$より、項をまとめる。


となる。適宜$\mathbf{x}_{*}$に影響しない項を$\mathrm{const.}$にまとめている。
 更に$\mathbf{m}_{\mathbf{x}_{*}}$に式(3.106)を代入すると

$$ \begin{align} \ln p(\mathbf{x}_{*}) &= - \frac{1}{2} \Bigl[ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} \Bigr.\\ &\qquad \Bigl. + 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) - \Bigl\{ \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Bigr\}^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Big] + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} \Bigr.\\ &\qquad \Bigl. + 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} + 2 \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} - ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} )^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - (\boldsymbol{\Lambda} \mathbf{x}_{*})^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda} \mathbf{x}_{*} - (\boldsymbol{\Lambda} \mathbf{x}_{*})^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} - (\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m})^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda} \mathbf{x}_{*} - (\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m})^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*} - \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda} \mathbf{x}_{*} - \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} - \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda} \mathbf{x}_{*} \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl\{ \mathbf{x}_{*}^{\top} \Bigl( \boldsymbol{\Lambda} - \boldsymbol{\Lambda} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda} \Bigr) \mathbf{x}_{*} - 2 \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \Bigr\} + \mathrm{const.} \\ &= - \frac{1}{2} \Bigl[ \mathbf{x}_{*}^{\top} \Bigl\{ \boldsymbol{\Lambda} - \boldsymbol{\Lambda} ( \boldsymbol{\Lambda} - \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} \boldsymbol{\Lambda} \Bigr\} \mathbf{x}_{*} - 2 \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} ( \boldsymbol{\Lambda} - \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \Bigr] + \mathrm{const.} \tag{3.107} \end{align} $$

【途中式の途中式】

  1. 後の因子について波括弧の項全体を$\mathbf{A}$、丸括弧の項全体を$\mathbf{B}$として考える。また式(A.4)より、$\boldsymbol{\Lambda}_{\mathbf{x}_{*}} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} = \mathbf{I}$である。

 $\mathbf{A}^{\top}$は$D$次元の横ベクトル、$\mathbf{B}$は$D$次元の縦ベクトルになるので、計算結果はスカラーとなる。よって式(1)のときと同様に、転置しても影響しない。
 更に式(A.2)より、$\mathbf{A}^{\top} \mathbf{B} = (\mathbf{A}^{\top} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}$となることから

$$ \begin{aligned} &\Bigl\{ \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Bigr\}^{\top} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \\ &= \Bigl[ \Bigl\{ \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Bigr\}^{\top} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \Bigr]^{\top} \\ &= ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} )^{\top} \boldsymbol{\Lambda}_{\mathbf{x}_{*}}^{-1} ( \boldsymbol{\Lambda} \mathbf{x}_{*} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} ) \end{aligned} $$

となる。

  1. 式(A.1)より、括弧を展開する。
  2. 式(A.2)より、括弧を展開する。また精度行列$\boldsymbol{\Lambda},\ \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$は対称行列なので式(A.3)より、転置できる。
    • よって、$(\boldsymbol{\Lambda} \mathbf{x}_{*})^{\top} = \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda}^{\top} = \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda}$と変形できる。
    • 同様に、$(\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m})^{\top} = \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{\top} = \mathbf{m}^{\top} \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$となる。
  3. 式を整理する。
  4. $\boldsymbol{\Lambda}_{\mathbf{x}_{*}} = \boldsymbol{\Lambda} - \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$を代入する。


となる。

 事後分布のときと同様に、この式を平均$\boldsymbol{\mu}_{*}$、精度$\boldsymbol{\Lambda}_{*}^{-1}$の多次元ガウス分布

$$ p(\mathbf{x}_{*}) = \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}_{*}, \boldsymbol{\Lambda}_{*}^{-1}) \tag{3.107} $$

とおき、対数をとって$\mathbf{x}_{*}$に関して整理すると

$$ \ln p(\mathbf{x}_{*}) = -\frac{1}{2} \{ \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda}_{*} \mathbf{x}_{*} - 2 \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda}_{*} \boldsymbol{\mu}_{*} \} + \mathrm{const.} \tag{3.101} $$

となる。

 従って式(3.107)と対応づけると、精度と平均は

$$ \begin{align} \boldsymbol{\Lambda}_{*} &= \boldsymbol{\Lambda} - \boldsymbol{\Lambda} ( \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} \boldsymbol{\Lambda} \\ &= ( \boldsymbol{\Lambda}^{-1} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1} )^{-1} \tag{3.109}\\ \boldsymbol{\mu}_{*} &= \boldsymbol{\Lambda}_{*}^{-1} \boldsymbol{\Lambda} ( \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} )^{-1} \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \\ &= \boldsymbol{\Lambda}_{*}^{-1} ( \boldsymbol{\Lambda}^{-1} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1} )^{-1} \mathbf{m} \\ &= \boldsymbol{\Lambda}_{*}^{-1} \boldsymbol{\Lambda}_{*} \mathbf{m} \\ &= \mathbf{m} \tag{3.110} \end{align} $$

【途中式の途中式】

  • $\boldsymbol{\Lambda}_{*}$について

 ウッドベリーの公式(A.7)の$\mathbf{U},\ \mathbf{V}$を単位行列$\mathbf{I}$として考えると

$$ \begin{aligned} (\mathbf{A} + \mathbf{U} \mathbf{B} \mathbf{V})^{-1} &= \mathbf{A}^{-1} - \mathbf{A}^{-1} \mathbf{U} ( \mathbf{B}^{-1} + \mathbf{V} \mathbf{A}^{-1} \mathbf{U} )^{-1} \mathbf{V} \mathbf{A}^{-1} \\ (\mathbf{A} + \mathbf{B})^{-1} &= \mathbf{A}^{-1} - \mathbf{A}^{-1} ( \mathbf{B}^{-1} + \mathbf{A}^{-1} )^{-1} \mathbf{A}^{-1} \end{aligned} $$

となる。
 $\boldsymbol{\Lambda}_{*}$の計算式がこの式の右辺の形をしていることから、$\mathbf{A} = \boldsymbol{\Lambda}^{-1},\ \mathbf{B} = \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1}$として左辺の形に変形する。

  • $\boldsymbol{\mu}_{*}$について

 「The Matrix Cookbook」の式(163)より

$$ (\mathbf{A}^{-1} + \mathbf{B}^{-1})^{-1} = \mathbf{A} (\mathbf{A} + \mathbf{B})^{-1} \mathbf{B} $$

この式の右辺の形から左辺の形に変形する。また式(3.109)より置き換える。



が求まる。

 事前分布パラメータ$\boldsymbol{\lambda}_{\boldsymbol{\mu}},\ \mathbf{m}$を事後分布のパラメータ$\hat{\boldsymbol{\lambda}}_{\boldsymbol{\mu}},\ \hat{\mathbf{m}}$に置き換えると、観測データによって学習を行った予測分布$p(\mathbf{x}_{*} | \mathbf{X})$のパラメータ

$$ \begin{aligned} \hat{\boldsymbol{\Lambda}}_{*} &= ( \boldsymbol{\Lambda}^{-1} + \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} )^{-1} \\ \hat{\boldsymbol{\mu}}_{*} &= \hat{\mathbf{m}} \end{aligned} $$

が得られる。

・Rでやってみよう

 多次元ガウス分布に従いランダムに生成したデータを用いて、パラメータを推定してみましょう。

 利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(mvtnorm)

 mvtnormパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvnorm()と確率密度計算関数dmvnorm()を使います。

 観測モデルのパラメータを指定します。この例では2次元のグラフで表現するため、$D = 2$のときのみ動作します。

# データ数を指定
N <- 100

# 観測モデルのパラメータを指定
mu_truth_d <- c(25, 50)
sigma_dd <- matrix(c(50, 30, 30, 50), nrow = 2, ncol = 2)
lambda_dd <- solve(sigma_dd^2)

 観測モデルの平均パラメータ$\boldsymbol{\mu}$をmu_truthとします。この項の例では未知の値であり、この値を推論するのが目的になります。
 $\boldsymbol{\Sigma}$ではなく、$(\sigma_{11}, \cdots, \sigma_{DD})$をsigma_ddとします。これは作図時に標準偏差を利用するためです。matrix()のデフォルトの仕様上、$(\sigma_{11}, \sigma_{12}, \sigma_{21}, \sigma_{22})$の順番に値を指定します。
 sigma_ddを使って精度パラメータ$\boldsymbol{\Lambda}$を計算します。solve()は逆行列を求める関数です。

 同様に事前分布のパラメータを指定します。

# 事前分布のパラメータを指定
m_d <- c(0, 0)
sigma_mu_dd <- matrix(c(100, 0, 0, 100), nrow = 2, ncol = 2)
lambda_mu_dd <- solve(sigma_mu_dd^2)

 事前分布の平均パラメータ$\mathbf{m}$をm_d、精度パラメータ$\boldsymbol{\Lambda}_{\boldsymbol{\mu}}$をlambda_mu_ddとします。

 gif画像作成用のコードでは、事後分布の初期値とも言えます。

 作図時に利用する格子状の点(データ)と真の平均値の点のデータフレームを作成しておきます。

# 作図用の点を生成
x_vec <- seq(mu_truth_d[1] - 2 * sigma_dd[1, 1], mu_truth_d[1] + 2 * sigma_dd[1, 1], by = 0.25)
y_vec <- seq(mu_truth_d[2] - 2 * sigma_dd[2, 2], mu_truth_d[2] + 2 * sigma_dd[2, 2], by = 0.25)
point_df <- tibble(
  x = rep(x_vec, times = length(y_vec)), 
  y = rep(y_vec, each = length(x_vec))
)
mu_df <- tibble(
  x = mu_truth_d[1], 
  y = mu_truth_d[2]
)

 平均から標準偏差の2倍の範囲をグラフ化することにします。

 観測モデルのパラメータに従って乱数を生成します。

# 2次元ガウス分布に従うデータを生成
x_nd <- mvtnorm::rmvnorm(n = N, mean = mu_truth_d, sigma = solve(lambda_dd))
summary(x_nd)
##        V1                V2        
##  Min.   :-137.83   Min.   :-59.34  
##  1st Qu.:  -3.13   1st Qu.: 10.83  
##  Median :  27.88   Median : 49.42  
##  Mean   :  24.28   Mean   : 46.81  
##  3rd Qu.:  52.11   3rd Qu.: 77.88  
##  Max.   : 133.55   Max.   :162.11

 solve()を使って、分散$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}^{-1}$に戻してから、引数sigmaに指定します。

 生成した値x_ndを観測データ$\mathbf{X}$として扱います。

 観測データを散布図で確認しましょう。

# 観測データの散布図を作成
tibble(
  x = x_nd[, 1], 
  y = x_nd[, 2]
) %>% 
  ggplot(data = ., aes(x = x, y = y)) + 
    geom_point() + 
    geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + 
    labs(title = "Multivariate Gaussian Distribution", 
         subtitle = paste0(
           "N=", N, ", mu=(", paste(round(mu_truth_d, 1), collapse = ", "), ")", 
           ", sigma=(", paste(round(sigma_dd, 1), collapse = ", "), ")"
         ), 
         x = expression(x[1]), y = expression(x[2]))

f:id:anemptyarchive:20200911230037p:plain
観測データの散布図


 では事後分布を求めます。

・事後分布

 観測データを使って、事後分布のパラメータを計算します。

# 事後分布のパラメータを計算
lambda_mu_hat_dd <- N * lambda_dd + lambda_mu_dd
m_hat_d <- solve(lambda_mu_hat_dd) %*% (lambda_dd %*% colSums(x_nd) + lambda_mu_dd %*% m_d) %>% 
  as.vector()

 事後分布の精度$\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}$をlambda_mu_hat_dd、平均$\hat{\mathbf{m}}$をm_hat_dとします。計算式(更新式)はそれぞれ次の式です。

$$ \begin{align} \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} &= N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \tag{3.102}\\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \left( \boldsymbol{\Lambda} \sum_{n=1}^N \mathbf{x}_n + \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m} \right) \tag{3.103} \end{align} $$


 求めたパラメータを使って、事後分布(の確率密度)を計算します。

# 事後分布を計算
posterior_df <- tibble(
  xy = point_df, 
  density = mvtnorm::dmvnorm(x = xy, mean = m_hat_d, sigma = solve(lambda_mu_hat_dd)) # 確率密度
) %>% 
  dplyr::select(density) %>% 
  cbind(point_df, .)
head(posterior_df)
##        x   y       density
## 1 -75.00 -50 7.674042e-126
## 2 -74.75 -50 1.609493e-125
## 3 -74.50 -50 3.365922e-125
## 4 -74.25 -50 7.018895e-125
## 5 -74.00 -50 1.459430e-124
## 6 -73.75 -50 3.025854e-124

 作図用にデータフレームとしておきます。

 事後分布を作図します。

# 作図
ggplot() + 
  geom_contour(data = posterior_df, aes(x, y, z = density, color = ..level..)) + # 等高線
  geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) + # 真の平均値
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0(
         "N=", N, ", m_hat=(", paste(round(m_hat_d, 1), collapse = ", "), ")", 
         ", sigma_mu_hat=(", paste(round(sqrt(solve(lambda_mu_hat_dd)), 1), collapse = ", "), ")"
       ), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density") # ラベル

 geom_contour()で、等高線グラフを描画します。格子状の点を渡す必要があります。
 また真の平均値の点も描画しておきます。

f:id:anemptyarchive:20200911230126p:plain
平均$\boldsymbol{\mu}$の事後分布


 続いて予測分布を求めます。

・予測分布

 事後分布パラメータを使って、予測分布のパラメータを計算します。

# 予測分布のパラメータを計算
lambda_star_hat_dd <- solve(solve(lambda_dd) + solve(lambda_mu_hat_dd))
mu_star_hat_d <- m_hat_d

 予測分布の精度$\hat{\boldsymbol{\Lambda}}_{*}$をlambda_star_hat_dd、平均$\hat{\mathbf{m}}$をm_star_hat_dとします。計算式(更新式)はそれぞれ次の式です。

$$ \begin{aligned} \hat{\boldsymbol{\Lambda}}_{*} &= ( \boldsymbol{\Lambda}^{-1} + \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} )^{-1} \\ \hat{\boldsymbol{\mu}}_{*} &= \hat{\mathbf{m}} \end{aligned} $$


 求めたパラメータを使って、予測分布(の確率密度)を計算します。

# 予測分布を計算
predict_df <- tibble(
  xy = point_df, 
  density = mvtnorm::dmvnorm(x = xy, mean = mu_star_hat_d, sigma = solve(lambda_star_hat_dd)) # 確率密度
) %>% 
  dplyr::select(density) %>% 
  cbind(point_df, .)
head(predict_df)
##        x   y      density
## 1 -75.00 -50 4.134909e-06
## 2 -74.75 -50 4.165242e-06
## 3 -74.50 -50 4.195678e-06
## 4 -74.25 -50 4.226216e-06
## 5 -74.00 -50 4.256855e-06
## 6 -73.75 -50 4.287595e-06

 作図用にデータフレームにしておきます。

 予測分布を作図します。

# 作図
ggplot() + 
  geom_contour(data = predict_df, aes(x, y, z = density, color = ..level..)) +  # 等高線
  geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) +  # 真の平均値
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = paste0(
         "N=", N, ", mu_star_hat=(", paste(round(mu_star_hat_d, 1), collapse = ", "), ")", 
         ", sigma_star_hat=(", paste(round(sqrt(solve(lambda_star_hat_dd)), 1), collapse = ", "), ")"
       ), 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density") # ラベル

 事後分布と同様に作図できます。

f:id:anemptyarchive:20200911230240p:plain
未知のデータ$\mathbf{x}_{*}$の予測分布


・おまけ

 gganimateパッケージを利用して、パラメータの推定値の推移を確認するためのgif画像を作成するコードです。

・コード(クリックで展開)

# 追加パッケージ
library(gganimate)

# データ数を指定
N <- 1

# 観測モデルのパラメータを指定
mu_truth_d <- c(25, 50)
sigma_dd <- matrix(c(50, 30, 30, 50), nrow = 2, ncol = 2)
lambda_dd <- solve(sigma_dd^2)

# 事前分布のパラメータを指定
m_d <- c(0, 0)
sigma_mu_dd <- matrix(c(100, 0, 0, 100), nrow = 2, ncol = 2)
lambda_mu_dd <- solve(sigma_mu_dd^2)

# 作図用の点を生成
x_vec <- seq(mu_truth_d[1] - 2 * sigma_dd[1, 1], mu_truth_d[1] + 2 * sigma_dd[1, 1], by = 0.25)
y_vec <- seq(mu_truth_d[2] - 2 * sigma_dd[2, 2], mu_truth_d[2] + 2 * sigma_dd[2, 2], by = 0.25)
point_df <- tibble(
  x = rep(x_vec, times = length(y_vec)), 
  y = rep(y_vec, each = length(x_vec))
)
mu_df <- tibble(
  x = mu_truth_d[1], 
  y = mu_truth_d[2]
)

# 事前分布を計算
posterior_df <- tidyr::tibble(
  xy = point_df, 
  density = mvtnorm::dmvnorm(x = xy, mean = m_d, sigma = solve(lambda_mu_dd)), # 確率密度
  iteration = 0 # 試行回数
) %>% 
  dplyr::select(density, iteration) %>% 
  cbind(point_df, .)


# 予測分布のパラメータを計算
lambda_star_dd <- solve(solve(lambda_dd) + solve(lambda_mu_dd))
mu_star_d <- m_d

# 初期値による予測分布を計算
predict_df <- tibble(
  xy = point_df, 
  density = mvtnorm::dmvnorm(x = xy, mean = mu_star_d, sigma = solve(lambda_star_dd)), # 確率密度
  iteration = 0 # 試行回数
) %>% 
  dplyr::select(density, iteration) %>% 
  cbind(point_df, .)

# 試行回数を指定
max_iter <- 50

# ベイズ推論
for(i in 1:max_iter) {
  
  # 2次元ガウス分布に従うデータを生成
  x_nd <- mvtnorm::rmvnorm(n = N, mean = mu_truth_d, sigma = solve(lambda_dd))
  
  
  # 事後分布のパラメータを更新
  old_lambda_mu_dd <- lambda_mu_dd
  lambda_mu_dd <- N * lambda_dd + lambda_mu_dd
  m_d <- solve(lambda_mu_dd) %*% (lambda_dd %*% colSums(x_nd) + old_lambda_mu_dd %*% m_d) %>% 
    as.vector()
  
  # 事後分布を計算
  tmp_posterior_df <- tidyr::tibble(
    xy = point_df, 
    density = mvtnorm::dmvnorm(x = xy, mean = m_d, sigma = solve(lambda_mu_dd)), # 確率密度
    iteration = i # 試行回数
  ) %>% 
    dplyr::select(density, iteration) %>% 
    cbind(point_df, .)
  
  
  # 予測分布のパラメータを更新
  lambda_star_dd <- solve(solve(lambda_dd) + solve(lambda_mu_dd))
  mu_star_d <- m_d
  
  # 予測分布を計算
  tmp_predict_df <- tibble(
    xy = point_df, 
    density = mvtnorm::dmvnorm(x = xy, mean = mu_star_d, sigma = solve(lambda_star_dd)), # 確率密度
    iteration = i # 試行回数
  ) %>% 
    dplyr::select(density, iteration) %>% 
    cbind(point_df, .)
  
  # 推論結果を結合
  posterior_df <- rbind(posterior_df, tmp_posterior_df)
  predict_df <- rbind(predict_df, tmp_predict_df)
  
  # 動作確認
  print(i)
}

# 事後分布を作図
posterior_graph <- ggplot() + 
  geom_contour(data = posterior_df, aes(x, y, z = density, color = ..level..)) +  # 等高線
  geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) +  # 平均の点
  transition_manual(iteration) +  # フレーム
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = "i={current_frame}", 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density") # ラベル

# gif画像を作成
animate(posterior_graph, nframes = max_iter + 1, fps = 10)

# 予測分布を作図
predict_graph <- ggplot() + 
  geom_contour(data = predict_df, aes(x, y, z = density, color = ..level..)) +  # 等高線
  geom_point(data = mu_df, aes(x = x, y = y), color = "red", shape = 3, size = 5) +  # 平均の点
  transition_manual(iteration) +  # フレーム
  labs(title = "Multivariate Gaussian Distribution", 
       subtitle = "i={current_frame}", 
       x = expression(x[1]), y = expression(x[2]), 
       color = "density") # ラベル

# gif画像を作成
animate(predict_graph, nframes = max_iter + 1, fps = 10)


f:id:anemptyarchive:20200911230357g:plain
平均$\boldsymbol{\mu}$の事後分布の推移

f:id:anemptyarchive:20200911230447g:plain
未知のデータ$\mathbf{x}_{*}$の予測分布の推移


参考文献

 最後の最後でちょろっとカンニングしちゃいました。Julia使いの方はこちらのブログを是非どうぞ。

  • K.B.Petersen, M.S.Pedersen. The Matrix Cookbook. Technical University of Denmark, 2012.

おわりに

 行列計算を見ても発作が起きなくなったので、飛ばしていた多次元ガウス分布の推論に取り組み始めました。そして本当に解けた!本当に成長したなぁと自画自賛しています。
 とはいえ、行列計算時の転置と逆行列の扱いにかなり難儀しました。何だか凄くアクロバットな変形をしましたね。

 近頃はPythonでDL入門してますが、MLSシリーズの新刊も出たことだし、年内にはこの本をやり切ってML入門を修了したいなぁ思ってますがはてさて。

 2020年9月11日は、ハロー!プロジェクトのグループ「Juice=Juice」のメジャーデビュー7周年記念日!

 おめでとうございます!!!!!!!!!!

【次節の内容】

www.anarchive-beta.com