混合正規分布

混合正規分布を定義します

MORE...

対数尤度関数

混合正規分布モデルの対数尤度関数を定義し、その最大化を考えます

MORE...

EM アルゴリズム

E-step と M-step を交互に反復する EM アルゴリズムにより、パラメータを最適化します

MORE...

混合正規分布モデル(GMM)とEMアルゴリズム

はじめに

混合正規分布モデル (GMM: Gaussian mixture model) は、凸凹した確率密度関数を複数の正規分布の和によって表現するもので、ここでは観測されたサンプル集合をもとに EM アルゴリズムで各正規分布のパラメータを推定します。

以下、太字はベクトルまたは行列を表し、行列 $\boldsymbol{A}$ の逆行列の転置を $\boldsymbol{A}^{-T}=\left(\boldsymbol{A}^{-1}\right)^T$ と略記することにします。

正規分布

$D$ 次元の正規分布は、平均ベクトル $\boldsymbol{\mu}\in\mathbb{R}^D$、非特異な分散共分散行列 $\boldsymbol{\Sigma}\in\mathbb{R}^{D\times D}$ をパラメータとして持ち、次のように表されます。 \begin{eqnarray} N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\} \label{Normal} \end{eqnarray}

混合正規分布

混合正規分布を、正規分布 $N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})$ の線形結合として、次のように定義します。 \begin{eqnarray} p(\boldsymbol{x}) &=& \sum_{k=1}^K \pi_k N(\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \label{GMM} \end{eqnarray} ただし、結合係数 $\pi_k$ は確率に対応するもので、 \begin{eqnarray} 0\leq \pi_k\leq 1 \end{eqnarray} かつ \begin{eqnarray} \sum_{k=1}^K \pi_k &=& 1 \label{sumpi1} \end{eqnarray} を満たすものとします。

混合正規分布モデルの対数尤度関数

観測されたサンプル $\boldsymbol{x}_n\in\mathbb{R}^D,\ n=1,2,3,\cdots N$ が式(\ref{GMM})の混合正規分布モデルから出てきたものと仮定し、その対数尤度関数 \begin{eqnarray} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& \log \prod_{n=1}^N p(\boldsymbol{x}_n) \\ &=& \sum_{n=1}^N \log p(\boldsymbol{x}_n) \\ &=& \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \end{eqnarray} を最大化 (正確には極大化) するパラメータ $\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}$ の組を推定します。 ここで $\boldsymbol{X}$ はサンプル集合 $\{\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots \boldsymbol{x}_N\}$ を行列で表したものです。 \begin{eqnarray} \boldsymbol{X} &=& \left( \begin{array}{ccc} & \boldsymbol{x}_1^T & \\ & \boldsymbol{x}_2^T & \\ & \vdots & \\ & \boldsymbol{x}_N^T & \\ \end{array} \right) \in \mathbb{R}^{N\times D} \end{eqnarray}

パラメータ推定の方針としては、対数尤度関数の極大点では、全ての $m=1,2,3,\cdots K$ で \begin{eqnarray} \left\{ \begin{array}{rcl} \displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0\\ \displaystyle\frac{\partial}{\partial\boldsymbol{\Sigma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0\\ \displaystyle\frac{\partial}{\partial\pi_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& 0 \end{array} \right.\nonumber \end{eqnarray} が成り立つことを利用します。

対数尤度関数を極大化する $\boldsymbol{\mu}_m$ を求める

対数尤度関数を $\boldsymbol{\mu}_m$ で偏微分すると \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& \frac{\partial}{\partial\boldsymbol{\mu}_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \\ &=& \sum_{n=1}^N \frac{\partial}{\partial\boldsymbol{\mu}_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) \\ && 一般に \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} ですから \nonumber\\ &=& \sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} \displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \\ && k\neq m の場合、\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0\ ですので \nonumber\\ &=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \label{dmu} \end{eqnarray} ここで式(\ref{Normal})から、正規分布 $N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)$ は \begin{eqnarray} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{x}_n-\boldsymbol{\mu}_m)\right\} \nonumber \\ && \boldsymbol{\mu}_m で偏微分することを考慮して \boldsymbol{x} と \boldsymbol{\mu} を入れ替え \nonumber\\ &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \end{eqnarray} これにより、式(\ref{dmu})の分子の $\pi_m$ の後ろは次のように書けます。 \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) &=& \frac{\partial}{\partial\boldsymbol{\mu}_m} (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\ &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \frac{\partial}{\partial\boldsymbol{\mu}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\ && 一般に \exp'(f(\boldsymbol{x}))=\exp(f(\boldsymbol{x}))f'(\boldsymbol{x})\ ですから \nonumber\\ &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}_m|^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\ && \hspace{32mm}\cdot\frac{\partial}{\partial\boldsymbol{\mu}_m} \left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\ &=& N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \frac{\partial}{\partial\boldsymbol{\mu}_m} \left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\ &=& -\frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \frac{\partial}{\partial\boldsymbol{\mu}_m} \left\{(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Sigma}_m^{-1}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \\ && 一般に \displaystyle\frac{\partial\boldsymbol{x}^T \boldsymbol{A} \boldsymbol{x}}{\partial\boldsymbol{x}}=(\boldsymbol{A}+\boldsymbol{A}^T)\boldsymbol{x}\ ですから \nonumber\\ &=& -\frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \left(\boldsymbol{\Sigma}_m^{-1}+\boldsymbol{\Sigma}_m^{-T}\right) (\boldsymbol{\mu}_m-\boldsymbol{x}_n) \\ && 分散共分散行列 \boldsymbol{\Sigma}_m は対称行列ですから \boldsymbol{\Sigma}_m^{-1} も対称行列なので \left(\boldsymbol{\Sigma}_m^{-1}+\boldsymbol{\Sigma}_m^{-T}\right)=2 \boldsymbol{\Sigma}_m^{-1} \nonumber\\ &=& - N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{\mu}_m-\boldsymbol{x}_n) \\ &=& N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m) \end{eqnarray} よって式(\ref{dmu})は \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\mu}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \nonumber\\ &=& \sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m) \\ \end{eqnarray} ここで \begin{eqnarray} \gamma_{n,m} &=& \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \label{gamma} \end{eqnarray} と置けば \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\mu}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} (\boldsymbol{x}_n-\boldsymbol{\mu}_m) \\ && 極大点では勾配が \boldsymbol{0} になりますので \nonumber\\ &=& \boldsymbol{0} \end{eqnarray} が成り立つには \begin{eqnarray} \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} \boldsymbol{x}_n &=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m^{-1} \boldsymbol{\mu}_m \end{eqnarray} 両辺に左から $\boldsymbol{\Sigma}_m$ を掛けて \begin{eqnarray} \sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n &=& \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\mu}_m \end{eqnarray} よって \begin{eqnarray} \boldsymbol{\mu}_m &=& \frac{\displaystyle\sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n}{\displaystyle\sum_{n=1}^N \gamma_{n,m}} \end{eqnarray}

負担率

式(\ref{gamma})の $\gamma_{n,m}$ は負担率と呼ばれるもので、サンプル $\boldsymbol{x}_n$ が $m$ 番目の正規分布 $N(\boldsymbol{x}|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)$ から出てきた確率を表します。

サンプル $\boldsymbol{x}_n$ は $K$ 個ある正規分布のどれかから出てきたわけですから、$m=1,2,3,\cdots K$ について $\gamma_{n,m}$ を足し合わせた全確率は 1 になります。 \begin{eqnarray} \require{cancel} \sum_{m=1}^K \gamma_{n,m} &=& \sum_{m=1}^K \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{i=1}^K \pi_i N(\boldsymbol{x}_n|\boldsymbol{\mu}_i,\boldsymbol{\Sigma}_i)} \\ &=& \frac{\cancel{\displaystyle\sum_{m=1}^K \pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}}{\cancel{\displaystyle\sum_{m=1}^K \pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}} \\ &=& 1 \label{sumgamma} \end{eqnarray}

対数尤度関数を極大化する $\boldsymbol{\Sigma}_m$ を求める

分散共分散行列 $\boldsymbol{\Sigma}_m$ の逆行列を \begin{eqnarray} \boldsymbol{\Gamma}_m=\boldsymbol{\Sigma}_m^{-1} \label{G} \end{eqnarray} と置くと \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Sigma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}_m)\ =\ \boldsymbol{0} \end{eqnarray} を満たす解は \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1})\ =\ \boldsymbol{0} \end{eqnarray} を解くことで計算できます。 \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1}) &=& \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1}) \\ &=& \sum_{n=1}^N \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1}) \\ && 一般に \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} ですから \nonumber\\ &=& \sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \\ && k\neq m の場合、\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0\ ですので \nonumber\\ &=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \label{dgamma} \end{eqnarray} 式(\ref{Normal})から、正規分布 $N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})$ は \begin{eqnarray} N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}) &=& (2\pi)^{-\frac{D}{2}} |\boldsymbol{\Sigma}|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\} \nonumber \end{eqnarray} ですから、式(\ref{dgamma})の分子の $\pi_m$ の後ろは次のように書けます。 \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1}) &=& (2\pi)^{-\frac{D}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m}|\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \label{dNdG} \end{eqnarray} 一般に $\{f(\boldsymbol{x})g(\boldsymbol{x})\}'=f'(\boldsymbol{x})g(\boldsymbol{x})+f(\boldsymbol{x})g'(\boldsymbol{x})$ であることから、$\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m}$ 以降は次のように書けます。 \begin{eqnarray} && \left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\right) \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} + |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \left[\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \right] \label{ddGm} \end{eqnarray} 前半の (  ) 内は
\begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} &=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{-\frac{1}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m| \\ && 一般に \frac{\partial|\boldsymbol{X}|}{\partial\boldsymbol{X}}=|\boldsymbol{X}|\boldsymbol{X}^{-T} ですので \nonumber\\ &=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{-\frac{1}{2}} |\boldsymbol{\boldsymbol{\Gamma}_m}|\boldsymbol{\boldsymbol{\Gamma}_m}^{-T} \\ &=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\boldsymbol{\boldsymbol{\Gamma}_m}^{-T} \\ \end{eqnarray} 後半の [  ] 内は
  $\displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\}$ \begin{eqnarray} && 一般に \exp'(f(\boldsymbol{x}))=\exp(f(\boldsymbol{x}))f'(\boldsymbol{x}) ですので \nonumber\\ &=& \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \label{dexpdGm} \\ && 一般に \displaystyle\frac{\partial\boldsymbol{a}^T \boldsymbol{X}^{-1} \boldsymbol{a}}{\partial\boldsymbol{X}}=\boldsymbol{a}\boldsymbol{a}^T ですから \nonumber\\ &=& -\frac{1}{2} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \end{eqnarray} よって式(\ref{ddGm})は \begin{eqnarray} && \hspace{-15mm}\left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\right) \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} + |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \left(\frac{\partial}{\partial\boldsymbol{\Gamma}_m} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \right) \\ &=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}}\boldsymbol{\boldsymbol{\Gamma}_m}^{-T} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\ && \hspace{6mm}-\frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \\ && 第1項,第2項共通に含まれる |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} と \exp\{\quad\} 部分は、スカラーなので括り出せて \nonumber\\ &=& \frac{1}{2} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \left\{ \boldsymbol{\boldsymbol{\Gamma}_m}^{-T} - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \end{eqnarray} 式(\ref{dNdG})は、$(2\pi)^{-\frac{D}{2}}$ を式(\ref{ddGm})に掛けたものですから \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1}) &=& (2\pi)^{-\frac{D}{2}} \frac{\partial}{\partial\boldsymbol{\Gamma}_m}|\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\} \nonumber\\ && 式(\ref{ddGm})を代入して \nonumber\\ &=& \frac{1}{2} \underbrace{(2\pi)^{-\frac{D}{2}} |\boldsymbol{\Gamma}_m|^{\frac{1}{2}} \exp\left\{-\frac{1}{2}(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T\boldsymbol{\Gamma}_m(\boldsymbol{\mu}_m-\boldsymbol{x}_n)\right\}}_{= N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)} \left\{ \boldsymbol{\boldsymbol{\Gamma}_m}^{-T} - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \\ &=& \frac{1}{2} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) \left\{ \boldsymbol{\boldsymbol{\Gamma}_m}^{-T} - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \label{ddGN} \end{eqnarray} よって、式(\ref{dgamma})は \begin{eqnarray} \frac{\partial}{\partial\boldsymbol{\Gamma}_m} \log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Gamma}_m^{-1}) &=& \sum_{n=1}^N \frac{\pi_m \displaystyle\frac{\partial}{\partial\boldsymbol{\Gamma}_m} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Gamma}_m^{-1})}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \nonumber\\ && 式(\ref{ddGN})を代入して \nonumber\\ &=& \frac{1}{2} \sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Gamma}_k^{-1})} \left\{ \boldsymbol{\boldsymbol{\Gamma}_m}^{-T} - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \\ && \boldsymbol{\Gamma}_m は、分散共分散行列 \boldsymbol{\Sigma}_m の逆行列でしたから \boldsymbol{\Gamma}_m^{-T}=\boldsymbol{\Sigma}_m^T \nonumber\\ &=& \frac{1}{2} \sum_{n=1}^N \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \left\{ \boldsymbol{\Sigma}_m^T - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \\ && 式(\ref{gamma}) の \gamma_{n,m} を使って書けば \nonumber\\ &=& \frac{1}{2} \sum_{n=1}^N \gamma_{n,m} \left\{ \boldsymbol{\Sigma}_m^T - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \\ && \boldsymbol{\Sigma}_m は対称行列ですから \boldsymbol{\Sigma}_m^T=\boldsymbol{\Sigma}_m なので \nonumber\\ &=& \frac{1}{2} \sum_{n=1}^N \gamma_{n,m} \left\{ \boldsymbol{\Sigma}_m - (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \right\} \\ && 極大点では勾配が \boldsymbol{0} になりますので \nonumber\\ &=& \boldsymbol{0} \end{eqnarray} とすると \begin{eqnarray} \sum_{n=1}^N \gamma_{n,m} \boldsymbol{\Sigma}_m &=& \sum_{n=1}^N \gamma_{n,m} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \end{eqnarray} より、以下を得ます。 \begin{eqnarray} \boldsymbol{\Sigma}_m &=& \frac{1}{\displaystyle\sum_{n=1}^N \gamma_{n,m}} \sum_{n=1}^N \gamma_{n,m} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T \end{eqnarray}

対数尤度関数を極大化する $\pi_m$ を求める

混合係数 $\pi_m$ は $\displaystyle\sum_{m=1}^K \pi_m=1$ を満たさねばならないという強い制約があるため、Lagrange の未定乗数法により、対数尤度関数に $\lambda\left(\displaystyle\sum_{j=1}^K \pi_k -1\right)$ を加えたものを極大化することにします。 \begin{eqnarray} \frac{\partial}{\partial\pi_m} \left\{\log p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}) + \lambda\left(\sum_{j=1}^K \pi_k -1\right)\right\} &=& \left\{\frac{\partial}{\partial\pi_m} \sum_{n=1}^N \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}+\lambda \\ &=& \left\{\sum_{n=1}^N \frac{\partial}{\partial\pi_m} \log \sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}+\lambda \\ && 一般に \log'(f(\boldsymbol{x}))=\frac{f'(\boldsymbol{x})}{f(\boldsymbol{x})} ですから \nonumber\\ &=& \sum_{n=1}^N \frac{\displaystyle\frac{\partial}{\partial\pi_m} \displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} +\lambda \\ && k\neq m の場合、\frac{\partial}{\partial\pi_m} \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)=0\ ですので \nonumber\\ &=& \sum_{n=1}^N \frac{N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)}+\lambda \\ && 式(\ref{gamma}) の \gamma_{n,m} を使って書けば \nonumber\\ &=& \frac{1}{\pi_m} \sum_{n=1}^N \gamma_{n,m}+\lambda \\ &=& 0 \end{eqnarray} とすると、両辺に $\pi_m$ を掛けて \begin{eqnarray} \sum_{n=1}^N \gamma_{n,m} + \pi_m \lambda &=& 0 \end{eqnarray} より \begin{eqnarray} \pi_m &=& -\frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} \label{pim} \end{eqnarray} ですが、式(\ref{sumpi1})より \begin{eqnarray} \sum_{m=1}^K \pi_m &=& 1 \end{eqnarray} であることを考えると \begin{eqnarray} \sum_{m=1}^K \pi_m &=& - \sum_{m=1}^K \frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} &=& 1 \end{eqnarray} でなければならず、変形してゆくと \begin{eqnarray} - \sum_{m=1}^K \frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} &=& - \frac{1}{\lambda} \sum_{m=1}^K \sum_{n=1}^N \gamma_{n,m} \\ &=& - \frac{1}{\lambda} \sum_{n=1}^N \sum_{m=1}^K \gamma_{n,m} \\ && 式(\ref{sumgamma}) より \sum_{m=1}^K \gamma_{n,m}=1 ですから \nonumber\\ &=& - \frac{1}{\lambda} \sum_{n=1}^N 1 \\ &=& - \frac{N}{\lambda} \\ &=& 1 \end{eqnarray} から \begin{eqnarray} \lambda &=& -N \end{eqnarray} 従って式(\ref{pim})は \begin{eqnarray} \pi_m &=& -\frac{1}{\lambda} \sum_{n=1}^N \gamma_{n,m} \nonumber\\ &=& \frac{1}{N} \sum_{n=1}^N \gamma_{n,m} \nonumber\\ \end{eqnarray} となります。

EM アルゴリズム

これまでに得た結果をまとめると \begin{eqnarray} \gamma_{n,m} &=& \frac{\pi_m N(\boldsymbol{x}_n|\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)}{\displaystyle\sum_{k=1}^K \pi_k N(\boldsymbol{x}_n|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)} \nonumber \end{eqnarray} の下で対数尤度関数を極大化するパラメータは \begin{eqnarray} \left\{ \begin{array}{rcl} \boldsymbol{\mu}_m &=& \frac{\displaystyle\sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n}{\displaystyle\sum_{n=1}^N \gamma_{n,m}}, \\ \boldsymbol{\Sigma}_m &=& \displaystyle\frac{1}{\displaystyle\sum_{n=1}^N \gamma_{n,m}} \displaystyle\sum_{n=1}^N \gamma_{n,m} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T, \\ \pi_m &=& \displaystyle\frac{1}{N} \sum_{n=1}^N \gamma_{n,m} \end{array} \right.\nonumber \end{eqnarray} ですが、$\displaystyle\sum_{n=1}^N \gamma_{n,m}$ は $m$ 番目の正規分布から出てきたサンプル数の推定値に相当しますので、これを $N_m$ と書くことにすると、もう少し簡単になります。 \begin{eqnarray} \left\{ \begin{array}{rcl} \boldsymbol{\mu}_m &=& \displaystyle\frac{1}{N_m} \sum_{n=1}^N \gamma_{n,m} \boldsymbol{x}_n,\\ \boldsymbol{\Sigma}_m &=& \displaystyle\frac{1}{N_m} \sum_{n=1}^N \gamma_{n,m} (\boldsymbol{\mu}_m-\boldsymbol{x}_n)(\boldsymbol{\mu}_m-\boldsymbol{x}_n)^T,\\ \pi_m &=& \displaystyle\frac{N_m}{N} \end{array} \right. \end{eqnarray} しかし $\gamma_{n,m}$ の計算式の中に $\boldsymbol{\mu}_m, \boldsymbol{\Sigma}_m, \pi_m$ が、 $\boldsymbol{\mu}_m, \boldsymbol{\Sigma}_m, \pi_m$ の計算式の中に $\gamma_{n,m}$ があるため、この連立方程式は非線形でスパッと解くことができません。 そこで次のように初期値から出発して、E-step と M-step を交互に反復する EM アルゴリズムにより、ジワジワと最適解に近づけてゆくようにします。

  1. 初期値 $\boldsymbol{\mu}_m^{(0)}, \boldsymbol{\Sigma}_m^{(0)}, \pi_m^{(0)}$ を適当に設定し、$i=0$ とする。
  2. $\boldsymbol{\mu}_m^{(i)}, \boldsymbol{\Sigma}_m^{(i)}, \pi_m^{(i)}$ を使い、$\gamma_{n,m}^{(i)}, N_m^{(i)}$ を計算する (E-step)。 \begin{eqnarray} \left\{ \begin{array}{rcl} \gamma_{n,m}^{(i)} &=& \frac{\pi_m^{(i)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_m^{(i)},\boldsymbol{\Sigma}_m^{(i)})}{\displaystyle\sum_{k=1}^K \pi_k^{(i)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k^{(i)},\boldsymbol{\Sigma}_k^{(i)})} \\ N_m^{(i)} &=& \displaystyle\sum_{n=1}^N \gamma_{n,m}^{(i)} \end{array} \right. \end{eqnarray}
  3. $\gamma_{n,m}^{(i)}, N_m^{(i)}$ を使い、パラメータ $\boldsymbol{\mu}_m^{(i+1)}, \boldsymbol{\Sigma}_m^{(i+1)}, \pi_m^{(i+1)}$ を更新する (M-step)。 \begin{eqnarray} \left\{ \begin{array}{rcl} \boldsymbol{\mu}_m^{(i+1)} &=& \displaystyle\frac{1}{N_m^{(i)}} \sum_{n=1}^N \gamma_{n,m}^{(i)} \boldsymbol{x}_n, \\ \boldsymbol{\Sigma}_m^{(i+1)} &=& \displaystyle\frac{1}{N_m^{(i)}} \sum_{n=1}^N \gamma_{n,m}^{(i)} (\boldsymbol{\mu}_m^{(i)}-\boldsymbol{x}_n)(\boldsymbol{\mu}_m^{(i)}-\boldsymbol{x}_n)^T, \\ \pi_m^{(i+1)} &=& \displaystyle\frac{N_m^{(i)}}{N} \end{array} \right. \end{eqnarray}
  4. パラメータ $\boldsymbol{\mu}_m^{(i+1)}, \boldsymbol{\Sigma}_m^{(i+1)}, \pi_m^{(i+1)}$ で対数尤度関数 \begin{eqnarray} \log p(\boldsymbol{X}|\boldsymbol{\pi}^{(i+1)},\boldsymbol{\mu}^{(i+1)},\boldsymbol{\Sigma}^{(i+1)}) &=& \sum_{n=1}^N \log \sum_{k=1}^K \pi_k^{(i+1)} N(\boldsymbol{x}_n|\boldsymbol{\mu}_k^{(i+1)},\boldsymbol{\Sigma}_k^{(i+1)}) \end{eqnarray} を計算し、十分収束したら終了。さもなくば $i\leftarrow i+1$ としてステップ 2 に戻る。
【注意】EM アルゴリズムは広域最適解に収束するとは限りません。