この記事では、PRMLの第9章で述べられている、混合ガウスモデルに対するEMアルゴリズムをpythonで実装します。
対応するjupyter notebookは筆者のgithubリポジトリにあります。連載全体の方針や、PRMLの他のアルゴリズムの実装については、連載のまとめページをご覧いただければと思います。
1 理論のおさらい
この節では、混合ガウスモデルに対するEMアルゴリズムの実装に必要な結果だけを、簡潔にまとめます。導出や、EMアルゴリズムの一般的な説明1などは、本の方をご参照ください。
なお、EMアルゴリズムの収束性の証明などはPRMLには記載はなく、私もフォローできていません2。
1.1 設定
今回は教師なし学習を考えます。
- $N \in \mathbb{N}$ : データ点の個数
- $d \in \mathbb{N}$ : データの次元
- $x_0, x_1, \dots , x_{N-1} \in \mathbb{R}^d$ : データ点
- データ点をまとめて、行列$X$と表します。ただし、$X_{n,i}$は$x_n$の$i$番目の要素とします。
1.2 モデル
$K \in \mathbb{N}$に対し、$K$個の混合要素を持つ混合ガウスモデルは次のように表されます:
$$
\begin{align}
p\left(X,Z \middle| \mu, \Sigma, \pi \right) =
\prod_{n=0}^{N-1} \prod_{k=0}^{K-1}
\left[\pi_k \mathcal{N}\left(x_n \middle| \mu_k , \Sigma_k \right) \right]^{z_{n,k}}
\end{align}
$$
ただし、
- $Z = (z_{n,k})_{ n \in \{ 0, 1, \dots, N-1 \} k \in \{ 0,1, \dots, K-1 \} }, \ z_{n,k}\in \{0,1 \}, \sum_{k=0}^{K-1}z_{n,k}=1$であり、$Z$は潜在変数を表す。
- $\pi_k \geq 0, \sum_{k=0}^{K-1} \pi_k = 1$
- $\mathcal{N}\left(x \middle| \mu_k , \Sigma_k \right)$は、期待値$\mu_k \in \mathbb{R}^d$、共分散行列$\Sigma_k$のガウス分布の確率密度。
$Z$についての和をとると、$X$についての分布が次のように得られます。
$$
\begin{align}
p \left( X \middle| \mu, \Sigma, \pi \right) =
\prod_{n=0}^{N-1} \left[ \sum_{k=0}^{K-1} \pi_k \mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right) \right]
\end{align}
$$
このモデルに対して、対数尤度$\log p(X|\mu, \Sigma, \pi)$を「最大化」するパラメーター$(\mu, \Sigma, \pi)$を求めます。
ただし、PRML9.2.1でも散々述べられているように、この最大化問題はill-posedです。したがって、求めるのはあくまで、対数尤度を極大にするパラメーターです3。
1.3 アルゴリズム
導出はPRMLにあるので全て省略しますが、EM algorithmでは次の計算を行い、対数尤度の「極大化」を行います。
- パラメーターの初期化を行う。
- E step : responsibility $\gamma_{n,k}$を次のように計算する:
$$
\begin{align}
\gamma_{n,k} := \frac{ \pi_{k} \mathcal{N}\left(x_n \middle| \mu_k , \Sigma_k \right) }{ \sum_{l=0}^{K-1} \pi_{l} \mathcal{N}\left(x_n \middle| \mu_l , \Sigma_l \right) }
\end{align}
$$ - M step : 上で計算したresponsibilityを用い、パラメーター$\mu, \Sigma, \pi$を次のように更新する
$$
\begin{align}
&{} N_k := \sum_{n=0}^{N-1} \gamma_{n,k} \
&{} \pi_k = \frac{N_k}{N}\
&{} \mu_k = \frac{1}{N_k} \sum_{n=0}^{N-1} \gamma_{n,k} x_n \
&{} \Sigma_k = \frac{1}{N_k} \sum_{n=0}^{N-1} \gamma_{n,k} (x_n-\mu_k)(x_n-\mu_k)^T
\end{align}
$$ - 収束基準が満たされるまで、2と3を繰り返す。
収束基準としては、
- パラメーターの変化が十分小さくなったら停止
- 対数尤度の変化が十分小さくなったら停止
等の基準が考えられますが、ここでは後者の対数尤度の値を基にした基準を用いてみることにします。
なお、対数尤度の表式は次で表されます:
$$
\begin{align}
\log p \left( X \middle| \mu, \Sigma, \pi \right) =
\sum_{n=0}^{N-1} \log \left[ \sum_{k=0}^{K-1} \pi_k \mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right) \right]
\end{align}
$$
2 数式からコードへ
ここまでの理論をもとにして、実装を行います。
混合ガウスモデルを表す、GaussianMixtureModel
クラスを作成します。
2.1 概要
GaussianMixtureModel
クラスに、以下のデータ属性とメソッドを持たせることにします:
- データ属性
-
K
: 混合要素の個数を表す整数($K$のこと) -
Pi
: ($K$, ) array,Pi[k]
= $\pi_k$ -
Mu
: ($K$, $d$) array,Mu[k,i]
= $\mu_{k,i}$ -
Sigma
: ($K$, $d$, $d$) array,Sigma[k,i,j]
= $\left( \Sigma_{k} \right)_{i,j}$
-
- メソッド
-
_init_params
: パラメーターPi
,Mu
, andSigma
の初期化を行うメソッド。入力データ点に依存させます。 -
_calc_nmat
: 次で定義される($N$, $K$) arrayNmat
を計算して返すメソッド:
Nmat[n,k]
= $\mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right)$。E-stepや対数尤度の計算などで用います。 -
_Estep
: 上記のE-stepの計算を行うメソッド。具体的にはresponsibility $\gamma$を計算して返します。 -
_Mstep
: 上記のM-stepの計算を行うメソッド。具体的には、$\gamma$を与えられたときに、Pi
,Mu
, andSigma
を計算して更新します。 -
calc_prob_density
: 与えられたデータ点$X$に対して、確率密度$p \left( X \middle| \mu, \Sigma, \pi \right) = \prod_{n=0}^{N-1} \left[ \sum_{k=0}^{K-1} \pi_k \mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right) \right]$を計算して返すメソッド。確率密度をプロットするために使います。 -
calc_log_likelihood
: 与えられたデータ点$X$に対して、対数尤度$\log p \left( X \middle| \mu, \Sigma, \pi \right) = \sum_{n=0}^{N-1} \log \left[ \sum_{k=0}^{K-1} \pi_k \mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right) \right]$を計算して返すメソッド。収束判定のために利用します。 -
fit
: fittingを行うメソッド。具体的には、_Estep
と_Mstep
を、収束基準が満たされるまで繰り返します。 -
predict_proba
: responsibility $\gamma$を計算して返すメソッド(_Estep
と同じです)。 -
predict
: 与えられたデータ点たちに対し、各点がどのクラスタに属すかを計算して返すメソッド。
-
以下に、各メソッドの詳細と対応するコードを書いていきます。
2.2 初期化: _init_params
_init_params
では、次のように初期化を行うことにします:
-
Pi
: 全ての要素が等しい値とします。 -
Mu
:K
個の点を、入力データ点から一様ランダムに選びます。 -
Sigma
:Sigma[k]
は、全て同一の共分散行列$\Sigma_0$とします。ただし、$\Sigma_0$は対角行列で、その$(i, i)$成分は、入力データの第$i$成分の分散にとります。
対応するコードは、以下の通りです。なお、結果の再現性を持たせるために、Mu
の初期化の際にはseedを指定できるようにしてあります。
def _init_params(self, X, random_state=None):
'''
Method for initializing model parameterse based on the size and variance of the input data array.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
'''
n_samples, n_features = np.shape(X)
rnd = np.random.RandomState(seed=random_state)
self.Pi = np.ones(self.K)/self.K
self.Mu = X[rnd.choice(n_samples, size=self.K, replace=False)]
self.Sigma = np.tile(np.diag(np.var(X, axis=0)), (self.K, 1, 1))
2.3 Nmat
の計算: _calc_nmat
_calc_nmat
では、Nmat[n,k]
= $\mathcal{N}\left(x_n \middle| \mu_k, \Sigma_k \right)$で定義されるarrayを計算して返します。
よりコードに直しやすいように書くと、次の計算を行います:
$$
\begin{align}
\verb| Nmat[n, k] | = \frac{1}{\sqrt{\det \verb|Sigma[k]|}} \exp\left( -\frac{1}{2} \sum_{i,j=1}^{d} \verb| Diff[n,k,i] L[k,i,j] Diff[n,k,j] | \right)
\end{align}
$$
ただし、
-
L
: (K, d, d) array,L[k,i,j]
= $\left( \Sigma_{k}^{-1} \right)_{i,j} $ -
Diff
: (N, K, d) array,Diff[n,k,i]
= $x_{n,i} - \mu_{k,i}$
対応するコードは次のように書けます。
ポイントとしては、
-
Diff
の計算にあたっては、X
とMu
をreshapeすることにより、numpy broadcastingを利用しています。 - 指数の部分の和には
numpy.einsum
を利用します(詳細は公式ドキュメント参照) - 複数の行列の逆行列や行列式は、
numpy.linalg.inv
やnumpy.linalg.det
を利用して一括して計算しています(こちらも詳細は公式ドキュメント参照。)。
def _calc_nmat(self, X):
'''
Method for calculating array corresponding $\mathcal{N}(x_n | \mu_k)$
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
Nmat : 2D numpy array
2-D numpy array representing probability density for each sample and each component,
where Nmat[n, k] = $\mathcal{N}(x_n | \mu_k)$.
'''
n_samples, n_features = np.shape(X)
Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )
L = np.linalg.inv(self.Sigma)
exponent = np.einsum("nkj,nkj->nk", np.einsum("nki,kij->nkj", Diff, L), Diff)
Nmat = np.exp(-0.5*exponent)/np.sqrt(np.linalg.det(self.Sigma)) / (2*np.pi)**(n_features/2)
return Nmat
2.4 E-step : _Estep
(N, K) arrayGam
を、responsibilityを表すarrayとします : Gam[n,k]
= $\gamma_{n,k}$。_Estep
では、Nmat
を用いて、次のようにGam
を計算します。
$$
\begin{align}
\verb|Gam[n,k]| = \frac{\verb|Nmat[n,k] Pi[k] |}{ \sum_{l=1} \verb| Nmat[n,l] Pi[l] |}
\end{align}
$$
これはシンプルですね。対応するコードも、そのまま翻訳するだけです。
def _Estep(self, X):
'''
Method for calculating the array corresponding to responsibility.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
Gam : 2D numpy array
2-D numpy array representing responsibility of each component for each sample in X,
where Gamt[n, k] = $\gamma_{n, k}$.
'''
n_samples, n_features = np.shape(X)
Nmat = self._calc_nmat(X)
tmp = Nmat * self.Pi
Gam = tmp/np.reshape(np.sum(tmp, axis=1), (n_samples, 1) )
return Gam
2.5 M-step : _Mstep
M-stepではいよいよ、パラメーターの更新を行います。arrayの表記で1節の数式を書き直すと、
$$
\begin{align}
\verb|Pi[k]| &= \frac{1}{N} \sum_{n=0}^{N-1} \verb| Gam[n,k] | \
\verb| Mu[k,i]| &= \frac{1}{N_k} \sum_{n=0}^{N-1} \verb| Gam[n,k] X_[n,i] | \
\verb| Sigma[k,i,j]| &= \frac{1}{N_k} \sum_{n=0}^{N-1} \verb| Gam[n,k] Diff[n,k,i] Diff[n,k,j]|
\end{align}
$$
となります。
-
Sigma
の計算は、numpy.einsum
を用いて処理します。 - $N_k$による割り算では、broadcastingを活用します。
def _Mstep(self, X, Gam):
'''
Method for calculating the model parameters based on the responsibility gamma.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Gam : 2D numpy array
2-D numpy array representing responsibility of each component for each sample in X,
where Gamt[n, k] = $\gamma_{n, k}$.
'''
n_samples, n_features = np.shape(X)
Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )
Nk = np.sum(Gam, axis=0)
self.Pi = Nk/n_samples
self.Mu = Gam.T @ X / np.reshape(Nk, (self.K, 1))
self.Sigma = np.einsum("nki,nkj->kij", np.einsum("nk,nki->nki", Gam, Diff), Diff) / np.reshape(Nk, (self.K, 1, 1))
2.6 確率密度 calc_prob_density
EMアルゴリズムの実行には不要ですが、後で確率密度の可視化をしたいので、実装しておきます。
def calc_prob_density(self, X):
'''
Method for calculating the probablity density $\sum_k \pi_k \mathcal{N}(x_n | \mu_k)$
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
prob_density : 2D numpy array
'''
prob_density = self._calc_nmat(X) @ self.Pi
return prob_density
2.7 対数尤度の計算: calc_log_likelihood
対数尤度を計算して返します。収束判定に用いたり、異なる結果を比べたりするのに用います(実験の節で後述)。
これは、先ほどのcalc_prob_density
を用いると容易に書けます。
def calc_log_likelihood(self, X):
'''
Method for calculating the log-likelihood for the input X and current model parameters.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
loglikelihood : float
The log-likelihood of the input data X with respect to current parameter set.
'''
log_likelihood = np.sum(np.log(self.calc_prob_density(X)))
return log_likelihood
2.8 fit
ここまで書いてきたもろもろを組合せて、fittingを行うメソッドです。
-
_Estep
と_Mstep
を繰り返し実行します。 - 繰り返し回数が
max_iter
に達するか、対数尤度の変化がtol
以下になったら停止します。 - EMアルゴリズムの各ステップで対数尤度は非減少であることが理論的に示されているので、収束基準は「更新後の対数尤度 - 更新前の対数尤度 <
tol
」としてあります。
def fit(self, X, max_iter, tol, disp_message, random_state=None):
'''
Method for performing learning.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
max_iter : int
Maximum number of iteration
tol : float, positive
Precision. If the change of parameter is below this value, the iteration is stopped
disp_message : Boolean
Whether or not to show the message about the number of iteration
'''
self._init_params(X, random_state=random_state)
log_likelihood = - np.float("inf")
for i in range(max_iter):
Gam = self._Estep(X)
self._Mstep(X, Gam)
log_likelihood_old = log_likelihood
log_likelihood = self.calc_log_likelihood(X)
if log_likelihood - log_likelihood_old < tol:
break
if disp_message:
print(f"n_iter : {i}")
print(f"log_likelihood : {log_likelihood}")
2.9 predict_proba
与えられたデータに対してresponsibilityを求めるメソッドです。実態は_Estep
を呼び出すだけなのですが、良く用いられる名前のメソッドなので別個に作ってみました。
コードはここでは省略します。
2.10 predict
最後に、与えられたデータがどのクラスタに属するかを返すメソッドを書きます。
def predict(self, X):
'''
Method for make prediction about which cluster input points are assigned to.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
pred : 1D numpy array
1D numpy array, with dtype=int, representing which class input points are assigned to.
'''
pred = np.argmax(self.predict_proba(X), axis=1)
return pred
2.11 全体のコード
GaussianMixtureModel
クラスのコードは、以下のようになります。
class GaussianMixtureModel:
def __init__(self, K):
self.K = K
self.Pi = None
self.Mu = None
self.Sigma = None
def _init_params(self, X, random_state=None):
'''
Method for initializing model parameterse based on the size and variance of the input data array.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
'''
n_samples, n_features = np.shape(X)
rnd = np.random.RandomState(seed=random_state)
self.Pi = np.ones(self.K)/self.K
self.Mu = X[rnd.choice(n_samples, size=self.K, replace=False)]
self.Sigma = np.tile(np.diag(np.var(X, axis=0)), (self.K, 1, 1))
def _calc_nmat(self, X):
'''
Method for calculating array corresponding $\mathcal{N}(x_n | \mu_k)$
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
Nmat : 2D numpy array
2-D numpy array representing probability density for each sample and each component,
where Nmat[n, k] = $\mathcal{N}(x_n | \mu_k)$.
'''
n_samples, n_features = np.shape(X)
Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )
L = np.linalg.inv(self.Sigma)
exponent = np.einsum("nkj,nkj->nk", np.einsum("nki,kij->nkj", Diff, L), Diff)
Nmat = np.exp(-0.5*exponent)/np.sqrt(np.linalg.det(self.Sigma)) / (2*np.pi)**(n_features/2)
return Nmat
def _Estep(self, X):
'''
Method for calculating the array corresponding to responsibility.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
Gam : 2D numpy array
2-D numpy array representing responsibility of each component for each sample in X,
where Gamt[n, k] = $\gamma_{n, k}$.
'''
n_samples, n_features = np.shape(X)
Nmat = self._calc_nmat(X)
tmp = Nmat * self.Pi
Gam = tmp/np.reshape(np.sum(tmp, axis=1), (n_samples, 1) )
return Gam
def _Mstep(self, X, Gam):
'''
Method for calculating the model parameters based on the responsibility gamma.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Gam : 2D numpy array
2-D numpy array representing responsibility of each component for each sample in X,
where Gamt[n, k] = $\gamma_{n, k}$.
'''
n_samples, n_features = np.shape(X)
Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )
Nk = np.sum(Gam, axis=0)
self.Pi = Nk/n_samples
self.Mu = Gam.T @ X / np.reshape(Nk, (self.K, 1))
self.Sigma = np.einsum("nki,nkj->kij", np.einsum("nk,nki->nki", Gam, Diff), Diff) / np.reshape(Nk, (self.K, 1, 1))
def calc_prob_density(self, X):
'''
Method for calculating the probablity density $\sum_k \pi_k \mathcal{N}(x_n | \mu_k)$
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
prob_density : 2D numpy array
'''
prob_density = self._calc_nmat(X) @ self.Pi
return prob_density
def calc_log_likelihood(self, X):
'''
Method for calculating the log-likelihood for the input X and current model parameters.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
loglikelihood : float
The log-likelihood of the input data X with respect to current parameter set.
'''
log_likelihood = np.sum(np.log(self.calc_prob_density(X)))
return log_likelihood
def fit(self, X, max_iter, tol, disp_message, random_state=None):
'''
Method for performing learning.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
max_iter : int
Maximum number of iteration
tol : float, positive
Precision. If the change of parameter is below this value, the iteration is stopped
disp_message : Boolean
Whether or not to show the message about the number of iteration
'''
self._init_params(X, random_state=random_state)
log_likelihood = - np.float("inf")
for i in range(max_iter):
Gam = self._Estep(X)
self._Mstep(X, Gam)
log_likelihood_old = log_likelihood
log_likelihood = self.calc_log_likelihood(X)
if log_likelihood - log_likelihood_old < tol:
break
if disp_message:
print(f"n_iter : {i}")
print(f"log_likelihood : {log_likelihood}")
def predict_proba(self, X):
'''
Method for calculating the array corresponding to responsibility. Just a different name for _Estep
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
Gam : 2D numpy array
2-D numpy array representing responsibility of each component for each sample in X,
where Gamt[n, k] = $\gamma_{n, k}$.
'''
Gam = self._Estep(X)
return Gam
def predict(self, X):
'''
Method for make prediction about which cluster input points are assigned to.
Parameters
----------
X : 2D numpy array
2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.
Returns
----------
pred : 1D numpy array
1D numpy array, with dtype=int, representing which class input points are assigned to.
'''
pred = np.argmax(self.predict_proba(X), axis=1)
return pred
3 実験
それでは、ここまで書いてきたコードを実際に動かしてみます。
3.1 データ
使うデータは、make_blobs
で作った、次の2次元のtoy dataです。
3つのblobを生成しました。データ点の個数は100にしてあります。
3.2 K = 3
の結果
K=3
のモデルを用いて学習した結果は以下の通りです。tol = 1e-4
としてあります。
n_iter : 23
log_likelihood : -515.3221954638727
左側がクラスの予測値、右側が混合ガウスモデルが与える確率密度になっています。
自然にクラスタリングできているのが見て取れるかと思います。
ただ、初期化の仕方によっては、次のように、大きく異なった結果が得られることもあります:
n_iter : 29
log_likelihood : -533.6831531127243
対数尤度の値を比べてみると前者の方が大きいので、直感的にも数値的にも、前者の結果の方が尤もらしいのではないかと思われます4 。
3.3 様々なK
についての結果
ついでに、$K$の値を変えた場合の結果も見ておきましょう。
4 まとめ
この記事では、混合ガウスモデルに対する、EMアルゴリズムを実装しました。やや数式とコードが長かったですが、その対応は明確にできたのではないかと思います。
一方で、まだ私の理解の足りていないところも多々残っています:
- ill-posedな最尤推定をどう理解すべきか。
- 勾配法と比べたときの長所/短所が理解できていない。
- 収束性の理論的証明をきちんと理解していない。
これらは、今後の宿題としたいと思います。
次回は、同じモデルに対して、Bayes的な扱い(変分推論)を行う10章の内容を実装します!
-
特に、9.4節のKL divergenceを用いた導出 ↩
-
恐らく、次の文献に詳しいことが載っているかと思います: Wu, C. F. Jeff (Mar 1983). "On the Convergence Properties of the EM Algorithm". Annals of Statistics. 11 (1): 95-103. https://projecteuclid.org/euclid.aos/1176346060 ↩
-
「対数尤度を極大化するパラメーターを求める意義とは?」という疑問があるのですが、私はまだきちんと理解できていません。 ↩
-
ただ、対数尤度の最大化はill-posedな問題なので、値が大きいものが本当に「良い」のかどうかは、正直なところ判断がつかないです。 ↩