続・わかりやすいパターン認識ひとり読書会 10 - 渡邊哲司（Satoshi Watanabe）の日記

今回は演習問題12.5を解いていきます。

続パタ読書会 1:ディリクレ過程混合モデル
 続パタ読書会 2:ディリクレ過程混合モデルのアルゴリズム
 続パタ読書会 3:実験
 続パタ読書会 4:実装
 続パタ読書会 5:事前確率P(s)の妥当性 1
続パタ読書会 6:事前確率P(s)の妥当性 2
続パタ読書会 7:ベル数 1
続パタ読書会 8:ベル数 2
続パタ読書会 9:P(s_k|x_k, s_-k, θ)の導出
続パタ読書会 10:演習問題12.5の計算
続パタ読書会 11:演習問題12.6の計算 1
続パタ読書会 12:演習問題12.6の計算 2

演習問題12.5の計算

演習問題12.5はテキスト中の(12.35)を導出せよとの問題です。具体的には

$\displaystyle\int{p(x_k\vert\theta_{new})p(\theta_{new})d\theta_{new}}=\left(\displaystyle\frac{\beta}{(1+\beta)\pi}\right)^{\frac{d}{2}}\cdot\displaystyle\frac{\vert \mathbf{S}_b\vert^{\frac{\nu+1}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1}{2}\right)}{\vert\mathbf{S}\vert^{\frac{\nu}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1-d}{2}\right)}$

ただし $\mathbf{S}_b$ は以下

$\mathbf{S}_b^{-1}=\mathbf{S}^{-1}+\displaystyle\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t$

を導出せよというものです。テキストでは上式の $p(\theta_{new})$ は $G_0(\theta_{new})$ となっています。

導出

以下の導出では積分消去される予定の $\theta_{new}$ は簡単のため $\theta$ で表します。初めに

$\begin{eqnarray} p(x_k\vert\theta) &=& \mathcal{N}(x_k;\mu, \mathbf{\Lambda}^{-1})\\ p(\theta) &=& \mathcal{N}(\mu;\mu_0, (\beta\mathbf{\Lambda})^{-1})\mathcal{W}(\mathbf{\Lambda};\nu, \mathbf{S}) \end{eqnarray}$

であったことを思い出しておきます。つまりパラメータの確率モデルを平均 $\mu$ 、分散共分散行列 $\mathbf{\Sigma}$ （式上では精度行列で表している）が共に未知の正規分布とし、その事前分布 $p(\theta)$ を共役事前分布の正規ウィシャート分布に設定していた場合の解が初めに提示していたものとなります。ではゴリゴリ計算していきます。まずは積分の中身と外側に出せる部分とに分割するまで計算していきます。

$\begin{eqnarray} \displaystyle\int{p(x_k\vert\theta)p(\theta)d\theta} &=& \displaystyle\int{\mathcal{N}(x_k;\mu, \mathbf{\Lambda}^{-1})\mathcal{N}(\mu;\mu_0, (\beta\mathbf{\Lambda})^{-1})\mathcal{W}(\mathbf{\Lambda};\nu, \mathbf{S})d\theta}\\ &=& \displaystyle\int\displaystyle\frac{|\mathbf{\Lambda}|^\frac{1}{2}}{(2\pi)^\frac{d}{2}}\mathrm{exp}\left(-\frac{1}{2}(x_k-\mu)^t\mathbf{\Lambda}(x_k-\mu)\right)\\ && \cdot\displaystyle\frac{|\beta\mathbf{\Lambda}|^\frac{1}{2}}{(2\pi)^\frac{d}{2}}\mathrm{exp}\left(-\frac{1}{2}(\mu-\mu_0)^t(\beta\mathbf{\Lambda})(\mu-\mu_0)\right)\\ && \cdot\displaystyle\frac{|\mathbf{\Lambda}|^{\frac{\nu-d-1}{2}}}{2^{\frac{\nu d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}|^{\frac{\nu}{2}}\cdot\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\mathrm{exp}\left(-\frac{1}{2}\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\right)d\theta\\ &=& \displaystyle\frac{1}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{\beta^{\frac{d}{2}}}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{1}{2^{\frac{\nu d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}|^{\frac{\nu}{2}}\cdot\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ && \cdot\displaystyle\int{Q(\theta)d\theta} \end{eqnarray}$

ここで $Q(\theta)$ は

$\begin{eqnarray} Q(\theta) &=& |\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{\nu-d-1}{2}}\\ && \cdot\mathrm{exp}\left(-\frac{1}{2}\left((x_k-\mu)^t\mathbf{\Lambda}(x_k-\mu)+(\mu-\mu_0)^t(\beta\mathbf{\Lambda})(\mu-\mu_0)+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\right)\right) \end{eqnarray}$

です。ここから積分の中身であるこの $Q(\theta)$ に絞って計算をしていきます。共役事前分布として正規ウィシャート分布を想定しているため、この $Q(\theta)$ も正規ウィシャート型になるであろうとわかっています。ですのでそれを見越して計算していくと見通しが良くなります。そこでまずは $\mathrm{exp(\cdot)}$ の外側にある $\mathbf{\Lambda}$ について正規分布型とウィシャート分布型に分かれるようにまとめてしまうと

$\begin{eqnarray} Q(\theta) &=& |\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{\nu-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}R(\theta)\right)\\ &=& |\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{(\nu+1)-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}R(\theta)\right)\\ &=& |\mathbf{\Lambda}|^{\frac{1}{2}}|\mathbf{\Lambda}|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}R(\theta)\right) \end{eqnarray}$

となります。ここで $R(\theta)$ は

$\begin{eqnarray} R(\theta) &=& (x_k-\mu)^t\mathbf{\Lambda}(x_k-\mu)+(\mu-\mu_0)^t(\beta\mathbf{\Lambda})(\mu-\mu_0)+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda}) \end{eqnarray}$

です。また $\nu_b=\nu+1$ と置きました。この $R(\theta)$ もやはり正規分布型の項とウィシャート分布型の項の和にまとまるように意識しながら計算していきます。ここで $\mathrm{exp}$ なので積ではなく和の形にすることに注意をします。

$\begin{eqnarray} R(\theta) &=& (x_k-\mu)^t\mathbf{\Lambda}(x_k-\mu)+(\mu-\mu_0)^t(\beta\mathbf{\Lambda})(\mu-\mu_0)+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& x_k^t\mathbf{\Lambda}x_k-x_k^t\mathbf{\Lambda}\mu-\mu^t\mathbf{\Lambda}x_k+\mu^t\mathbf{\Lambda}\mu\\ && +\beta\mu_0^t\mathbf{\Lambda}\mu_0-\beta\mu_0^t\mathbf{\Lambda}\mu-\mu^t\mathbf{\Lambda}\beta\mu_0+\beta\mu^t\mathbf{\Lambda}\mu\\ && +\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& (1+\beta)\mu^t\mathbf{\Lambda}\mu-\mu^t\mathbf{\Lambda}(x_k+\beta\mu_0)-(x_k+\beta\mu_0)^t\mathbf{\Lambda}\mu\\ && +x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& (1+\beta)\left(\mu^t\mathbf{\Lambda}\mu-\mu^t\mathbf{\Lambda}\left(\frac{x_k+\beta\mu_0}{1+\beta}\right)-\left(\frac{x_k+\beta\mu_0}{1+\beta}\right)^t\mathbf{\Lambda}\mu\right)\\ && +x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& (1+\beta)\left(\left(\mu-\frac{x_k+\beta\mu_0}{1+\beta}\right)^t\mathbf{\Lambda}\left(\mu-\frac{x_k+\beta\mu_0}{1+\beta}\right)\right)\\ && -\frac{1}{1+\beta}(x_k+\beta\mu_0)^t\mathbf{\Lambda}(x_k+\beta\mu_0)\\ && +x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& \left(\left(\mu-\frac{x_k+\beta\mu_0}{1+\beta}\right)^t\left((1+\beta)\mathbf{\Lambda}\right)\left(\mu-\frac{x_k+\beta\mu_0}{1+\beta}\right)\right)\\ && -\frac{1}{1+\beta}(x_k+\beta\mu_0)^t\mathbf{\Lambda}(x_k+\beta\mu_0)\\ && +x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& N(\theta)+W(\theta) \end{eqnarray}$

ここで

$\begin{eqnarray} N(\theta) &=& \left(\mu-\mu_b\right)^t\left((1+\beta)\mathbf{\Lambda}\right)\left(\mu-\mu_b\right)\\ W(\theta) &=& -\frac{1}{1+\beta}(x_k+\beta\mu_0)^t\mathbf{\Lambda}(x_k+\beta\mu_0)+x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda}) \end{eqnarray}$

と置きました。ただし

$\mu_b=\displaystyle\frac{x_k+\beta\mu_0}{1+\beta}$

です。ここまでですでに正規分布型 $N(\theta)$ は形が整いましたので $W(\theta)$ に絞ってウィシャート分布型を目指して計算していきます。

$\begin{eqnarray} W(\theta) &=& -\frac{1}{1+\beta}(x_k+\beta\mu_0)^t\mathbf{\Lambda}(x_k+\beta\mu_0)+x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0\\ &&+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& -\frac{1}{1+\beta}\left(x_k^t\mathbf{\Lambda}x_k+x_k^t\mathbf{\Lambda}\beta\mu_0+\beta\mu_0^t\mathbf{\Lambda}x_k+\beta^2\mu_0^t\mathbf{\Lambda}\mu_0\right)+x_k^t\mathbf{\Lambda}x_k+\beta\mu_0^t\mathbf{\Lambda}\mu_0\\ &&+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& \frac{\beta}{1+\beta}\left(\left(\frac{1+\beta}{\beta}-\frac{1}{\beta}\right)x_k^t\mathbf{\Lambda}x_k-x_k^t\mathbf{\Lambda}\mu_0-\mu_0^t\mathbf{\Lambda}x_k+\left((1+\beta)-\beta\right)\mu_0^t\mathbf{\Lambda}\mu_0\right)\\ &&+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& \frac{\beta}{1+\beta}\left(x_k^t\mathbf{\Lambda}x_k-x_k^t\mathbf{\Lambda}\mu_0-\mu_0^t\mathbf{\Lambda}x_k+\mu_0^t\mathbf{\Lambda}\mu_0\right)\\ &&+\mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})\\ &=& \mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})+\frac{\beta}{1+\beta}(x_k-\mu_0)^t\mathbf{\Lambda}(x_k-\mu_0) \end{eqnarray}$

ここで最後の式第二項の $\displaystyle\frac{\beta}{1+\beta}(x_k-\mu_0)^t\mathbf{\Lambda}(x_k-\mu_0)$ はスカラーのため $\mathrm{tr}$ でも同じ値なので

$\begin{eqnarray} W(\theta) &=& \mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})+\mathrm{tr}\left(\frac{\beta}{1+\beta}(x_k-\mu_0)^t\mathbf{\Lambda}(x_k-\mu_0)\right)\\ &=& \mathrm{tr}(\mathbf{S}^{-1}\mathbf{\Lambda})+\mathrm{tr}\left(\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t\mathbf{\Lambda}\right) && \because\mathrm{tr}(ABC)=\mathrm{tr}(CAB)\\ &=& \mathrm{tr}\left(\mathbf{S}^{-1}\mathbf{\Lambda}+\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t\mathbf{\Lambda}\right) && \because\mathrm{tr}(A)+\mathrm{tr}(B)=\mathrm{tr}(A+B)\\ &=& \mathrm{tr}\left(\left(\mathbf{S}^{-1}+\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t\right)\mathbf{\Lambda}\right)\\ &=& \mathrm{tr}(\mathbf{S}_b^{-1}\mathbf{\Lambda}) \end{eqnarray}$

となります。ただし

$\mathbf{S}_b^{-1}=\mathbf{S}^{-1}+\displaystyle\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t$

です。これでウィシャート分布型もまとまりましたので $Q(\theta)$ 全体をまとめると

$\begin{eqnarray} Q(\theta) &=& |\Lambda|^{\frac{1}{2}}|\Lambda|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}\left(N(\theta)+W(\theta)\right)\right)\\ &=& |\Lambda|^{\frac{1}{2}}\mathrm{exp}\left(-\frac{1}{2}N(\theta)\right)\cdot|\Lambda|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}W(\theta)\right)\\ &=& |\Lambda|^{\frac{1}{2}}\mathrm{exp}\left(-\frac{1}{2}\left(\mu-\mu_b\right)^t\left((1+\beta)\mathbf{\Lambda}\right)\left(\mu-\mu_b\right)\right)\\ && \cdot|\Lambda|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}\mathrm{tr}(\mathbf{S}_b^{-1}\mathbf{\Lambda})\right)\\ &=& \frac{1}{(1+\beta)^{\frac{d}{2}}}|(1+\beta)\Lambda|^{\frac{1}{2}}\mathrm{exp}\left(-\frac{1}{2}\left(\mu-\mu_b\right)^t\left((1+\beta)\mathbf{\Lambda}\right)\left(\mu-\mu_b\right)\right)\\ && \cdot|\Lambda|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}\mathrm{tr}(\mathbf{S}_b^{-1}\mathbf{\Lambda})\right)\\ \end{eqnarray}$

となります。よってこの $Q(\theta)$ の積分は

$\begin{eqnarray} \displaystyle\int{Q(\theta)d\theta} &=& \frac{1}{(1+\beta)^{\frac{d}{2}}}\displaystyle\int|(1+\beta)\Lambda|^{\frac{1}{2}}\mathrm{exp}\left(-\frac{1}{2}\left(\mu-\mu_b\right)^t\left((1+\beta)\mathbf{\Lambda}\right)\left(\mu-\mu_b\right)\right)\\ && \cdot|\Lambda|^{\frac{\nu_b-d-1}{2}}\mathrm{exp}\left(-\frac{1}{2}\mathrm{tr}(\mathbf{S}_b^{-1}\mathbf{\Lambda})\right)d\theta \end{eqnarray}$

となります。ここで積分の中を見ると $\mathcal{N}(\mu;\mu_b, ((1+\beta)\mathbf{\Lambda})^{-1})\cdot\mathcal{W}(\mathbf{\Lambda};\nu_b, \mathbf{S}_b)$ の形をしているので積分計算をするまでもなく正規化項の逆数となり、

$\displaystyle\int{Q(\theta)d\theta}=\frac{1}{(1+\beta)^{\frac{d}{2}}}\cdot(2\pi)^{\frac{d}{2}}\cdot 2^{\frac{\nu_b d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}_b|^{\frac{\nu_b}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu_b+1-j}{2}\right)}$

と計算できます。結果的に初めの式は

$\begin{eqnarray} \displaystyle\int{p(x_k\vert\theta)p(\theta)d\theta} &=& \displaystyle\frac{1}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{\beta^{\frac{d}{2}}}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{1}{2^{\frac{\nu d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}|^{\frac{\nu}{2}}\cdot\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ && \cdot\displaystyle\int{Q(\theta)d\theta}\\ &=& \displaystyle\frac{1}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{\beta^{\frac{d}{2}}}{(2\pi)^{\frac{d}{2}}}\cdot\displaystyle\frac{1}{2^{\frac{\nu d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}|^{\frac{\nu}{2}}\cdot\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ && \cdot\frac{1}{(1+\beta)^{\frac{d}{2}}}\cdot(2\pi)^{\frac{d}{2}}\cdot 2^{\frac{\nu_b d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}\cdot|\mathbf{S}_b|^{\frac{\nu_b}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu_b+1-j}{2}\right)}\\ &=& \frac{(2\pi)^{\frac{d}{2}}\cdot 2^{\frac{(\nu+1)d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}}{(2\pi)^{\frac{d}{2}}\cdot(2\pi)^{\frac{d}{2}}\cdot 2^{\frac{\nu d}{2}}\cdot\pi^{\frac{d(d-1)}{4}}}\cdot\left(\frac{\beta}{1+\beta}\right)^{\frac{d}{2}}\\ && \cdot\frac{|\mathbf{S}_b|^{\frac{\nu+1}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+2-j}{2}\right)}}{|\mathbf{S}|^{\frac{\nu}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ &=& \frac{1}{\pi^{\frac{d}{2}}}\cdot\left(\frac{\beta}{1+\beta}\right)^{\frac{d}{2}}\cdot\frac{|\mathbf{S}_b|^{\frac{\nu+1}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+2-j}{2}\right)}}{|\mathbf{S}|^{\frac{\nu}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ &=& \left(\frac{\beta}{(1+\beta)\pi}\right)^{\frac{d}{2}}\cdot\frac{|\mathbf{S}_b|^{\frac{\nu+1}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+2-j}{2}\right)}}{|\mathbf{S}|^{\frac{\nu}{2}}\displaystyle\prod_{j=1}^{d}{\Gamma\left(\frac{\nu+1-j}{2}\right)}}\\ &=& \left(\frac{\beta}{(1+\beta)\pi}\right)^{\frac{d}{2}}\cdot\frac{|\mathbf{S}_b|^{\frac{\nu+1}{2}}\Gamma\left(\displaystyle\frac{\nu+1}{2}\right)}{|\mathbf{S}|^{\frac{\nu}{2}}\Gamma\left(\displaystyle\frac{\nu+1-d}{2}\right)} \end{eqnarray}$

となり、導出が出来ました。

今回はこのくらいにして次回は演習問題12.6の導出をしようと思います。