t分布のパラメータをEMアルゴリズムで推定する

t分布のパラメータをEMアルゴリズムで推定するメモ．

計算する値が多いのでまとめておく．平均$\vec{\mu}$, 分散 $\mathbf{\Sigma}$ は解析的に求まる．ただし，自由度$\nu$は求まらないのでグリッドサーチ，ランダムサーチ，他の反復法などが必要になる．

Estep:
- 各データに対して，事後分布になるガンマ分布$q_{i}(h_{i})$を求める．
Mstep:
- 平均パラメータ$\vec{\mu}$と共分散行列パラメータ $\mathbf{\Sigma}$ の推定
  - 事後分布の期待値$E[ q_{i}(h_{i}) ]$を使う．
- 自由度$\nu$の推定
  - 条件付期待値$\sum_{i=1}^{I} \int \hat{q}_{i}(h_{i}) \log{\left[ Pr(\vec{x}_{i}, h_{i}|\mathbf{\theta}) \right] } d h_{i}$を計算（グリッドサーチなどで自由度を変えながらこの期待値が大きくなっているか確認）

参考：

混合分布とt分布

潜在変数と混合分布の考え方を使うとt分布を導出できる．

混合正規分(MoG)では，潜在変数がcategorical分布（試行回数1の多項分布）に従うとし，重みをcategorical分布の確率とした混合分布であった．

t分布におけるEMアルゴリズム

EMアルゴリズムの概要

通常のMLでは，データに対して対数尤度を最大化するときのパラメータを求める．

$\hat{\vec{\theta}} = \mathop{ \text{argmax}}_{\vec{\theta}} \left[ \sum_{i=1}^{I} \log{ \left[ Pr(\vec{x}_{i} | \vec{\theta}) \right] } \right]$

$Pr(x)$ を容易に計算できない場合，EMアルゴリズムでは潜在変数を導入して同時分布を周辺化して$Pr(x)$を求め尤度の最大化を図る．

$\hat{\vec{\theta}} = \mathop{\text{argmax}}_{\vec{\theta}} \left[ \sum_{i=1}^{I} \log{ \left[ \int Pr(\vec{x}_{i}, \vec{h}_{i} | \vec{\theta}) d\vec{h}_{i} \right] } \right]$

実際には上式でなく，次の下界を最大化する．

$\begin{eqnarray} \cal{B}\left[ \left\{ q_{i}(\vec{h}_{i}) \right\}, \vec{\theta} \right] &=& \sum_{i=1}^{I} \int q_{i}(\vec{h}_{i}) \log{ \left[ \frac{Pr(\vec{x}_{i}, \vec{h}_{i} | \vec{\theta})}{q_{i}(\vec{h}_{i})} \right] }d\vec{h}_{i} \\\\ &\leq& \sum_{i=1}^{I} \log{ \left[ \int q_{i}(\vec{h_{i}}) \frac{Pr(\vec{x}_{i}, \vec{h}_{i} | \vec{\theta})}{ q_{i}(\vec{h_{i}}) } d\vec{h}_{i} \right] } = \sum_{i=1}^{I} \log{ \left[ \int Pr(\vec{x}_{i}, \vec{h}_{i} | \vec{\theta}) d\vec{h}_{i} \right] } \end{eqnarray}$

GOAL: t分布のパラメータ $\mathbf{\theta} = \left{ \vec{\mu}, \mathbf{\Sigma}, \nu \right}$ をデータ$\mathbf{x}_{1 \cdots I}$から推定する．

下界最大化にはE-stepとM-stepを使って反復させる．

E-step: 前回のM-stepで推定したパラメータ$\mathbf{\theta}$で固定して，$q(h_{i})$を変えて下界を最大化する．（＝潜在変数$h$の事後分布$q(h_{i})$を求める．）

$q_{i}(h_{i}) = Pr(h_{i} | \vec{x}_{i}, \mathbf{\theta}^{[ t ] }) = \frac{Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}^{[ t ] } ) Pr(h_{i})}{ Pr(\vec{x}_{i} | \mathbf{\theta}^{ [ t ] })}$

M-step: E-stepで求めた事後分布$q_{i}(h_{i})$で固定して，データと潜在変数の同時分布密度の条件付期待値（下界）を最大化するときのパラメータ$\mathbf{\theta}$を求める．

$\hat{\mathbf{\theta}}^{[ t+1 ] } = \mathop{argmax}_{\mathbf{\theta}} \sum_{ i=1 }^{ I } \int \hat{ q }_{ i }( h_{ i } ) \log{\left[ Pr(\vec{x}_{i}, h_{i} | \mathbf{\theta}) \right] } d h_{i}$

E-Mstepのイメージ

E-step

E-stepで求める事後分布 $q_{i}(h_{i})$ がどんな式になるか導出する．

結果としては以下のガンマ分布になる．これは，事前分布をガンマ分布，尤度を正規分布としたとき”分散の逆数”（精度）にかかるスケールパラメータに対して共役であるため．

$\begin{eqnarray} q_{i}( h_{i} ) = Pr( h_{i} | \vec{x}_{i}, \mathbf{\theta}^{[ t ] }) &=& \frac{Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}^{ [ t ] })Pr( h_{i} )}{ Pr(\vec{x}_{i} | \mathbf{\theta}^{ [ t ] })} \\\\ &=& \frac{ \text{Norm}_{\vec{x}_{i}} \left[ \vec{\mu}, \frac{\mathbf{\Sigma}}{ h_{i}} \right] \text{Gam}_{h_{i}}\left[ \frac{\nu}{2}, \frac{\nu}{2}\right] }{Pr(\vec{x}_{i})} \\\\ &=& \text{Gam}_{h_{i}}\left[ \frac{\nu+D}{2}, \frac{(\vec{x}_{i} - \mu)^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})}{2} + \frac{\nu}{2} \right] \end{eqnarray}$

導出過程は，以下に載せる．

分子の式展開：

$\begin{eqnarray} \text{Norm}_{\vec{x}_{i}} \left[ \vec{\mu}, \frac{\mathbf{\Sigma}}{ h_{i}} \right] \text{Gam}_{h_{i}}\left[ \frac{\nu}{2}, \frac{\nu}{2}\right] &=& \frac{1}{(2\pi)^{\frac{D}{2}} |\frac{\mathbf{\Sigma}}{h_{i}}|^{\frac{1}{2}} }\exp{\left[ - \frac{1}{2}(\vec{x}_{i}-\vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right]} \cdot \frac{ (\frac{\nu}{2})^{\frac{\nu}{2}}}{\Gamma\left[ \frac{\nu}{2} \right] } \exp{\left[ -\frac{\nu}{2} h_{i} \right]} h_{i}^{\frac{\nu}{2} - 1} \end{eqnarray}$

分母の式展開：分母の周辺化はt分布になっている．

$\begin{eqnarray} Pr(\vec{x}_{i}) &=& \int Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}^{ [ t ] }) Pr(h_{i}) dh_{i} \\\\ &=& \text{Stud}_{x}\left[ \mu, \sigma^{2}, \nu \right] \\\\ &=& \frac{\Gamma\left[ \frac{\nu + D}{2} \right] }{(\nu \pi)^{\frac{D}{2}} |\Sigma|^{\frac{1}{2}} \Gamma\left[ \frac{\nu}{2} \right] } \left( 1 + \frac{(\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1} (\vec{x}_{i} - \vec{\mu})}{\nu}\right)^{-\frac{\nu + D}{2}} \end{eqnarray}$

分子に出てくる定数倍の行列式は，$|a A| = a^{D} |A|$，このように展開できるので

$|\frac{\mathbf{\Sigma}}{h_{i}}|^{\frac{1}{2}} = (\frac{1}{h_{i}^{D}} |\mathbf{\Sigma}|)^{\frac{1}{2}} = \frac{1}{h_{i}^{\frac{D}{2}}} |\mathbf{\Sigma}|^{\frac{1}{2}}$

$\begin{eqnarray} \frac{Pr(\vec{x}_{i}|h_{i})Pr(h_{i})}{Pr(\vec{x}_{i})} &=& \frac{ \frac{1}{(2\pi)^{\frac{D}{2}} \frac{1}{h_{i}^{\frac{D}{2}}}|\mathbf{\Sigma}|^{\frac{1}{2}} }\exp{\left[ - \frac{1}{2}(\vec{x}_{i}-\vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right] } \cdot \frac{ (\frac{\nu}{2})^{\frac{\nu}{2}}}{\Gamma\left[ \frac{\nu}{2} \right] } \exp\left[ -\frac{\nu}{2} h_{i} \right] h_{i}^{\frac{\nu}{2} - 1} }{ \frac{\Gamma\left[ \frac{\nu + D}{2} \right] }{(\nu \pi)^{\frac{D}{2}} |\Sigma|^{\frac{1}{2}} \Gamma\left[ \frac{\nu}{2} \right] } \left( 1 + \frac{(\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1} (\vec{x}_{i} - \vec{\mu})}{\nu}\right)^{-\frac{\nu + D}{2}} } \\\\ &=& \frac{ \frac{1}{2^{\frac{D}{2}} }\exp{\left[ - \frac{1}{2}(\vec{x}_{i}-\vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right]} \cdot \frac{ (\frac{\nu}{2})^{\frac{\nu}{2}}}{1} \exp\left[ -\frac{\nu}{2} h_{i} \right] h_{i}^{\frac{\nu+D}{2} - 1} } { \frac{\Gamma\left[ \frac{\nu + D}{2} \right] }{\nu^{\frac{D}{2}} } \left( 1 + \frac{(\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1} (\vec{x}_{i} - \vec{\mu})}{\nu}\right)^{-\frac{\nu + D}{2}} } \\\\ &=& \frac{(\frac{\nu}{2})^{\frac{\nu+D}{2}} \left( 1 + \frac{(\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1} (\vec{x}_{i} - \vec{\mu})}{\nu}\right)^{\frac{\nu + D}{2}} \exp{ \left[ - \frac{1}{2} \left\{ \nu + (\vec{x}_{i}-\vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right\} h_{i} \right] } h_{i}^{\frac{\nu+D}{2} - 1} } { \Gamma\left[ \frac{\nu + D}{2} \right] } \\\\ &=& \frac{ \left\{ \frac{1}{2} \left( \nu + (\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1} (\vec{x}_{i} - \vec{\mu})\right)\right\}^{\frac{\nu + D}{2}} \exp{ \left[ - \frac{1}{2} \left\{ \nu + (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right\} h_{i} \right] } h_{i}^{\frac{\nu+D}{2} - 1} } { \Gamma\left[ \frac{\nu + D}{2} \right] } \\\\ &=& \text{Gamma}_{h_{i}}\left[ \frac{\nu + D}{2}, \frac{1}{2} \left\{ \nu + (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}( \vec{x}_{i} - \vec{\mu} ) \right\}\right] \end{eqnarray}$

M-Step

E-stepの事後分布で固定して条件付期待値の対数（下界）を最大にするときのパラメータ$\mathbf{\theta}$を求める．

式展開すると以下の期待値になる．

$\begin{eqnarray} \hat{\mathbf{\theta}}^{[ t+1 ] } &=& \mathop{\text{argmax}}_{\mathbf{\theta}} \left[ \sum_{i=1}^{I} \int \hat{q}_{i}(h_{i}) \log{\left[ Pr(\vec{x}_{i}, h_{i}|\mathbf{\theta}) \right] } d h_{i} \right] \\\\ &=& \mathop{\text{argmax}}_{\mathbf{\theta}} \left[ \sum_{i=1}^{I} \int \hat{q}_{i}(h_{i}) \left( \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } + \log{\left[ Pr(h_{i}) \right] }\right) dh_{i} \right] \\\\ &=& \mathop{\text{argmax}}_{\mathbf{\theta}} \left[ \sum_{i=1}^{I} \int Pr(h_{i} | \vec{x}_{i}, \mathbf{\theta}^{[ t ] }) \left( \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } + \log{\left[ Pr(h_{i}) \right] }\right) dh_{i} \right] \\\\ &=& \mathop{\text{argmax}}_{\mathbf{\theta}} \left[ \sum_{i=1}^{I}\left\{ E\left[ \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } \right] + E\left[ \log{\left[ Pr( h_{ i } ) \right] } \right] \right\} \right] \\ \end{eqnarray}$

この2つの期待値がどんな式になるのか分かればいい．

ガンマ分布における対数変換した正規分布密度の期待値$E[\log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] }]$

$\begin{eqnarray} \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } &=& \log{\left[ \text{Norm}_{\vec{x}_{i}}(\vec{\mu}, \mathbf{\Sigma}/h_{i}) \right] }\\ &=& - \frac{D}{2} \log{2\pi} - \frac{1}{2}\log{|\mathbf{\Sigma}/h_{i}|} - \frac{1}{2} (\vec{x}- \vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}- \vec{\mu}) \\\\ &=& - \frac{D}{2} \log{2\pi} - \frac{1}{2}\log{(|\mathbf{\Sigma}|\frac{1}{h_{i}^{D}})} - \frac{1}{2} (\vec{x}- \vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}- \vec{\mu}) \\\\ &=& -\frac{D}{2}\log{2\pi} - \frac{1}{2}\log{|\mathbf{\Sigma}|} + \frac{D}{2}\log{h_{i}} - \frac{1}{2} (\vec{x}- \vec{\mu})^{T}h_{i}\mathbf{\Sigma}^{-1}(\vec{x}- \vec{\mu}) \end{eqnarray}$

よって，期待値を取ると，

（xは観測データで，今は定数である．）

$\begin{eqnarray} E[ \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } ] &=& - \frac{D}{2}\log{2\pi} - \frac{1}{2} \log{|\mathbf{\Sigma}|} + \frac{D}{2} E[ \log{h_{i}} ] - \frac{1}{2} (\vec{x}- \vec{\mu})^{T} \mathbf{\Sigma}^{-1}(\vec{x} - \vec{\mu})E[ h_{i} ] \\\\ &=& \frac{D \cdot E[ \log{h_{i}} ] - D \log{2\pi} - \log{|\mathbf{\Sigma}|} - ( \vec{x} - \vec{\mu} )^{T} \mathbf{\Sigma}^{-1}( \vec{x}- \vec{\mu} ) E[ h_{i} ] }{2} \end{eqnarray}$

ガンマ分布における事前ガンマ分布密度の対数の期待値 $E[\log{Pr(h_{i})} ]$

$\begin{eqnarray} \log{Pr(h_{i})} &=& \log{ \text{Gam}_{h_{i}}\left[ \frac{\nu}{2}, \frac{\nu}{2} \right]} \\\\ &=& \log{ \left[ \frac{(\frac{\nu}{2})^{\frac{\nu}{2}}}{\Gamma\left[ \frac{\nu}{2} \right] } \exp{ \left[ -\frac{\nu}{2} h_{i} \right] } h_{i}^{\frac{\nu}{2} - 1} \right] } \\\\ &=& \frac{\nu}{2}\log{[ \frac{\nu}{2} ] } - \log{\Gamma\left[ \frac{\nu}{2} \right]} - \frac{\nu}{2} h_{i} + (\frac{\nu}{2} - 1)\log{h_{i}} \end{eqnarray}$

よって，期待値は，

$E[ \log{Pr(h_{i})} ] = \frac{\nu}{2}\log{[ \frac{\nu}{2} ] } - \log{\Gamma\left[ \frac{\nu}{2} \right] } - \frac{\nu}{2} E[ h_{ i } ] + (\frac{\nu}{2} - 1) E[ \log{h_{i}} ]$

ここで，条件付期待値Eを求めるには厳密には$E[h_{i} ], E[ \log{h_{i}} ]$を求める必要がある．

ガンマ分布の期待値$E[ h_{i} ]$を求める

まず，一般的なガンマ分布$\text{Gamma}_{x}(\alpha, \beta)$の期待値を導出する．

参考：https://to-kei.net/distribution/gamma-distribution/g-parameter-derivation/

$\begin{eqnarray} E[ x ] &=& \int_{0}^{\infty} x \text{Gam}_{x}\left[ \alpha, \beta \right] dx \\\\ &=& \int_{0}^{\infty} x \frac{\beta^{\alpha}}{\Gamma\left[ \alpha \right] } \exp\left[ -\beta x \right] x^{\alpha - 1} dx \\\\ &=& \frac{1}{\Gamma\left[ \alpha \right] } \int_{0}^{\infty} \beta^{\alpha}x^{\alpha} \exp\left[ -\beta x \right] dx \\\\ &=& \frac{1}{\Gamma\left[ \alpha \right] } \int_{0}^{\infty} (\beta x)^{\alpha} \exp\left[ -\beta x \right] dx \\\\ &=& \frac{\beta^{-1}}{\Gamma\left[ \alpha \right] } \int_{0}^{\infty} (\beta x)^{\alpha} \exp\left[ -\beta x \right] \beta dx \end{eqnarray}$

$\beta x = t$と置いて変数変換する. $dt = \beta dx$

$E[ x ] = \frac{\beta^{-1}}{\Gamma\left[ \alpha \right] } \int_{0}^{\infty} t^{\alpha} \exp\left[ -t \right] dt$

この積分はガンマ積分になっている．

$\begin{eqnarray} E[ x ] &=& \frac{\beta^{-1}}{\Gamma\left[ \alpha \right] } \Gamma[ \alpha + 1 ] \\ &=& \frac{\beta^{-1}}{\Gamma\left[ \alpha \right] } \alpha \Gamma[ \alpha ] \\ &=& \frac{\alpha}{\beta} \end{eqnarray}$

t分布のパラメータ推定に話を戻すと，$E[h_{i} ]$はどうなるか．

E-stepで求まる事後分布は以下なので，

$\text{Gamma}_{ h_{ i } }\left[ \frac{\nu + D}{2}, \frac{1}{2} \left\{ \nu + (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right\}\right]$

事後分布の期待値は，以下になる．

$E[ h_{i} ] = \frac{\nu + D}{\nu + (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})}$

ガンマ分布の確率変数を対数変換した期待値$E[ \log{ h_{i} } ]$を求める

参考：https://math.stackexchange.com/questions/138252/expected-value-of-ln-x-if-x-is-gammaa-b-distributed

まず，一般的なガンマ分布$\text{Gamma}_{x}(\alpha, \beta)$で，対数変換した変数の期待値を求める．

$\begin{eqnarray} E[ \log{x} ] &=& \int_{0}^{\infty} \log{x} \text{Gam}_{x}\left[ \alpha, \beta \right] dx \\\\ &=& \int_{0}^{\infty} \log{x} \frac{\beta^{\alpha}}{\Gamma\left[ \alpha \right]} \exp{\left[ -\beta x \right]} x^{\alpha - 1} dx \end{eqnarray}$

$\beta x = y$として

$\begin{eqnarray} E[ \log{x} ] &=& \frac{1}{\Gamma\left[ \alpha \right]} \int_{0}^{\infty} \log{x} {(\beta x)}^{\alpha - 1} \exp\left[ -\beta x \right] \beta dx \\\\ &=& \frac{1}{\Gamma\left[ \alpha \right]} \int_{0}^{\infty} \log{(\frac{y}{\beta})} y^{\alpha - 1} \exp{ \left[ -y \right]} dy \\\\ &=& \frac{1}{\Gamma\left[ \alpha \right]} \int_{0}^{\infty} (\log{y} - \log{(\beta)}) y^{\alpha - 1} \exp{\left[ -y \right]} dy \\\\ &=& -\log{\beta} + \frac{1}{\Gamma\left[ \alpha \right] } \int_{0}^{\infty} \log{(y)} y^{\alpha - 1} \exp\left[ -y \right] dy \end{eqnarray}$

ここで，ガンマ関数の微分を考える．

$$ \Gamma\left[ \alpha \right] = \int_{0}^{\infty} x^{\alpha-1} e^{-x} dx $$

$\begin{eqnarray} \Gamma[ \alpha ]' = \frac{d \Gamma\left[ \alpha \right] }{d \alpha} &=& \int_{0}^{\infty} \frac{d}{d \alpha} \exp{(\log{x^{\alpha-1}} )} e^{-x} dx \\\\ &=& \int_{0}^{\infty} \frac{d}{d \alpha} \exp{\{ (\alpha-1)(\log{x} ) \}} e^{-x } dx \\\\ &=& \int_{0}^{\infty} \exp{\{ (\alpha-1)(\log{x} ) \}} \log{(x)} e^{-x } dx \\\\ &=& \int_{0}^{\infty} \log{(x)} x^{\alpha-1} e^{-x } dx \end{eqnarray}$

$\begin{eqnarray} E[ \log{x} ] &=& \frac{\Gamma\left[ \alpha \right]' }{\Gamma\left[ \alpha \right] } - \log{\beta} \end{eqnarray}$

$\frac{d \log{ f( x ) } }{ d x } = \frac{ f'(x) }{ f(x) }$より，

$\begin{eqnarray} E[ \log{x} ] &=& \frac{d \log{ \Gamma[ \alpha ] }}{d \alpha} - \log{\beta} \end{eqnarray}$

$\frac{d \log{ \Gamma[ \alpha ] }}{d \alpha}$ をdigmma関数$\psi$ とする．つまり，対数の期待値は以下になる．

$\begin{eqnarray} E[ \log{x} ] = \frac{d \log{ \Gamma[\alpha ] }}{d \alpha} - \log{\beta} = \psi(\alpha) - \log{\beta} \end{eqnarray}$

よって，t分布のパラメータ推定に話を戻すと，事後分布のガンマ分布は次より，

$\text{Gamma}_{h_{i}}\left[ \frac{\nu + D}{2}, \frac{1}{2} \left\{ \nu + (\vec{x}_{i}-\vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu}) \right\}\right]$

各パラメータ$\mathbf{\theta}$を推定

M-stepの条件付き期待値をt分布の各パラメータで微分して0とおけば推定できる．ただし，自由度$\nu$は解析的に求められない．

結果，以下の式になる．

$\begin{eqnarray} \vec{\mu}^{[ t+1 ] } &=& \frac{\sum_{i=1}^{ I } E[ h_{i} ] \vec{x}_{i}}{\sum_{i=1}^{I} E[ h_{i} ]} \\\\ \mathbf{\Sigma}^{[ t+1 ] } &=& \frac{\sum_{i=1}^{ I } E[ h_{i} ] (\vec{x}_{i} - \vec{\mu}^{[ t+1 ] }) (\vec{x}_{i} - \vec{\mu}^{ [ t+1 ] } )^{T}}{\sum_{i=1}^{ I } E[ h_{i} ] } \end{eqnarray}$

平均パラメータ$\mu$の更新式

M-stepの条件付き期待値をパラメータ$\mu$で微分すると，

$\frac{d E}{d \vec{\mu}} = \frac{d}{d \vec{\mu} } \sum_{i=1}^{I} \frac{D \cdot E[ \log{h_{i}} ] - D\log{2 \pi} - \log{|\mathbf{\Sigma}|} - (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})E[ h_{i} ] }{2} = 0 \\\\ \sum_{i=1}^{I} - \frac{1}{2} (- 2 \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})E[ h_{i} ] ) = 0\\\\ \sum_{i=1}^{I} \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})E[ h_{i} ] = 0 \\\\ \sum_{i=1}^{I} (\vec{x}_{i} - \vec{\mu})E[ h_{i} ] = 0 \\\\ \sum_{i=1}^{I} \vec{x}_{i} E[ h_{i} ] - \vec{\mu}\sum_{i=1}^{I} E[ h_{i} ] = 0\\ \vec{\mu}\sum_{i=1}^{I} E[ h_{i} ] = \sum_{i=1}^{I} \vec{x}_{i} E[ h_{i} ] \\\\ \vec{\mu} = \frac{\sum_{i=1}^{I} \vec{x}_{i} E[ h_{i} ] }{\sum_{i=1}^{I} E[ h_{i} ] }$

共分散行列パラメータ $\mathbf{\Sigma}$ の更新式

M-stepの条件付き期待値をパラメータ $\mathbf{\Sigma}$ で微分すると，

$\frac{d \sum_{i=1}^{I} E[\log{ Pr(h_{i})} ] }{d \mathbf{\Sigma}} = 0, \frac{d \log{|A|} }{d A} = (A^{-1})^{T}, \frac{d \vec{x}^{T} \mathbf{A} \vec{x}}{d \mathbf{A}} = - (\mathbf{A}^{-1})^{T}\vec{x}\vec{x}^{T}(\mathbf{A}^{-1})^{T}$

より，

$\begin{eqnarray} \frac{d \sum_{i=1}^{I} E[\log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] ] }}{d \mathbf{\Sigma}} = \frac{d }{d \mathbf{\Sigma}} \sum_{i=1}^{I} \left(-\frac{1}{2}\log{|\mathbf{\Sigma}|} -\frac{1}{2}(\vec{x}_{i} - \vec{\mu})^{T} \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})E[ h_{i} ] \right) &=& 0 \\\\ \sum_{i=1}^{I} \left\{ -\frac{1}{2}(\mathbf{\Sigma}^{-1})^{T} - \frac{1}{2}\left(- \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})(\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1} E[ h_{i} ] \right) \right\} &=& 0 \end{eqnarray}$

逆行列のみの項を左辺に移行

$\begin{eqnarray} \sum_{i=1}^{I} \mathbf{\Sigma}^{-1} &=& \sum_{i=1}^{I} \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})(\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1} E[ h_{i} ] \\\\ I \mathbf{\Sigma}^{-1} &=& \sum_{i=1}^{I} \mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})(\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1} E[ h_{i} ] \\\\ I \mathbf{\Sigma} &=& \sum_{i=1}^{I} (\vec{x}_{i} - \vec{\mu})(\vec{x}_{i} - \vec{\mu})^{T} E[ h_{i} ] \\\\ \mathbf{\Sigma} &=& \frac{\sum_{i=1}^{I} (\vec{x}_{i} - \vec{\mu})(\vec{x}_{i} - \vec{\mu})^{T} E[ h_{i} ] }{I} \\\\ \end{eqnarray}$

1次元の正規分布データからt分布のパラメータを推定するcode

各ステップの式をまとめる．

E-Step：

$\begin{eqnarray} q_{i}(h_{i}) = Pr(h_{i} | \vec{x}_{i}, \mathbf{\theta}^{[ t ]}) &=& \frac{Pr( \vec{x}_{i} | h_{i}, \mathbf{\theta}^{[ t ] })Pr(h_{i})}{Pr(\vec{x}_{i} | \mathbf{\theta}^{[ t ]})} \\\\ &=& \frac{ \text{Norm}_{x_{i}} \left[ \mu, \frac{\sigma^{2}}{ h_{i}} \right] \text{Gam}_{h_{i}}\left[ \frac{\nu}{2}, \frac{\nu}{2} \right] }{Pr(\vec{x}_{i})} \\\\ &=& \text{Gam}_{h_{i}}\left[ \frac{\nu+1}{2}, \frac{(x_{i} - \mu)^{2}(\sigma^{2})^{-1}}{2} + \frac{\nu}{2} \right] \\ \end{eqnarray}$

M-step:

$\mathop{\text{argmax}}_{\mathbf{\theta}} \left[ \sum_{i=1}^{I}\left\{ E\left[ \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } \right] + E\left[ \log{\left[ Pr(h_{i}) \right] } \right] \right\} \right]$

$\mu = \frac{\sum_{i=1}^{I} x_{i} E[ h_{i} ] }{\sum_{i=1}^{I} E[ h_{i} ] } \\\\ \sigma^{2} = \frac{\sum_{i=1}^{I} (x_{i} - \mu)^{2} E[ h_{i} ] }{I}$

$\begin{eqnarray} E\left[ \log{\left[ Pr(\vec{x}_{i} | h_{i}, \mathbf{\theta}) \right] } \right] &=& \frac{D \cdot E[ \log{h_{i}} ] - D\log{2 \pi} - \log{|\mathbf{\Sigma}|} - (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})E[ h_{i} ] }{2} \\\\ &=& \frac{1 \cdot E[ \log{h_{i}} ] - 1 \cdot \log{2 \pi} - \log{\sigma^{2}} - (x_{i} - \mu)^{2}(\sigma^{2})^{-1}E[ h_{i} ] }{2} \end{eqnarray}$

$\begin{eqnarray} E\left[ \log{\left[ Pr(h_{i}) \right] } \right] &=& \frac{\nu}{2}\log{\left[ \frac{\nu}{2} \right] } - \log{\Gamma\left[ \frac{\nu}{2}\right] } + \left( \frac{\nu}{2} - 1\right) E[ \log{h_{i}} ] - \frac{\nu}{2}E[ h_{i} ] \end{eqnarray}$

$\begin{eqnarray} E[ \log{h_{i}} ] &=& \psi \left[ \frac{\nu +D}{2} \right] - \log{\left[ \frac{\nu + (\vec{x}_{i} - \vec{\mu})^{T}\mathbf{\Sigma}^{-1}(\vec{x}_{i} - \vec{\mu})}{2} \right] } \\\\ &=& \psi \left[ \frac{\nu + 1}{2} \right] - \log{\left[ \frac{\nu + (x_{i} - \mu)^{2}(\sigma^{2})^{-1}}{2} \right] } \end{eqnarray}$

真の分布を標準正規分布$N(0, 1)$とし，
そこから観測データを100個サンプリングする．
初期分布を$t(-2, 2, \nu = 1)$からEMアルゴリズムを使ってパラメータを推定する．
自由度$\nu$はグリッドサーチで変えながら条件付き期待値が最大になるパラメータを求める．

import scipy as sp
from scipy import stats
from scipy import special
import matplotlib.pyplot as plt

%matplotlib inline

def EStep(data_x, mu, sigma2, nu):
    '''
    return posterior Gamma each data
    '''
    q_h = []
    alpha = (nu+1)/2
    for i in sp.arange(len(data_x)):
        beta = ((data_x[i] - mu)**2/sigma2 + nu)*.5
        q_h.append(stats.gamma(a=alpha, scale=1/beta))
    return q_h


def calc_cond_expectation(data_x, expectations_posterior, mu, sigma2, nu):
    '''
    E[hi] = expectation_posterior
    E[log hi] = digamma[(nu+D)/2] - log[  ]
    '''

    E_hi = expectations_posterior
    E_log_hi = special.digamma((nu+1) / 2) - \
        sp.log((nu + ((data_x - mu)**2)/sigma2) * .5)

    E_posterior_logNormpdf = (E_log_hi - sp.log(2 * sp.pi) -
                              sp.log(sigma2) - (data_x - mu)**2 / sigma2 * E_hi) / 2
    E_posterior_priorGamma = nu * .5 * sp.log(nu * .5) - sp.log(
        special.gamma(nu * .5)) + (.5 * nu - 1) * E_log_hi - nu * .5 * E_hi

    return sp.sum(E_posterior_logNormpdf + E_posterior_priorGamma)


def MStep(data_x, q_h, nu):
    '''
    return
    '''
    I = len(data_x)
    expectations_posterior = sp.empty_like(q_h)
    for i in sp.arange(len(q_h)):
        dist = q_h[i]
        expectations_posterior[i] = dist.mean()

    mu = sp.sum(data_x * expectations_posterior) / sp.sum(expectations_posterior)
    sigma2 = sp.sum((data_x - mu)**2 * expectations_posterior) / I

    cond_expectation = calc_cond_expectation(
        data_x, expectations_posterior, mu, sigma2, nu)
    return cond_expectation, mu, sigma2, nu


n = 100
EM_iterations = 15
nus = sp.arange(1, n+5, 5)

truth_dist = stats.norm(loc=0, scale=sp.sqrt(1))
data_x = truth_dist.rvs(size=n)
# ini_para
mu, sigma2 = -2., 2.

# EM
max_para = {"cond_expectation": -1e20, "mu": mu, "sigma2": sigma2}

for i in sp.arange(EM_iterations):
    mu, sigma2 = max_para["mu"], max_para["sigma2"]

    for nu in nus:
        q_h = EStep(data_x, mu, sigma2, nu)
        c_e, mu, sigma2, nu = MStep(data_x, q_h, nu)

        if max_para["cond_expectation"] < c_e:
            max_para = {
                "cond_expectation": c_e,
                "mu": mu,
                "sigma2": sigma2,
                "nu": nu
            }
            print("i:", i, end=" ")
            print("max:", max_para)


plt.hist(data_x, bins=20, density=True)

x = sp.linspace(-7, 7, 100)
t_dist = stats.t(loc=max_para["mu"], scale=sp.sqrt(max_para["sigma2"]), df=max_para["nu"])
plt.plot(x, t_dist.pdf(x))
plt.tight_layout()
plt.show()

i: 0 max: {'cond_expectation': -262.2832259124987, 'mu': -0.534967953023884, 'sigma2': 0.821448924989282, 'nu': 1}
i: 0 max: {'cond_expectation': -201.08196951409755, 'mu': -0.02662036249833737, 'sigma2': 0.8928237170918831, 'nu': 6}
i: 0 max: {'cond_expectation': -189.87233960822854, 'mu': 0.11100223434262527, 'sigma2': 1.0040720996170245, 'nu': 11}
i: 0 max: {'cond_expectation': -178.58292159707602, 'mu': 0.1302403218320351, 'sigma2': 1.0637091848820932, 'nu': 16}
i: 0 max: {'cond_expectation': -168.18094371837404, 'mu': 0.13009946345391527, 'sigma2': 1.0943402299966092, 'nu': 21}
i: 0 max: {'cond_expectation': -159.1255565412171, 'mu': 0.1286775435488971, 'sigma2': 1.1123900510434954, 'nu': 26}
i: 0 max: {'cond_expectation': -151.29765752987493, 'mu': 0.12763401158513304, 'sigma2': 1.1244249220900722, 'nu': 31}
i: 0 max: {'cond_expectation': -144.46941473104542, 'mu': 0.12689199040076732, 'sigma2': 1.1331130219771408, 'nu': 36}
i: 0 max: {'cond_expectation': -138.435638740986, 'mu': 0.12633866659507098, 'sigma2': 1.1397130086509437, 'nu': 41}
i: 0 max: {'cond_expectation': -133.03832544436656, 'mu': 0.12590917369284998, 'sigma2': 1.1449086549309175, 'nu': 46}
i: 0 max: {'cond_expectation': -128.15924392689732, 'mu': 0.12556560605155995, 'sigma2': 1.1491098777985396, 'nu': 51}
i: 0 max: {'cond_expectation': -123.70919358531602, 'mu': 0.12528427748478907, 'sigma2': 1.152579494972972, 'nu': 56}
i: 0 max: {'cond_expectation': -119.61970292853879, 'mu': 0.12504955453819885, 'sigma2': 1.1554944861073908, 'nu': 61}
i: 0 max: {'cond_expectation': -115.83724776098929, 'mu': 0.1248506723402878, 'sigma2': 1.15797866234179, 'nu': 66}
i: 0 max: {'cond_expectation': -112.31927348944568, 'mu': 0.12467996330386541, 'sigma2': 1.1601213558893626, 'nu': 71}
i: 0 max: {'cond_expectation': -109.03142710599255, 'mu': 0.12453181180519048, 'sigma2': 1.1619886979279859, 'nu': 76}
i: 0 max: {'cond_expectation': -105.94559681225746, 'mu': 0.12440200768289753, 'sigma2': 1.16363072102759, 'nu': 81}
i: 0 max: {'cond_expectation': -103.03849772306813, 'mu': 0.12428733074330173, 'sigma2': 1.1650859930446327, 'nu': 86}
i: 0 max: {'cond_expectation': -100.2906332539778, 'mu': 0.12418527498002875, 'sigma2': 1.166384732376521, 'nu': 91}
i: 0 max: {'cond_expectation': -97.68551969630089, 'mu': 0.12409386041163888, 'sigma2': 1.1675509569732803, 'nu': 96}
i: 0 max: {'cond_expectation': -95.20909843640226, 'mu': 0.12401150159472255, 'sigma2': 1.1686040005059892, 'nu': 101}