はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」によって、「数式」と「プログラム」から理解するのが目標です。
省略してある内容等ありますので、本とあわせて読んでください。
この記事では、「尤度関数を平均が未知の多次元ガウス分布(多変量正規分布)」、「事前分布を多次元ガウス分布」とした場合の「平均パラメータの事後分布」と「未観測値の予測分布」を導出します。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【前の節の内容】
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
3.4.1 平均が未知の場合
尤度関数を多次元ガウス分布(Maltivariate Gaussian Distribution)・多変量正規分布(Maltivariate Normal Distribution)、事前分布を多次元ガウス分布とするモデルに対するベイズ推論を導出する。
多次元ガウス分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照のこと。
尤度関数の確認
$D$次元ガウス分布に従うと仮定する$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$、$\mathbf{x}_n = (x_{n,1}, x_{n,2}, \cdots, x_{n,D})^{\top}$を用いて、平均パラメータ$\boldsymbol{\mu}$の事後分布と未観測データ$\mathbf{x}_{*}$の予測分布を求めていく。
$\mathbf{X}$の各データが独立に発生しているとの仮定の下で、尤度$p(\mathbf{X} | \mu)$は
$$
\begin{aligned}
p(\mathbf{X} | \boldsymbol{\mu})
&= p(\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N | \boldsymbol{\mu})
\\
&= p(\mathbf{x}_1 | \boldsymbol{\mu})
p(\mathbf{x}_2 | \boldsymbol{\mu})
\cdots
p(\mathbf{x}_N | \boldsymbol{\mu})
\\
&= \prod_{n=1}^N
p(\mathbf{x}_n | \boldsymbol{\mu})
\end{aligned}
$$
と分解できる。
ここで、観測データ$\mathbf{X}$、未知の平均ベクトル(平均パラメータ)$\boldsymbol{\mu}$、精度行列(精度パラメータ)$\boldsymbol{\Lambda}$をそれぞれ
$$
\mathbf{X}
= \begin{bmatrix}
x_{1,1} & x_{1,2} & \cdots & x_{1,D} \\
x_{2,1} & x_{2,2} & \cdots & x_{2,D} \\
\vdots & \vdots & \ddots & \vdots \\
x_{N,1} & x_{N,2} & \cdots & x_{N,D}
\end{bmatrix}
,\
\boldsymbol{\mu}
= \begin{bmatrix}
\mu_1 \\
\mu_2 \\
\vdots \\
\mu_D
\end{bmatrix}
,\
\boldsymbol{\Lambda}
= \begin{bmatrix}
\lambda_{1,1} & \lambda_{1,2} & \cdots & \lambda_{1,D} \\
\lambda_{2,1} & \lambda_{2,2} & \cdots & \lambda_{2,D} \\
\vdots & \vdots & \ddots & \vdots \\
\lambda_{D,1} & \lambda_{D,2} & \cdots & \lambda_{D,D}
\end{bmatrix}
$$
とする。精度行列は、分散共分散行列$\boldsymbol{\Sigma}$の逆行列$\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}$である($\lambda_{i,j}$は固有値のことじゃないよ)。
事後分布の計算
次に、平均パラメータ$\boldsymbol{\mu}$の事後分布$p(\boldsymbol{\mu} | \mathbf{X})$を導出する。
観測データ$\mathbf{X}$が与えられた下での$\boldsymbol{\mu}$の事後分布は、ベイズの定理より
$$
\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.99}
\end{align}
$$
【途中式の途中式】(クリックで展開)
- 尤度と事前分布に、それぞれ具体的な(対数をとった)確率分布の式を代入する。
- $(\mathbf{A} + \mathbf{B})^{\top} = (\mathbf{A}^{\top} + \mathbf{B}^{\top})$より、括弧を展開する。また、$\boldsymbol{\mu}$と無関係な項を$\mathrm{const.}$にまとめる。
- 次の関係から式を整理する。また、$\boldsymbol{\mu}$と無関係な項を$\mathrm{const.}$にまとめる。
1行目の2つ目の項について
$$
\begin{aligned}
\mathbf{x}_n^{\top}
\boldsymbol{\Lambda}
\boldsymbol{\mu}
&= \begin{bmatrix}
x_{n,1} & \cdots & x_{n,D}
\end{bmatrix}
\begin{bmatrix}
\lambda_{1,1} & \cdots & \lambda_{1,D} \\
\vdots & \ddots & \vdots \\
\lambda_{D,1} & \cdots & \lambda_{D,D}
\end{bmatrix}
\begin{bmatrix}
\mu_1 \\
\vdots \\
\mu_D
\end{bmatrix}
\\
&= \begin{bmatrix}
\sum_{d=1}^D x_{n,d} \lambda_{d,1} &
\cdots &
\sum_{d=1}^D x_{n,d} \lambda_{d,D}
\end{bmatrix}
\begin{bmatrix}
\mu_1 \\
\vdots \\
\mu_D
\end{bmatrix}
\\
&= \sum_{d=1}^D x_{n,d} \lambda_{d,1} \mu_1
+ \cdots
+ \sum_{d=1}^D x_{n,d} \lambda_{d,D} \mu_D
\\
&= \sum_{d=1}^D \sum_{d'=1}^D
x_{n,d} \lambda_{d,d'} \mu_{d'}
\end{aligned}
$$
であり、また3つ目の項は
$$
\begin{aligned}
\boldsymbol{\mu}^{\top}
\boldsymbol{\Lambda}
\mathbf{x}_n
&= \begin{bmatrix}
\mu_1 & \cdots & \mu_D
\end{bmatrix}
\begin{bmatrix}
\lambda_{1,1} & \cdots & \lambda_{1,D} \\
\vdots & \ddots & \vdots \\
\lambda_{D,1} & \cdots & \lambda_{D,D}
\end{bmatrix}
\begin{bmatrix}
x_{n,1} \\
\vdots \\
x_{n,D}
\end{bmatrix}
\\
&= \begin{bmatrix}
\sum_{d=1}^D \mu_d \lambda_{d,1} &
\cdots &
\sum_{d=1}^D \mu_d \lambda_{d,D}
\end{bmatrix}
\begin{bmatrix}
x_{n,1} \\
\vdots \\
x_{n,D}
\end{bmatrix}
\\
&= \sum_{d=1}^D \mu_d \lambda_{d,1} x_{n,1}
+ \cdots
+ \sum_{d=1}^D \mu_d \lambda_{d,D} x_{n,D}
\\
&= \sum_{d=1}^D \sum_{d'=1}^D
\mu_d \lambda_{d,d'} x_{n,d'}
\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$の行列なので、その行列の積$\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}$はスカラになる。スカラは転置しても影響しないので、$(\mathbf{A} \mathbf{B} \mathbf{C})^{\top} = \mathbf{C}^{\top} \mathbf{B}^{\top} \mathbf{A}^{\top}$より
$$
\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}$は対称行列なので、転置しても影響せず$\boldsymbol{\Lambda} = \boldsymbol{\Lambda}^{\top}$である。
2行目の項についても同様に
$$
\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}
$$
である。
- $\boldsymbol{\mu}^{\top} \boldsymbol{\mu}$と$\boldsymbol{\mu}^{\top}$の項をそれぞれまとめて式を整理する。
ちなみに$\sum_{n=1}^N \mathbf{x}_n$は次のような$D$次元ベクトルである。
$$
\sum_{n=1}^N \mathbf{x}_n
= \begin{bmatrix}
\sum_{n=1}^N x_{n,1} \\
\vdots \\
\sum_{n=1}^N x_{n,D}
\end{bmatrix}
$$
となる。
事後分布はガウス分布になることから、平均$\hat{\mathbf{m}}$、精度$\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}$の$D$次元ガウス分布
$$
p(\boldsymbol{\mu} | \mathbf{X})
= \mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})
\tag{3.100}
$$
とおく。この式の対数をとり、$\boldsymbol{\mu}$に関して整理すると
$$
\begin{align}
\ln p(\boldsymbol{\mu} | \mathbf{X})
&= \ln \mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})
\\
&= - \frac{1}{2} \Bigl\{
(\boldsymbol{\mu} - \hat{\mathbf{m}})^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
(\boldsymbol{\mu} - \hat{\mathbf{m}})
+ \ln |\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}|
+ D \ln 2 \pi
\Bigr\}
\\
&= - \frac{1}{2} \Bigl\{
\boldsymbol{\mu}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\boldsymbol{\mu}
- \boldsymbol{\mu}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\hat{\mathbf{m}}
- \hat{\mathbf{m}}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\boldsymbol{\mu}
+ \hat{\mathbf{m}}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\hat{\mathbf{m}}
\Biggr\}
+ \mathrm{const.}
\\
&= -\frac{1}{2} \Bigl\{
\boldsymbol{\mu}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\boldsymbol{\mu}
- 2 \boldsymbol{\mu}^{\top}
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\hat{\mathbf{m}}
\Bigr\}
+ \mathrm{const.}
\tag{3.101}
\end{align}
$$
となる。
式(3.99)と式(3.101)との対応関係から、事後分布の精度パラメータは
$$
\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}
\\
\Rightarrow
\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}
$$
と求められる。
式(3.103)、(3.102)が超パラメータ$\hat{\mathbf{m}}, \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}$の計算式(更新式)である。
予測分布の計算
続いて、多次元ガウス分布に従う未観測データ$\mathbf{x}_{*} = (x_{*,1}, x_{*,2}, \cdots, x_{*,D})^{\top}$に対する予測分布を導出する。
事前分布による予測分布
事前分布(観測データによる学習を行っていない分布)$p(\boldsymbol{\mu})$を用いて、未学習の予測分布$p(\mathbf{x}_{*})$を求める。その結果を用いて、学習後の予測分布$p(\mathbf{x}_{*} | \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}_{*})
}
$$
予測分布を求める。
この式の両辺の対数をとり
$$
\ln p(\boldsymbol{\mu} | \mathbf{x}_{*})
= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu})
+ p(\boldsymbol{\mu})
- \ln p(\mathbf{x}_{*})
$$
予測分布に関して式を整理すると
$$
\ln p(\mathbf{x}_{*})
= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu})
- \ln p(\boldsymbol{\mu} | \mathbf{x}_{*})
+ \mathrm{const.}
\tag{3.104}
$$
となる。ただし、$\mathbf{x}_{*}$に影響しない$\ln p(\boldsymbol{\mu})$を$\mathrm{const.}$とおいた。
この式から予測分布の具体的な式を計算する。
$p(\boldsymbol{\mu} | \mathbf{x}_{*})$は、1つのデータ$\mathbf{x}_{*}$が与えられた下での$\boldsymbol{\mu}$の条件付き分布である。つまり$p(\boldsymbol{\mu} | \mathbf{x}_{*})$は、$N$個の観測データが与えられた下での事後分布$p(\boldsymbol{\mu} | \mathbf{X})$と同様の手順で求められる(同様のパラメータになる)。
したがって、事後分布のパラメータ(3.102)、(3.103)を用いると、$p(\boldsymbol{\mu} | x_{*})$は$N = 1$より
$$
\ln p(\boldsymbol{\mu} | \mathbf{x}_{*})
= - \frac{1}{2} \Bigl\{
(\boldsymbol{\mu} - \tilde{\mathbf{m}})^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
(\boldsymbol{\mu} - \tilde{\mathbf{m}})
+ \ln |\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}|
+ D \ln 2 \pi
\Bigr\}
= \ln \mathcal{N}(\boldsymbol{\mu} | \tilde{\mathbf{m}}, \tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})
\tag{3.105}
$$
平均$\tilde{\mathbf{m}}$、精度$\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}$の$D$次元ガウス分布となる。ただし
$$
\begin{align}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
&= \boldsymbol{\Lambda}
+ \boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\\
\tilde{\mathbf{m}}
&= \tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
$$
とおく。
尤度(3.96)と式(3.105)を式(3.104)に代入して、$\mathbf{x}_{*}$に関して式を整理する。
$$
\begin{aligned}
\ln p(\mathbf{x}_{*})
&= \ln \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
- \ln \mathcal{N}(\boldsymbol{\mu} | \tilde{\mathbf{m}}, \tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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} - \tilde{\mathbf{m}})^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
(\boldsymbol{\mu} - \tilde{\mathbf{m}})
+ \ln |\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\boldsymbol{\mu}
+ \boldsymbol{\mu}^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\mathbf{m}}
+ \tilde{\mathbf{m}}^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\boldsymbol{\mu}
- \tilde{\mathbf{m}}^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\mathbf{m}}
\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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\mathbf{m}}
- \tilde{\mathbf{m}}^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\mathbf{m}}
\Bigr\}
+ \mathrm{const.}
\end{aligned}
$$
適宜$\mathbf{x}_{*}$に影響しない項を$\mathrm{const.}$にまとめている。さらに、$\tilde{\mathbf{m}}$に式(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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} (
\boldsymbol{\Lambda}
\mathbf{x}_{*}
+ \boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\mathbf{m}
)
- \Bigl\{
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} (
\boldsymbol{\Lambda}
\mathbf{x}_{*}
+ \boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m}
)
\Bigr\}^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}
\mathbf{x}_{*}
- (\boldsymbol{\Lambda} \mathbf{x}_{*})^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\mathbf{m}
- (\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m})^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}
\mathbf{x}_{*}
- (\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m})^{\top}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}
\mathbf{x}_{*}
- \mathbf{x}_{*}^{\top}
\boldsymbol{\Lambda}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}_{\boldsymbol{\mu}} \mathbf{m}
- \mathbf{m}^{\top}
\boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}
\mathbf{x}_{*}
\Bigr\}
+ \mathrm{const.}
\\
&= - \frac{1}{2} \Bigl\{
\mathbf{x}_{*}^{\top} \Bigl(
\boldsymbol{\Lambda}
- \boldsymbol{\Lambda}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
\boldsymbol{\Lambda}
\Bigr)
\mathbf{x}_{*}
- 2 \mathbf{x}_{*}^{\top}
\boldsymbol{\Lambda}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
$$
【途中式の途中式】
- 波括弧を$\mathbf{a}$、丸括弧を$\mathbf{b}$として次のように考える。また逆行列の性質より、$\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} \tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} = \mathbf{I}_D$である。
$\mathbf{a}^{\top}$は$D$次元の横ベクトル、$\mathbf{b}$は$D$次元の縦ベクトルになるので、$\mathbf{a}^{\top} \mathbf{b}$はスカラとなる。よって、転置しても影響せず$\mathbf{a}^{\top} \mathbf{b} = (\mathbf{a}^{\top} \mathbf{b})^{\top}$である。さらに、転置の性質$(\mathbf{a}^{\top} \mathbf{b})^{\top} = \mathbf{b}^{\top} \mathbf{a}$より、後の因子は
$$
\begin{aligned}
\Bigl\{
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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\{
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-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}
\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} (
\boldsymbol{\Lambda}
\mathbf{x}_{*}
+ \boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\mathbf{m}
)
\end{aligned}
$$
となる。
あるいは、波括弧に関して${\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1} \mathbf{b}}^{\top} = \mathbf{b}^{\top} (\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})^{\top}$であり、また精度行列は対称行列なので$(\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1})^{\top} = \tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}$であることから、上の式となる。
- 後の因子の括弧を展開し、式を整理する。
- 同様に、$(\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}}$と変形できる。
- $\mathbf{x}_{*}^{\top} \mathbf{x}_{*}$と$\mathbf{x}_{*}$の項をそれぞれまとめて式を整理する。
- $\tilde{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}} = \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$を代入する。
となる。
事後分布のときと同様に、予測分布(3.107)を平均$\boldsymbol{\mu}_{*}$、精度$\boldsymbol{\Lambda}_{*}$の$D$次元ガウス分布
$$
p(\mathbf{x}_{*})
= \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}_{*}, \boldsymbol{\Lambda}_{*}^{-1})
\tag{3.108}
$$
とおき、対数をとって$\mathbf{x}_{*}$に関して整理すると
$$
\begin{aligned}
\ln p(\mathbf{x}_{*})
&= \ln \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}_{*}, \boldsymbol{\Lambda}_{*}^{-1})
\\
&= - \frac{1}{2} \Bigl\{
(\mathbf{x}_{*} - \boldsymbol{\mu}_{*})^{\top}
\boldsymbol{\Lambda}_{*}
(\mathbf{x}_{*} - \boldsymbol{\mu}_{*})
+ \ln |\boldsymbol{\Lambda}_{*}^{-1}|
+ D \ln 2 \pi
\Bigr\}
\\
&= - \frac{1}{2} \Biggl\{
\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\}
+ \mathrm{const.}
\\
&= - \frac{1}{2} \Bigl\{
\mathbf{x}_{*}^{\top}
\boldsymbol{\Lambda}_{*}
\mathbf{x}_{*}
- 2 \mathbf{x}_{*}^{\top}
\boldsymbol{\Lambda}_{*}
\boldsymbol{\mu}_{*}
\Bigr\}
+ \mathrm{const.}
\end{aligned}
$$
となる。
したがって、式(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}
\end{align}
$$
【途中式の途中式】
ウッドベリーの公式(A.7)の$\mathbf{U}, \mathbf{V}$を単位行列$\mathbf{I}$として考えると次のようになる。
$$
\begin{align}
(\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}
\tag{A.7}\\
\Rightarrow
(\mathbf{A} + \mathbf{B})^{-1}
&= \mathbf{A}^{-1}
- \mathbf{A}^{-1} (
\mathbf{B}^{-1}
+ \mathbf{A}^{-1}
)^{-1}
\mathbf{A}^{-1}
\end{align}
$$
$\boldsymbol{\Lambda}_{*}$の計算式が、この式の右辺の形をしているので、$\mathbf{A} = \boldsymbol{\Lambda}^{-1}, \mathbf{B} = \boldsymbol{\Lambda}_{\boldsymbol{\mu}}^{-1}$として左辺の形に変形する。
となる。
また、予測分布の平均パラメータは
$$
\begin{align}
\boldsymbol{\Lambda}_{*} \boldsymbol{\mu}_{*}
&= \boldsymbol{\Lambda} (
\boldsymbol{\Lambda}
+ \boldsymbol{\Lambda}_{\boldsymbol{\mu}}
)^{-1}
\boldsymbol{\Lambda}_{\boldsymbol{\mu}}
\mathbf{m}
\\
\Rightarrow
\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}
$$
【途中式の途中式】
「The Matrix Cookbook」の式(163)を用いる。
$$
(\mathbf{A}^{-1} + \mathbf{B}^{-1})^{-1}
= \mathbf{A} (\mathbf{A} + \mathbf{B})^{-1} \mathbf{B}
\tag{163}
$$
2行目の式の一部が、この式の右辺の形をしているので、$\mathbf{A} = \boldsymbol{\Lambda}$、$\mathbf{B} = \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$として左辺の形に変形する。また式(3.109)より、$\boldsymbol{\Lambda}_{*}$となる。
となる。
事後分布による予測分布
予測分布の計算に事前分布$p(\boldsymbol{\mu})$を用いることで、観測データによる学習を行っていない予測分布$p(\mathbf{x}_{*})$(のパラメータ$\boldsymbol{\mu}_{*}, \boldsymbol{\Lambda}_{*}$)を求めた。事後分布$p(\boldsymbol{\mu} | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(\mathbf{x}_{*} | \mathbf{X})$(のパラメータ$\hat{\boldsymbol{\mu}}_{*}, \hat{\boldsymbol{\Lambda}}_{*}$)を求められる。
よって、$\boldsymbol{\mu}_{*}, \boldsymbol{\lambda}_{*}$を構成する事前分布のパラメータ$\mathbf{m}, \boldsymbol{\Lambda}_{\boldsymbol{\mu}}$を事後分布のパラメータ(3.103)、(3.102)に置き換えると
$$
\begin{aligned}
\hat{\boldsymbol{\Lambda}}_{*}
&= (
\boldsymbol{\Lambda}^{-1}
+ \hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}^{-1}
)^{-1}
\\
&= \{
\boldsymbol{\Lambda}^{-1}
+ (
N \boldsymbol{\Lambda} + \boldsymbol{\Lambda}_{\boldsymbol{\mu}}
)^{-1}
\}^{-1}
\\
\hat{\boldsymbol{\mu}}_{*}
&= \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)
\end{aligned}
$$
が得られる。したがって、予測分布は平均$\hat{\boldsymbol{\mu}}_{*}$、精度$\hat{\boldsymbol{\Lambda}}_{*}$の$D$次元ガウス分布となる。
$$
p(\mathbf{x}_{*} | \mathbf{X})
= \mathcal{N}(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_{*}, \hat{\boldsymbol{\Lambda}}_{*}^{-1})
$$
また、上の式が予測分布のパラメータの計算式(更新式)である。
参考文献
最後の最後でちょろっとカンニングしちゃいました。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周年記念日！
おめでとうございます！！！！！！！！！！
- 2021/04/10：加筆修正しました。その際にRで実装編と記事を分割しました。
【次の節の内容】
www.anarchive-beta.com
