PRML 読書会 #13 「10.2 変分混合ガウス分布」資料(1)

「パターン認識と機械学習」(PRML)読書会 #13 で担当する 10.2「変分混合ガウス分布」と 10.2.1「変分事後分布」の資料です。その２はこちら。

10.1 で述べた変分ベイズ近似(Variational Bayesian)を混合ガウス分布の推論に適応する例を見る。

X = {x_n} : 観測変数
Z = {z_n} : 対応する潜在変数、 1-of-K
N: データ件数
D: x の次元

★ PRML 上巻の付録Ｂは「公式」として使います

変分ベイズの手順(PRML 10.1 の復習)

(1) パラメータに事前分布を入れる(ベイズ！)
(2) 事後分布を「変分近似」する
- 事後分布を q(Z) とおく
- 適当な分解を仮定(真には分解できない=近似)
(3) 仮定した分解に基づき、変分法で個々の因子を導出する
- パラメータの更新式を求め、繰り返しによって解く

VB の計算のキモ(PRML 10.1 復習)

「変分近似」 q(Z) = Πq_i(Z_i) に対して、q_j の最適解(＝真の事後分布とのKLを最小にするもの)である q_j^*(Z_j) は (10.9) 式で求められる

$\ln q_j^*(Z_j) = \mathbb{E}_{i \ne j}[\ln p(X, Z)] + \rm{const}$

パラメータ Z_i について最適化したい！
同時分布の対数を「 Z_i 以外の全てのパラメータ」(観測変数除く)について期待値を取る
同時分布の対数をバラして、Z_i に依存しない項は全部 const. に押し込む
残った項をがんばって計算

このパターンをパラメータ(群)の数だけ繰り返し。

モデル(条件付き分布)

潜在変数 Z の条件付き分布

$p(\boldsymbol{Z}|\boldsymbol{\pi})=\prod_{n=1}^N \prod_{k=1}^K \pi_k^{z_{nk}}$

多項分布(1-of-K)
π:混合比

観測変数 X の条件付き分布

$p(\boldsymbol{X}|\boldsymbol{Z},\boldsymbol{\mu},\boldsymbol{\Lambda})=\prod_{n=1}^N \prod_{k=1}^K \mathcal{N}(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k^{-1})^{z_{nk}}$

混合ガウス分布
z_n により混合要素を指定

事前分布

πは多項分布のハイパーパラメータなので、その共役事前分布であるディリクレ分布を入れる

$p(\boldsymbol{\pi})=\rm{Dir}(\boldsymbol{\pi}|\boldsymbol{\alpha}_0) = C(\boldsymbol{\alpha}_0) \prod_{k=1}^K \pi_k^{\alpha_0-1}$

where $C(\boldsymbol{\alpha}) = \frac{\Gamma(\alpha_1 +...+ \alpha_K )}{\Gamma(\alpha_1)...\Gamma(\alpha_K)}$

ディリクレ分布のパラメータ α_0 は、各混合要素に対する事前の観測回数に相当

★PRML では対称にするように書かれているけど、傾斜させれば k 番目を最大クラスタになりやすくしたい、とかできるんだろうか？

μ、Λはガウス分布の平均と精度なので、その共役事前分布であるガウス−ウィシャート分布を入れる

$\begin{eqnarray} &&p(\boldsymbol{\mu},\boldsymbol{\Lambda})\\ &=& p(\boldsymbol{\mu}|\boldsymbol{\Lambda})p(\boldsymbol{\Lambda}) \\ &=& \prod_{k=1}^K \mathcal{N}(\boldsymbol{\mu}_k|\boldsymbol{m}_0,(\beta\boldsymbol{\Lambda}_k)^{-1})\mathcal{W}(\boldsymbol{\Lambda}_k|\boldsymbol{W}_0,\nu_0)\end{eqnarray}$

where $\mathcal{W}(\boldsymbol{\Lambda}|\boldsymbol{W},\nu) = B(\boldsymbol{W},\nu)|\boldsymbol{\Lambda}|^{(\nu-D-1)/2}\exp\left(-\frac{1}{2}\rm{Tr}(\boldsymbol{W}^{-1}\boldsymbol{\Lambda})\right)$ ,
$B(\boldsymbol{W},\nu) = |\boldsymbol{W}|^{-\nu/2}\left\{ 2^{\nu D/2} \pi^{D(D-1)/4} \prod_{i=1}^D \Gamma\left( \frac{\nu+1-i}{2} \right) \right\}^{-1}$

同時分布

$p(\boldsymbol{X},\boldsymbol{Z},\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})=p(\boldsymbol{X}|\boldsymbol{Z},\boldsymbol{\mu},\boldsymbol{\Lambda})p(\boldsymbol{Z}|\boldsymbol{\pi})p(\boldsymbol{\pi})p(\boldsymbol{\mu}|\boldsymbol{\Lambda})p(\boldsymbol{\Lambda})$

事後分布

変分近似: (真の)事後分布を潜在変数 Z とパラメータ π, μ, Λ に分解する。

$q(\boldsymbol{Z},\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})=q(\boldsymbol{Z})q(\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})$

「ベイズ混合モデルについて実際に計算可能な解を得る上で、必要な唯一の仮定」
なにげにこの２つの式が一番重要。

q(Z) を推定する

$\begin{eqnarray}\ln q^*(\boldsymbol{Z}) &=& \mathbb{E}_{\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda}}[\ln p(\boldsymbol{X},\boldsymbol{Z},\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})] + \rm{const} \\&=& \mathbb{E}_{\boldsymbol{\pi}}[\ln p(\boldsymbol{Z}|\boldsymbol{\pi})] + \mathbb{E}_{\boldsymbol{\mu},\boldsymbol{\Lambda}}[\ln p(\boldsymbol{X}|\boldsymbol{Z},\boldsymbol{\mu},\boldsymbol{\Lambda})] + \rm{const} \\&=& \sum_{n=1}^N \sum_{k=1}^K z_{nk}\ln \rho_{nk} + \rm{const}\end{eqnarray}$

where $\begin{eqnarray} \ln \rho_{nk} &=& \mathbb{E}[\ln\pi_k]+\frac{1}{2}\mathbb{E}[\ln|\boldsymbol{\Lambda}_k|]-\frac{D}{2}\ln(2\pi)-\frac{1}{2}\mathbb{E}_{\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k}[(\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T \boldsymbol{\Lambda}_k (\boldsymbol{x}_n-\boldsymbol{\mu}_k)] \end{eqnarray}$

ρの各項を求めるには π,μ,Λ の事後分布が必要 → 後回し

「q^*(Z) のモーメントは他の変数の事後分布によるモーメントに相互に依存」
「繰り返しで解く必要あり」

両辺の exponential を取ると
$q^*(\boldsymbol{Z}) = \frac{1}{C} \prod_{n=1}^N \prod_{k=1}^K \rho_{nk}^{z_{nk}}$

正規化定数 C は、Z_n が 1-of-K であることを用いると
$C = \sum_{\boldsymbol{Z}} \prod_{n=1}^N \prod_{k=1}^K \rho_{nk}^{z_{nk}} = \sum_{z_1}...\sum_{z_N} \prod_{n=1}^N \prod_{k=1}^K \rho_{nk}^{z_{nk}} = \prod_{n=1}^N \sum_{j=1}^K \rho_{nj}$

$\sum_{k=1} z_{nk} = 1$ より
$\sum_{j=1}^K \rho_{nj} = \left(\sum_{j=1}^K \rho_{nj}\right)^{\sum_{k=1} z_{nk}} = \prod_{k=1}^K \left(\sum_{j=1}^K \rho_{nj}\right)^{z_{nk}}$

したがって

$\begin{eqnarray}q^*(\boldsymbol{Z}) &=& \frac{1}{C} \prod_{n=1}^N \prod_{k=1}^K \rho_{nk}^{z_{nk}}\\&=& \prod_{n=1}^N \prod_{k=1}^K \left(\frac{\rho_{nk}}{\textstyle\sum_{j=1}^K \rho_{nj}}\right)^{z_{nk}}\\&=& \prod_{n=1}^N \prod_{k=1}^K r_{nk}^{z_{nk}}\end{eqnarray}$

where $r_{nk} = \frac{\rho_{nk}}{\textstyle\sum_{j=1}^K \rho_{nj}}$

このとき

$\mathbb{E}[z_{nk}] = \sum_{\boldsymbol{Z}} z_{nk}q^*(\boldsymbol{Z}) = r_{nk}$

が言えるので、r_nk は「負担率」を表していることがわかる

q(π) を推定する

$\begin{eqnarray}&&\ln q^*(\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})\\ &=& \mathbb{E}_{\boldsymbol{Z}}[\ln p(\boldsymbol{X},\boldsymbol{Z},\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda})] + \rm{const} \\&=& \mathbb{E}_{\boldsymbol{Z}}[\ln p(\boldsymbol{X}|\boldsymbol{Z},\boldsymbol{\mu},\boldsymbol{\Lambda})] + \mathbb{E}_{\boldsymbol{Z}}[\ln p(\boldsymbol{Z}|\boldsymbol{\pi})] +\ln p(\boldsymbol{\pi}) + \sum_{k=1}^K \ln p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) + \rm{cnst}\end{eqnarray}$

この右辺は π のみを含む項と μ_k、Λ_k のみを含む項に分かれることから、
q(π,μ,Λ) は以下の「導出される分解」(10.2.5) を持つことがわかる

$q(\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Lambda}) = q(\boldsymbol{\pi})\prod_{k=1}^K q(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k)$

したがってπに依存する項を取り出すことで q^*(π) を得る。

$\begin{eqnarray} \ln q^*(\boldsymbol{\pi}) &=& \mathbb{E}_{\boldsymbol{Z}}[\ln p(\boldsymbol{Z}|\boldsymbol{\pi})] +\ln p(\boldsymbol{\pi}) + \rm{const} \\ &=& \sum_{k=1}^K \sum_{n=1}^N r_{nk} \ln \pi_k + (\alpha_0-1)\sum_{k=1}^K \ln \pi_k + \rm{const} \\ &=& \sum_{k=1}^K (\alpha_0+N_k-1) \ln \pi_k + \rm{const} \end{eqnarray}$

where $N_k = \sum_{n=1}^N r_{nk}$

この式から q^*(π) がディリクレ分布であることがわかる。

$q^*(\boldsymbol{\pi})= \rm{Dir}(\boldsymbol{\pi}|\boldsymbol{\alpha})$ , where $\boldsymbol{\alpha} = (\alpha_k), \; \alpha_k = \alpha_0 + N_k$

q(μ|Λ) を推定する

同様に μ_k,Λ_k に依存する項を取り出すと、

$\begin{eqnarray} && \ln q^*(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) = \ln q^*(\boldsymbol{\mu}_k|\boldsymbol{\Lambda}_k)+ \ln q^*(\boldsymbol{\Lambda}_k) \\ &=& \mathbb{E}_{\boldsymbol{Z}}[\ln p(\boldsymbol{X}|\boldsymbol{Z},\boldsymbol{\mu},\boldsymbol{\Lambda})] + \ln p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) + \rm{const} \\ &=& \sum_{n=1}^N r_{nk} \ln \mathcal{N}(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k^{-1}) + \ln \mathcal{N}(\boldsymbol{\mu}_k|\boldsymbol{m}_0,(\beta\boldsymbol{\Lambda}_k)^{-1}) + \ln \mathcal{W}(\boldsymbol{\Lambda}_k|\boldsymbol{W}_0,\nu_0) + \rm{const} \\ &=& \sum_{n=1}^N \frac{r_{nk}}{2} \{ \ln|\boldsymbol{\Lambda}_k| - D \ln(2\pi) - (\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T \boldsymbol{\Lambda}_k (\boldsymbol{x}_n-\boldsymbol{\mu}_k)\} \\ && + \frac{1}{2} \{ \ln (\beta_0|\boldsymbol{\Lambda}_k|) - D \ln(2\pi) - \beta_0 (\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T \boldsymbol{\Lambda}_k (\boldsymbol{\mu}_k-\boldsymbol{m}_0)\} \\ && - \frac{\nu_0}{2}\ln|\boldsymbol{W}_0| - \frac{\nu_0 D}{2}\ln 2 - \frac{D(D-1)}{4}\ln \pi - \sum_{i=1}^D \ln \Gamma\left( \frac{\nu_0+1+i}{2} \right) \\ && \; + \frac{\nu_0-D-1}{2}\ln |\Lambda_k| - \frac{1}{2}\rm{Tr}(\boldsymbol{W}_0^{-1}\Lambda_k) + \rm{const} \end{eqnarray}$

μ_k を含む項のみ取り出すと、

$\begin{eqnarray} && \ln q^*(\boldsymbol{\mu}_k|\boldsymbol{\Lambda}_k) \\&=& - \sum_{n=1}^N \frac{r_{nk}}{2} (\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T \boldsymbol{\Lambda}_k (\boldsymbol{x}_n-\boldsymbol{\mu}_k) - \frac{\beta_0}{2} (\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T \boldsymbol{\Lambda}_k (\boldsymbol{\mu}_k-\boldsymbol{m}_0) + \rm{const} \\ &=& - \frac{1}{2} \boldsymbol{\mu}_k^T (N_k \boldsymbol{\Lambda}_k + \beta_0 \boldsymbol{\Lambda}_k) \boldsymbol{\mu}_k + \boldsymbol{\mu}_k^T( N_k \boldsymbol{\Lambda}_k \bar{\boldsymbol{x}}_k + \beta_0 \boldsymbol{\Lambda}_k \boldsymbol{m}_0) + \rm{const}\\ &=& - \frac{1}{2}(N_k + \beta_0) \boldsymbol{\mu}_k^T \boldsymbol{\Lambda}_k \boldsymbol{\mu}_k + \boldsymbol{\mu}_k^T \boldsymbol{\Lambda}_k( N_k \bar{\boldsymbol{x}}_k + \beta_0 \boldsymbol{m}_0) + \rm{const}\end{eqnarray}$

where $N_k = \sum_{n=1}^N r_{nk}, \; \bar{\boldsymbol{x}}_k=\frac{1}{N_k}\sum_{n=1}^N r_{nk} \boldsymbol{x}_n$

これは q^*(μ_k|Λ_k) の exponential に対応するので、q^* がガウス分布 $\mathcal{N}(\boldsymbol{\mu}_k|\boldsymbol{m}_k, (\beta_k \boldsymbol{\Lambda}_k)^{-1})$ であることがわかる。ここで

$\beta_k = N_k + \beta_0$ ,
$\boldsymbol{m}_k= \frac{N_k \bar{\boldsymbol{x}}_k + \beta_0 \boldsymbol{m}_0}{N_k + \beta_0}$

q(Λ) を推定する

$\ln q^*(\boldsymbol{\Lambda}_k) = \ln q^*(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) - \ln q^*(\boldsymbol{\mu}_k|\boldsymbol{\Lambda}_k)$ から Λ_k に関する項のみ取り出す。

$\begin{eqnarray} && \ln q^*(\boldsymbol{\Lambda}_k) \\ &=& \sum_{n=1}^N \frac{r_{nk}}{2} \{ \ln|\boldsymbol{\Lambda}_k| - (\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T \boldsymbol{\Lambda}_k (\boldsymbol{x}_n-\boldsymbol{\mu}_k)\} \\ && + \frac{1}{2} \{ \ln (\beta_0|\boldsymbol{\Lambda}_k|) - \beta_0 (\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T \boldsymbol{\Lambda}_k (\boldsymbol{\mu}_k-\boldsymbol{m}_0)\} \\ && + \frac{\nu_0-D-1}{2}\ln |\boldsymbol{\Lambda}_k| - \frac{1}{2}\rm{Tr}(\boldsymbol{W}_0^{-1}\boldsymbol{\Lambda}_k) \\ && - \frac{1}{2}\left\{ \ln (\beta_k|\boldsymbol{\Lambda}_k|) - \beta_k (\boldsymbol{\mu}_k-\boldsymbol{m}_k)^T \boldsymbol{\Lambda}_k (\boldsymbol{\mu}_k-\boldsymbol{m}_k)\} \right\} + \rm{const} \end{eqnarray}$

ここでベクトル x に対し、x^TΛx はスカラーゆえトレースをとっても同じ値なので、(C.8) とあわせて次のような変形ができる。
$\boldsymbol{x}^T\boldsymbol{\Lambda}\boldsymbol{x} = \rm{Tr}(\boldsymbol{x}^T\boldsymbol{\Lambda}\boldsymbol{x}) = \rm{Tr}(\boldsymbol{x}\boldsymbol{x}^T\boldsymbol{\Lambda})$

これを用いてさらに ln q^*(Λ_k) を変形しつつ、ln|Λ| と Tr(*Λ) に整理すると

$\begin{eqnarray}&& \ln q^*(\boldsymbol{\Lambda}_k) \\ &=&\{\boldsymbol{W}_0^{-1} + \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\boldsymbol{\mu}_k)(\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T+ \beta_0(\boldsymbol{\mu}_k-\boldsymbol{m}_0)(\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T - \beta_k (\boldsymbol{\mu}_k-\boldsymbol{m}_k)(\boldsymbol{\mu}_k-\boldsymbol{m}_k)^T\} \boldsymbol{\Lambda}_k \\ && + \frac{\textstyle\sum_{n=1}^N r_{nk} + \nu_0-D-1}{2} \ln|\boldsymbol{\Lambda}_k| + \rm{const}\end{eqnarray}$

よって q^*(Λ_k) はウィシャート分布 $\mathcal{W}(\boldsymbol{\Lambda}_k|\boldsymbol{W}_k,\nu_k)$ であることがわかる。ここで

$\boldsymbol{W}_k^{-1}=\boldsymbol{W}_0^{-1} + \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\boldsymbol{\mu}_k)(\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T+ \beta_0(\boldsymbol{\mu}_k-\boldsymbol{m}_0)(\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T - \beta_k (\boldsymbol{\mu}_k-\boldsymbol{m}_k)(\boldsymbol{\mu}_k-\boldsymbol{m}_k)^T$ ,
$\nu_k = \sum_{n=1}^N r_{nk} + \nu_0 = \nu_0 + N_k$

さらに W_k を変形していく。

$\begin{eqnarray}&& \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\boldsymbol{\mu}_k)(\boldsymbol{x}_n-\boldsymbol{\mu}_k)^T \\&=& \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k+\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k+\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)^T \\&=& \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)^T \\&& + \sum_{n=1}^N r_{nk}\{(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)^T +(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)^T +(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)^T \}\\ &=& \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)^T +N_k(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)^T \end{eqnarray}$

$\begin{eqnarray}&& N_k(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)(\bar{\boldsymbol{x}}_k-\boldsymbol{\mu}_k)^T+ \beta_0(\boldsymbol{\mu}_k-\boldsymbol{m}_0)(\boldsymbol{\mu}_k-\boldsymbol{m}_0)^T- \beta_k (\boldsymbol{\mu}_k-\boldsymbol{m}_k)(\boldsymbol{\mu}_k-\boldsymbol{m}_k)^T \\&=& (N_k+\beta_0-\beta_k)\boldsymbol{\mu}_k\boldsymbol{\mu}_k^T + (-N_k\bar{\boldsymbol{x}}_k-\beta_0\boldsymbol{m}_0+\beta_k\boldsymbol{m}_k)\boldsymbol{\mu}_k^T \\&&+ \boldsymbol{\mu}_k(-N_k\bar{\boldsymbol{x}}_k-\beta_0\boldsymbol{m}_0+\beta_k\boldsymbol{m}_k)^T + N_k\bar{\boldsymbol{x}}_k \bar{\boldsymbol{x}}_k^T + \beta_0\boldsymbol{m}_0\boldsymbol{m}_0^T - \beta_k\boldsymbol{m}_k\boldsymbol{m}_k\\&=& N_k\bar{\boldsymbol{x}}_k \bar{\boldsymbol{x}}_k^T + \beta_0\boldsymbol{m}_0\boldsymbol{m}_0^T - \beta_k\left(\frac{N_k \bar{\boldsymbol{x}}_k + \beta_0 \boldsymbol{m}_0}{\beta_k}\right)\left(\frac{N_k \bar{\boldsymbol{x}}_k + \beta_0 \boldsymbol{m}_0}{\beta_k}\right)^T\\&=&\beta_k^{-1}\{ (N_k\beta_k-N_k^2)\bar{\boldsymbol{x}}_k\bar{\boldsymbol{x}}_k^T -\beta_0N_k\boldsymbol{m}_0\bar{\boldsymbol{x}}_k^T -\beta_0N_k\bar{\boldsymbol{x}}_k\boldsymbol{m}_0^T + (\beta_0\beta_k-\beta_0^2)\boldsymbol{m}_0\boldsymbol{m}_0^T\}\\&=&\frac{\beta_0N_k}{\beta_k}(\bar{\boldsymbol{x}}_k-\boldsymbol{m}_0)(\bar{\boldsymbol{x}}_k-\boldsymbol{m}_0)^T\end{eqnarray}$

これらをあわせると

$\boldsymbol{W}_k^{-1}=\boldsymbol{W}_0^{-1} + \sum_{n=1}^N r_{nk}(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)(\boldsymbol{x}_n-\bar{\boldsymbol{x}}_k)^T + \frac{\beta_0N_k}{\beta_k}(\bar{\boldsymbol{x}}_k-\boldsymbol{m}_0)(\bar{\boldsymbol{x}}_k-\boldsymbol{m}_0)^T$

【(2) に続く】