はじめに
『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。
この記事は、3.4.3項の内容です。「尤度関数を平均と精度が未知の多次元ガウス分布(多変量正規分布)」、「事前分布をガウス・ウィシャート分布」とした場合の「パラメータの事後分布」と「未観測値の予測分布」を導出します。
省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【前の節の内容】
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
3.4.3 平均・精度が未知の場合
多次元ガウス分布に従うと仮定する$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$を用いて、平均パラメータ$\boldsymbol{\mu}$と精度パラメータ$\boldsymbol{\Lambda}$の事後分布と未観測データ$\mathbf{x}_{*}$の予測分布を求めていく。
・観測モデルの設定
観測データ$\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}$である(固有値のことじゃないよ)。
尤度$p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\Lambda})$を$D$次元のガウス分布
$$
p(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
= \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
$$
とし、$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の事前分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$を、ガウス・ウィシャート分布
$$
\begin{align}
p(\boldsymbol{\mu}, \boldsymbol{\Lambda})
&= \mathrm{NW}(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{m}, \beta, \nu, \mathbf{W})
\\
&= \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
\tag{3.125}
\end{align}
$$
とする。ここで、$\mathbf{m},\ \beta$はガウス分布の平均パラメータと精度パラメータの係数であり、$\nu,\ $はウィシャート分布の自由度とスケール行列？である。自由度は$\nu > D - 1$の値をとり、$\mathbf{W}$は正定値行列である。
・事後分布の計算
観測データ$\mathbf{X}$が与えられた下での$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})$は、ベイズの定理を用いて
$$
\begin{align}
p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})
&= \frac{
p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
p(\boldsymbol{\boldsymbol{\mu}, \Lambda})
}{
p(\mathbf{X})
}
\tag{3.126}\\
&\propto
p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
p(\boldsymbol{\mu}, \boldsymbol{\Lambda})
\\
&= \left\{
\prod_{n=1}^N
p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda})
\right\}
p(\boldsymbol{\boldsymbol{\mu}, \Lambda})
\\
&= \left\{
\prod_{n=1}^N
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
\right\}
\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
\end{align}
$$
となる。分母の$p(\mathbf{X})$は$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$に影響しないため省略して、比例関係のみに注目する。省略した部分については、最後に正規化することで対応できる。
またパラメータの依存関係から、左辺の事後分布は
$$
p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})
= p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X})
p(\boldsymbol{\Lambda} | \mathbf{X})
\tag{3.127}
$$
と分解できる。
・平均の事後分布の計算
まずは、$\boldsymbol{\mu}$の事後分布$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X})$を導出する。
$\mathbf{X}$が与えられた下での$\boldsymbol{\mu}$の事後分布は、式(3.126)、(3.127)より
$$
\begin{aligned}
p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X})
&= \frac{
\left\{
\prod_{n=1}^N
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
\right\}
\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
}{
p(\boldsymbol{\Lambda} | \mathbf{X})
p(\mathbf{X})
}
\\
&\propto
\left\{
\prod_{n=1}^N
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
\right\}
\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\end{aligned}
$$
となる。$\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})$と$p(\boldsymbol{\Lambda} | \mathbf{X})$も$\boldsymbol{\mu}$に影響しないため省く。
$p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X})$と3.4.1項「平均が未知の場合」の式(3.98)を比べると、事前分布の精度パラメータが$\boldsymbol{\Lambda}_{\boldsymbol{\mu}}$から$\beta \boldsymbol{\Lambda}$に置き換わっただけで同じ形状である。よって同様の手順で求められるので、$\boldsymbol{\Lambda}_{\boldsymbol{\mu}}$と$\beta \boldsymbol{\Lambda}$を置き換えればよい。
平均が未知の場合の$\boldsymbol{\mu}$の事後分布の精度パラメータ
$$
\hat{\boldsymbol{\Lambda}}_{\boldsymbol{\mu}}
= N \boldsymbol{\Lambda}
+ \beta \boldsymbol{\Lambda}
\tag{3.102}
$$
について$\boldsymbol{\Lambda}_{\boldsymbol{\mu}}$を$\beta \boldsymbol{\Lambda}$に置き換えると
$$
\begin{align}
\hat{\beta} \boldsymbol{\Lambda}
&= N \boldsymbol{\Lambda}
+ \beta \boldsymbol{\Lambda}
\\
\hat{\beta}
&= N + \beta
\tag{3.129.a}
\end{align}
$$
精度パラメータの係数$\hat{\beta}$の計算式(更新式)が得られる。
同様に、$\boldsymbol{\mu}$の事後分布の平均パラメータ
$$
\hat{\mathbf{m}}
= (\hat{\beta} \boldsymbol{\Lambda})^{-1} \left(
\boldsymbol{\Lambda}
\sum_{n=1}^N \mathbf{x}_n
+ \beta \boldsymbol{\Lambda} \mathbf{m}
\right)
\tag{3.103}
$$
についても置き換えると
$$
\begin{align}
\hat{\mathbf{m}}
&= (\hat{\beta} \boldsymbol{\Lambda})^{-1} \left(
\boldsymbol{\Lambda}
\sum_{n=1}^N \mathbf{x}_n
+ \beta \boldsymbol{\Lambda} \mathbf{m}
\right)
\\
&= \frac{1}{\hat{\beta}}
\boldsymbol{\Lambda}^{-1} \left(
\boldsymbol{\Lambda}
\sum_{n=1}^N \mathbf{x}_n
+ \beta \boldsymbol{\Lambda} \mathbf{m}
\right)
\\
&= \frac{1}{\hat{\beta}} \left(
\sum_{n=1}^N \mathbf{x}_n
+ \beta \mathbf{m}
\right)
\tag{3.129.b}
\end{align}
$$
平均パラメータ$\hat{\mathbf{m}}$の計算式(更新式)が得られる。
したがって、$\boldsymbol{\mu}$の事後分布は
$$
p(\boldsymbol{\mu} | \mathbf{X})
= \mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, (\hat{\beta} \boldsymbol{\Lambda})^{-1})
$$
平均$\hat{\mathbf{m}}$、精度$\hat{\beta} \boldsymbol{\Lambda}$の$D$次元ガウス分布となることが分かる。
・精度の事後分布の計算
続いて、$\boldsymbol{\Lambda}$の事後分布$p(\boldsymbol{\Lambda} | \mathbf{X})$を導出する。
$\mathbf{X}$が与えられた下での$\boldsymbol{\Lambda}$の事後分布も、式(3.126)、(3.127)より
$$
\begin{aligned}
p(\boldsymbol{\Lambda} | \mathbf{X})
&= \frac{
\left\{
\prod_{n=1}^N
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
\right\}
\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
}{
p(\boldsymbol{\mu} | \boldsymbol{\Lambda}, \mathbf{X})
p(\mathbf{X})
}
\\
&\propto
\frac{
\left\{
\prod_{n=1}^N
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
\right\}
\mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
\mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
}{
\mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, (\hat{\beta} \boldsymbol{\Lambda})^{-1})
}
\end{aligned}
$$
となる。
この分布の具体的な形状を明らかにしていく。対数をとって指数部分の計算を分かりやすくして、$\boldsymbol{\Lambda}$に関して整理する。
$$
\begin{aligned}
\ln p(\boldsymbol{\Lambda} | \mathbf{X})
&= \sum_{n=1}^N
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})
+ \ln \mathcal{N}(\boldsymbol{\mu} | \mathbf{m}, (\beta \boldsymbol{\Lambda})^{-1})
- \ln \mathcal{N}(\boldsymbol{\mu} | \hat{\mathbf{m}}, (\hat{\beta} \boldsymbol{\Lambda})^{-1})
+ \ln \mathcal{W}(\boldsymbol{\Lambda} | \nu, \mathbf{W})
+ \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}
\beta \boldsymbol{\Lambda}
(\boldsymbol{\mu} - \mathbf{m})
+ \ln |(\beta \boldsymbol{\Lambda})^{-1}|
+ D \ln 2 \pi
\Bigr\} \\
&\qquad
+ \frac{1}{2} \Bigl\{
(\boldsymbol{\mu} - \hat{\mathbf{m}})^{\top}
\hat{\beta} \boldsymbol{\Lambda}
(\boldsymbol{\mu} - \hat{\mathbf{m}})
+ \ln |(\hat{\beta} \boldsymbol{\Lambda})^{-1}|
+ D \ln 2 \pi
\Bigr\} \\
&\qquad
+ \frac{\nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
+ \ln C_{\mathcal{W}}(\nu, \mathbf{W})
+ \mathrm{const.}
\\
&= - \frac{1}{2} \left\{
\sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n
- 2 \sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ N \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- N \ln |\boldsymbol{\Lambda}|
\right\} \\
&\qquad
- \frac{1}{2} \Bigl\{
\beta
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- 2 \beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ \beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
- \ln |\beta \boldsymbol{\Lambda}|
\Bigr\} \\
&\qquad
+ \frac{1}{2} \Bigl\{
\hat{\beta}
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- 2 \hat{\beta}
\hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ \hat{\beta}
\hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \hat{\mathbf{m}}
- \ln |\hat{\beta} \boldsymbol{\Lambda}|
\Bigr\} \\
&\qquad
+ \frac{\nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
+ \mathrm{const.}
\end{aligned}
$$
適宜$\boldsymbol{\Lambda}$に影響しない項を$\mathrm{const.}$にまとめる。また、3つの精度行列の行列式の対数について、行列式の性質$|\mathbf{A}^{-1}| = |\mathbf{A}|^{-1}$と自然対数の性質$\ln x^{-1} = - \ln x$より、$\ln |\boldsymbol{\Lambda}^{-1}| = - \ln |\boldsymbol{\Lambda}|$となる。さらに、3行目の1つ目と2つ目の項の$\hat{\beta},\ \hat{\mathbf{m}}$に式(3.129)を代入すると
$$
\begin{align}
\ln p(\boldsymbol{\Lambda} | \mathbf{X})
&= - \frac{1}{2} \left\{
\sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n
- 2 \sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ N \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- N \ln |\boldsymbol{\Lambda}|
\right\} \\
&\qquad
- \frac{1}{2} \Bigl\{
\beta
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- 2 \beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ \beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
- \ln \beta^D
- \ln |\boldsymbol{\Lambda}|
\Bigr\} \\
&\qquad
+ \frac{1}{2} \left\{
(N + \beta)
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- 2 \Bigl(
\sum_{n=1}^N \mathbf{x}_n + \beta \mathbf{m}
\Bigr)^{\top}
\boldsymbol{\Lambda} \boldsymbol{\mu}
+ \hat{\beta}
\hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \hat{\mathbf{m}}
- \ln \hat{\beta}^D
- \ln |\boldsymbol{\Lambda}|
\right\} \\
&\qquad
+ \frac{\nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
+ \mathrm{const.}
\\
&= \frac{N}{2} \ln |\boldsymbol{\Lambda}|
- \frac{1}{2} \left\{
\sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n
+ \beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
- \hat{\beta}
\hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \hat{\mathbf{m}}
\right\} \\
&\qquad
+ \frac{\nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
+ \mathrm{const.}
\\
&= \frac{N + \nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr} \left[
\left\{
\sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top}
+ \beta \mathbf{m} \mathbf{m}^{\top}
+ \hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top}
+ \mathbf{W}^{^-1}
\right\}
\boldsymbol{\Lambda}
\right]
+ \mathrm{const.}
\tag{3.131}
\end{align}
$$
【途中式の途中式】(クリックで展開)
- 行列式の性質$|c \mathbf{A}| = c^D |\mathbf{A}|$と自然対数の性質$\ln x y = \ln x + \ln y$より、$\ln |\beta \boldsymbol{\Lambda}| = \ln \beta^D + \ln |\boldsymbol{\Lambda}|$となる。
- 式を整理する。
- 波括弧内の項を$\mathrm{Tr}(\cdot)$に置き換え、$\ln |\boldsymbol{\Lambda}|$と$\boldsymbol{\Lambda}$の項をそれぞれまとめて式を整理する。
3.4.2項の事後分布の導出時と同様に、波括弧内の1つ目の項は
$$
\begin{aligned}
\sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n
&= \sum_{n=1}^N
\sum_{d=1}^D \sum_{d'=1}^D
x_{n,d'} \lambda_{d,d'} x_{n,d}
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\begin{bmatrix}
\sum_{d=1}^D x_{n,1} x_{n,d} \lambda_{d,1} & \cdots & \sum_{d=1}^D x_{n,1} x_{n,d} \lambda_{d,D} \\
\vdots & \ddots & \vdots \\
\sum_{d=1}^D x_{n,D} x_{n,d} \lambda_{d,1} & \cdots & \sum_{d=1}^D x_{n,D} x_{n,d} \lambda_{d,D}
\end{bmatrix}
\right)
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\begin{bmatrix}
x_{n,1} x_{n,1} & \cdots & x_{n,1} x_{n,D} \\
\vdots & \ddots & \vdots \\
x_{n,D} x_{n,1} & \cdots & x_{n,D} 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}
\right)
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\begin{bmatrix}
x_{n,1} \\
\vdots \\
x_{n,D}
\end{bmatrix}
\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}
\right)
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\mathbf{x}_n \mathbf{x}_n^{\top} \boldsymbol{\Lambda}
\right)
\end{aligned}
$$
と置き換えられる。同様に、2つ目の項は
$$
\begin{aligned}
\beta
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
&= \beta
\sum_{d=1}^D \sum_{d'=1}^D
m_{d'} \lambda_{dd'} m_d
\\
&= \mathrm{Tr} \left(
\beta
\begin{bmatrix}
\sum_{d=1}^D m_1 m_d \lambda_{d,1} & \cdots & \sum_{d=1}^D m_1 m_d \lambda_{d,D} \\
\vdots & \ddots & \vdots \\
\sum_{d=1}^D m_D m_d \lambda_{d,1} & \cdots & \sum_{d=1}^D m_D m_d \lambda_{d,D}
\end{bmatrix}
\right)
\\
&= \mathrm{Tr}(
\beta
\mathbf{m} \mathbf{m}^{\top} \boldsymbol{\Lambda}
)
\end{aligned}
$$
と変形でき、3つ目の項も$\hat{\beta} \hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \hat{\mathbf{m}} = \mathrm{Tr}(\hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda})$と変形できる。よってトレースの性質$\mathrm{Tr}(\mathbf{A}) + \mathrm{Tr}(\mathbf{B}) = \mathrm{Tr}(\mathbf{A} + \mathbf{B})$より、波括弧内の項は
$$
\begin{aligned}
&\sum_{n=1}^N
\mathbf{x}_n^{\top} \boldsymbol{\Lambda} \mathbf{x}_n
+ \beta \mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
- \hat{\beta} \hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \hat{\mathbf{m}}
+ \mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\mathbf{x}_n \mathbf{x}_n^{\top} \boldsymbol{\Lambda}
\right)
+ \mathrm{Tr}(\beta \mathbf{m} \mathbf{m}^{\top} \boldsymbol{\Lambda})
+ \mathrm{Tr}(\hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda})
+ \mathrm{Tr}(\mathbf{W}^{-1} \Lambda)
\\
&= \mathrm{Tr} \left(
\sum_{n=1}^N
\mathbf{x}_n \mathbf{x}_n^{\top} \boldsymbol{\Lambda}
+ \beta \mathbf{m} \mathbf{m}^{\top} \boldsymbol{\Lambda}
+ \hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top} \boldsymbol{\Lambda}
+ \mathbf{W}^{^-1} \boldsymbol{\Lambda}
\right)
\end{aligned}
$$
で置き換えられる。
となる。
式(3.131)について
$$
\begin{aligned}
\hat{\mathbf{W}}^{-1}
&= \sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top}
+ \beta \mathbf{m} \mathbf{m}^{\top}
- \hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top}
+ \mathbf{W}^{-1}
\\
\hat{\nu}
&= N + \nu
\end{aligned}
\tag{3.133}
$$
とおき
$$
\ln p(\boldsymbol{\Lambda} | \mathbf{X})
= \frac{\hat{\nu} - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
- \frac{1}{2}
\mathrm{Tr}(\hat{\mathbf{W}}^{-1} \boldsymbol{\Lambda})
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化する)と
$$
p(\boldsymbol{\Lambda} | \mathbf{X})
= \mathcal{W}(\boldsymbol{\Lambda} | \hat{\nu}, \hat{\mathbf{W}})
\tag{3.132}
$$
$\boldsymbol{\Lambda}$の事後分布は式の形から、パラメータ$\hat{\nu},\ \hat{\mathbf{W}}$を持つはウィシャート分布となることが分かる。
また、式(1.133)が超パラメータ$\hat{\nu},\ \hat{\mathbf{W}}$の計算式(更新式)である。
・予測分布の計算
最後に、多次元ガウス分布に従う未観測データ$\mathbf{x}_{*} = (x_{*,1}, x_{*,2}, \cdots, x_{*,D})^{\top}$に対する予測分布を導出する。
先に、事前分布(観測データによる学習を行っていない分布)$p(\mu)$を用いて、未学習の予測分布$p(\mathbf{x}_{*})$を求める。その結果を用いて、学習後の予測分布$p(\mathbf{x}_{*} | \mathbf{X})$を求める。
1次元ガウス分布のときと同様に、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Lambda}$を周辺化(積分計算)して
$$
p(\mathbf{x}_{*})
= \iint
p(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
p(\boldsymbol{\mu}, \boldsymbol{\Lambda})
d\boldsymbol{\mu} d\boldsymbol{\Lambda}
\tag{3.134}
$$
予測分布を求めることは避け、ベイズの定理を用いて
$$
p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})
= \frac{
p(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
p(\boldsymbol{\mu}, \boldsymbol{\Lambda})
}{
p(\mathbf{x}_{*})
}
$$
予測分布を求める。
この式の両辺の対数をとり
$$
\ln p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})
= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
+ \ln p(\boldsymbol{\mu}, \boldsymbol{\Lambda})
- \ln p(\mathbf{x}_{*})
$$
予測分布に関して式を整理すると
$$
\ln p(\mathbf{x}_{*})
= \ln p(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
- \ln p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})
+ \mathrm{const.}
\tag{3.135}
$$
となる。ただし、$\mathbf{x}_{*}$に影響しない$\ln p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$を$\mathrm{const.}$とおいた。
この式から予測分布の具体的な式を計算する。
$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})$は、1つのデータ$\mathbf{x}_{*}$が与えられた下での$\boldsymbol{\mu}, \boldsymbol{\Lambda}$の条件付き分布である。つまり$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})$は、$N$個の観測データが与えられた下での事後分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})$と同様の手順で求められる(同様のパラメータになる)。
したがって、事後分布のパラメータ(3.129)、(3.133)を用いると、$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | x_{*})$は$N = 1$より
$$
\begin{align}
p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{x}_{*})
&= \mathrm{NW}(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \tilde{\mathbf{m}}, \{(1 + \beta), 1 + \nu, \tilde{\mathbf{W}})
\\
&= \mathcal{N}(\boldsymbol{\mu} | \tilde{\mathbf{m}}, \{(1 + \beta) \boldsymbol{\Lambda}\}^{-1})
\mathcal{W}(\boldsymbol{\lambda} | 1 + \nu, \tilde{\mathbf{W}})
\tag{1.136}
\end{align}
$$
となる。ただし
$$
\begin{align}
\tilde{\mathbf{m}}
&= \frac{1}{1 + \beta}
(\mathbf{x}_{*} + \beta \mathbf{m})
\tag{3.137.a}\\
\tilde{\mathbf{W}}^{-1}
&= \mathbf{x}_{*} \mathbf{x}_{*}^{\top}
+ \beta \mathbf{m} \mathbf{m}^{\top}
- \hat{\beta} \hat{\mathbf{m}} \hat{\mathbf{m}}^{\top}
+ \mathbf{W}^{-1}
\\
&= \mathbf{x}_{*} \mathbf{x}_{*}^{\top}
+ \beta \mathbf{m} \mathbf{m}^{\top}
- \frac{1}{1 + \beta}
(\mathbf{x}_{*} + \beta \mathbf{m})
(\mathbf{x}_{*} + \beta \mathbf{m})^{\top}
+ \mathbf{W}^{-1}
\\
&= \mathbf{x}_{*} \mathbf{x}_{*}^{\top}
+ \beta \mathbf{m} \mathbf{m}^{\top}
- \frac{1}{1 + \beta}
\mathbf{x}_{*} \mathbf{x}_{*}^{\top}
- \frac{2 \beta}{1 + \beta}
\mathbf{x}_{*} \mathbf{m}^{\top}
- \frac{\beta^2}{1 + \beta}
\mathbf{m} \mathbf{m}^{\top}
+ \mathbf{W}^{-1}
\\
&= \frac{\beta}{1 + \beta}
\mathbf{x}_{*} \mathbf{x}_{*}^{\top}
- \frac{2 \beta}{1 + \beta}
\mathbf{x}_{*} \mathbf{m}^{\top}
+ \frac{\beta}{1 + \beta}
\mathbf{m} \mathbf{m}^{\top}
+ \mathbf{W}^{-1}
\\
&= \frac{\beta}{1 + \beta}
(\mathbf{x}_{*} - \mathbf{m})
(\mathbf{x}_{*} - \mathbf{m})^{\top}
+ \mathbf{W}^{-1}
\tag{3.137.b}
\end{align}
$$
とおく。分かりにくければ、さらに$\tilde{\beta} = 1 + \beta$、$\tilde{\nu} = 1 + \nu$とおいてもよい。
尤度と式(3.136)を式(3.135)に代入して、$\mathbf{x}_{*}$に関して式を整理する。
$$
\begin{aligned}
\ln p(\mathbf{x}_{*})
&= \ln \mathcal{N}(\mathbf{x}_{*} | \boldsymbol{\mu}, \boldsymbol{\Lambda})
- \ln \mathcal{N}(\boldsymbol{\mu} | \tilde{\mathbf{m}}, \{(1 + \beta) \boldsymbol{\Lambda}\}^{-1})
- \ln \mathcal{W}(\boldsymbol{\Lambda} | 1 + \nu, \tilde{\mathbf{W}})
+ \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}
(1 + \beta) \boldsymbol{\Lambda}
(\boldsymbol{\mu} - \tilde{\mathbf{m}})
+ \ln |\{(1 + \beta) \boldsymbol{\Lambda}\}^{-1}|
+ D \ln 2 \pi
\Bigr\} \\
&\qquad
- \frac{1 + \nu - D - 1}{2}
\ln |\boldsymbol{\Lambda}|
+ \frac{1}{2}
\mathrm{Tr}(\tilde{\mathbf{W}}^{-1} \boldsymbol{\Lambda}) \\
&\qquad
+ \frac{1 + \nu}{2}
\ln |\tilde{\mathbf{W}}|
+ \frac{(1 + \nu) D}{2} \ln 2
+ \frac{D (D - 1)}{4} \ln \pi
+ \sum_{d=1}^D
\ln \Gamma \left(\frac{1 + \nu + 1 - d}{2} \right)
+ \mathrm{const.}
\\
&= - \frac{1}{2} \Bigl\{
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*}
- 2 \mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
+ \boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
\Bigr\} \\
&\qquad
+ \frac{1}{2} \Bigl\{
(1 + \beta)
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \boldsymbol{\mu}
- 2 (1 + \beta)
\boldsymbol{\mu}^{\top} \boldsymbol{\Lambda} \tilde{\mathbf{m}}
+ (1 + \beta)
\tilde{\mathbf{m}}^{\top} \boldsymbol{\Lambda} \tilde{\mathbf{m}}
\Bigr\} \\
&\qquad
+ \frac{1}{2}
\mathrm{Tr} \left(
\frac{\beta}{1 + \beta}
(\mathbf{x}_{*} - \mathbf{m})
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\boldsymbol{\Lambda}
+ \mathbf{W}^{-1}
\boldsymbol{\Lambda}
\right) \\
&\qquad
+ \frac{1 + \nu}{2}
\ln \Bigl|
\mathbf{W}^{-1}
+ \frac{\beta}{1 + \beta}
(\mathbf{x}_{*} - \mathbf{m})
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\Bigr|^{-1}
+ \mathrm{const.}
\end{aligned}
$$
適宜$\mathbf{x}_{*}$に影響しない項を$\mathrm{const.}$にまとめる。また2つ目から3つ目の式変形では、$|\tilde{\mathbf{W}}| = |\tilde{\mathbf{W}}^{-1}|^{-1}$として、$\tilde{\mathbf{W}}^{-1}$に式(3.137.b)を代入した。この式をさらに整理すると
$$
\begin{align}
\ln p(\mathbf{x}_{*})
&= - \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}_{*} + \beta \mathbf{m})
- \frac{1}{1 + \beta}
(\mathbf{x}_{*} + \beta \mathbf{m})^{\top}
\boldsymbol{\Lambda}
(\mathbf{x}_{*} + \beta \mathbf{m})
\Bigr\} \\
&\qquad
+ \frac{1}{2}
\frac{\beta}{1 + \beta}
\mathrm{Tr} \Bigl(
(\mathbf{x}_{*} - \mathbf{m})
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\boldsymbol{\Lambda}
\Bigr)
+ \frac{1}{2}
\mathrm{Tr}(\mathbf{W}^{-1} \boldsymbol{\Lambda}) \\
&\qquad
- \frac{1 + \nu}{2}
\ln |\mathbf{W}^{-1}| \left|
\mathbf{I}_D
+ \frac{\beta}{1 + \beta}
\mathbf{W}
(\mathbf{x}_{*} - \mathbf{m})
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\right|
+ \mathrm{const.}
\\
&= - \frac{1}{2} \Bigl\{
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*}
- \frac{1}{1 + \beta}
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*}
- \frac{2 \beta}{1 + \beta}
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{m}
- \frac{\beta^2}{1 + \beta}
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
\Bigr\} \\
&\qquad
+ \frac{1}{2} \left\{
\frac{\beta}{1 + \beta}
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{x}_{*}
- \frac{2 \beta}{1 + \beta}
\mathbf{x}_{*}^{\top} \boldsymbol{\Lambda} \mathbf{m}
+ \frac{\beta}{1 + \beta}
\mathbf{m}^{\top} \boldsymbol{\Lambda} \mathbf{m}
\right\} \\
&\qquad
- \frac{1 + \nu}{2}
\ln \left|
\mathbf{I}_1
+ \frac{\beta}{1 + \beta}
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\mathbf{W}^{\top}
(\mathbf{x}_{*} - \mathbf{m})
\right|
- \frac{1 + \nu}{2}
\ln |\mathbf{W}^{-1}|
+ \mathrm{const.}
\\
&= - \frac{1 + \nu}{2}
\ln \left\{
1 + \frac{\beta}{1 + \beta}
(\mathbf{x}_{*} - \mathbf{m})^{\top}
\mathbf{W}
(\mathbf{x}_{*} - \mathbf{m})
\right\}
+ \mathrm{const.}
\tag{3.138}
\end{align}
$$
【途中式の途中式】
- 項を展開または分割する。
- $\tilde{\mathbf{m}}$に式(3.137.b)を代入する。
- トレースの項は、トレースの性質$\mathrm{Tr}(\mathbf{A}) + \mathrm{Tr}(\mathbf{B}) = \mathrm{Tr}(\mathbf{A} + \mathbf{B})$を用いる。
- 行列式の項は、行列式の性質$|\mathbf{A} \mathbf{B}| = |\mathbf{A}| |\mathbf{B}|$と逆行列の性質$\mathbf{A} \mathbf{A}^{-1} = \mathbf{I}$を用いる。
- 項を展開または変形する。
- 事後分布の導出時に確認した関係より、$\mathrm{Tr}( (\mathbf{x}_{*} - \mathbf{m}) (\mathbf{x}_{*} - \mathbf{m})^{\top} \boldsymbol{\Lambda} ) = (\mathbf{x}_{*} - \mathbf{m})^{\top} \boldsymbol{\Lambda} (\mathbf{x}_{*} - \mathbf{m})$である。
- $\mathbf{W} (\mathbf{x}_{*} - \boldsymbol{\mu})$を1つの行列とみると、$\mathbf{W} (\mathbf{x}_{*} - \mu) (\mathbf{x}_{*} - \mu)^{\top}$は$D \times 1$、$1 \times D$の行列の積である。そこで、行列式の性質$|\mathbf{I}_N + \mathbf{C}^{\top} \mathbf{D}| = |\mathbf{I}_M + \mathbf{C} \mathbf{D}^{\top}|$の変形を行う。
- $a \ln x y = a (\ln x + \ln y) = a \ln x + a \ln y$の変形を行う。
- 式を整理する。
- 行列式の中がスカラになるので、行列式の定義より$|a| = a$である。
- $(\mathbf{x}_{*} - \boldsymbol{\mu})^{\top} \mathbf{W}^{\top} (\mathbf{x}_{*} - \boldsymbol{\mu})$はスカラになるので転置できるため、$(\mathbf{A} \mathbf{B} \mathbf{C})^{\top} = \mathbf{C}^{\top} \mathbf{B}^{\top} \mathbf{A}^{\top}$の変形を行う。
となる。
式(3.138)について
$$
\begin{aligned}
\boldsymbol{\mu}_s
&= \mathbf{m}
\\
\boldsymbol{\Lambda}_s
&= \frac{
\nu_s \beta
}{
1 + \beta
}
\mathbf{W}
\\
&= \frac{
(1 - D + \nu) \beta
}{
1 + \beta
}
\mathbf{W}
\\
\nu_s
&= 1 - D + \nu
\end{aligned}
\tag{3.140}
$$
とおき
$$
\ln \mathrm{St}(\mathbf{x} | \boldsymbol{\mu}_s, \boldsymbol{\Lambda}_s, \nu_s)
= - \frac{\nu_s + D}{2}
\ln \left\{
1
+ \frac{1}{\nu_s}
(\mathbf{x}_{*} - \boldsymbol{\mu}_s)^{\top}
\boldsymbol{\Lambda}_s
(\mathbf{x}_{*} - \boldsymbol{\mu}_s)
\right\}
+ \mathrm{const.}
$$
さらに$\ln$を外し、$\mathrm{const.}$を正規化項に置き換える(正規化項する)と
$$
p(\mathbf{x}_{*})
= \mathrm{St}(\mathbf{x}_{*} | \boldsymbol{\mu}_s, \boldsymbol{\Lambda}_s, \nu_s)
\tag{3.139}
$$
予測分布は式の形状から、パラメータ$\boldsymbol{\mu}_s,\ \boldsymbol{\Lambda}_s,\ \nu_s$を持つ多次元のスチューデントのt分布となることが分かる。
予測分布の計算に事前分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda})$を用いることで、観測データによる学習を行っていない予測分布$p(\mathbf{x}_{*})$(のパラメータ$\boldsymbol{\mu}_s,\ \boldsymbol{\Lambda}_s,\ \nu_s$)を求めた。事後分布$p(\boldsymbol{\mu}, \boldsymbol{\Lambda} | \mathbf{X})$を用いると、同様の手順で観測データ$\mathbf{X}$によって学習した予測分布$p(\mathbf{x}_{*} | \mathbf{X})$(のパラメータ$\hat{\boldsymbol{\mu}}_s,\ \hat{\boldsymbol{\Lambda}}_s,\ \hat{\nu}_s$)を求められる。
よって、$\boldsymbol{\mu}_s,\ \boldsymbol{\Lambda}_s,\ \nu_s$を構成する事前分布のパラメータ$\beta,\ \mathbf{m},\ \nu,\ \mathbf{W}$を事後分布のパラメータ(3.129)、(1.33)に置き換えると
$$
\begin{aligned}
\hat{\boldsymbol{\mu}}_s
&= \hat{\mathbf{m}}
\\
&= \frac{1}{N + \beta} \left(
\sum_{n=1}^N \mathbf{x}_n
+ \beta \mathbf{m}
\right)
\\
\hat{\boldsymbol{\Lambda}}_s
&= \frac{
(1 - D + \hat{\nu}) \hat{\beta}
}{
1 + \hat{\beta}
}
\hat{\mathbf{W}}
\\
&= \frac{
(1 - D + \hat{\nu}) (N + \beta)
}{
N + 1 + \beta
}
\hat{\mathbf{W}}
\\
\hat{\nu}_s
&= 1 - D + \hat{\nu}
\\
&= N + 1 - D + \nu
\end{aligned}
$$
が得られる。したがって、予測分布はパラメータ$\hat{\boldsymbol{\mu}}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s$を持つ多次元のスチューデントのt分布となる。
$$
p(\mathbf{x}_{*})
= \mathrm{St}(\mathbf{x}_{*} | \hat{\boldsymbol{\mu}}_s, \hat{\boldsymbol{\Lambda}}_s, \hat{\nu}_s)
$$
また、上の式が予測分布のパラメータの計算式 (更新式) である。
参考文献
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
多次元ガウス分布のベイズ推論完了！ただただ愚直に解きました。もっとスマートに導出できたりするのでしょうか？
そろそろPRMLに移ってもいいもんでしょうかね。次の線形回帰ができたら試しに読んでみる(2度目の正直)。
- 2021/04/12：加筆修正しました。その際にRで実装編と記事を分割しました。
現在PRML三度目の挑戦に詰まって放置中、、、
【次の節の内容】
www.anarchive-beta.com
