この記事では、PRMLの第9章で述べられている、混合ガウスモデルに対するEMアルゴリズムをpythonで実装します。
対応するjupyter notebookは筆者のgithubリポジトリにあります。連載全体の方針や、PRMLの他のアルゴリズムの実装については、連載のまとめページをご覧いただければと思います。

1 理論のおさらい

この節では、混合ガウスモデルに対するEMアルゴリズムの実装に必要な結果だけを、簡潔にまとめます。導出や、EMアルゴリズムの一般的な説明¹などは、本の方をご参照ください。
なお、EMアルゴリズムの収束性の証明などはPRMLには記載はなく、私もフォローできていません²。

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です。したがって、求めるのはあくまで、対数尤度を極大にするパラメーターです³。

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, and Sigmaの初期化を行うメソッド。入力データ点に依存させます。
- _calc_nmat : 次で定義される($N$, $K$) array Nmatを計算して返すメソッド: 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, and Sigmaを計算して更新します。
- 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

対数尤度の値を比べてみると前者の方が大きいので、直感的にも数値的にも、前者の結果の方が尤もらしいのではないかと思われます⁴ 。

3.3 様々な`K`についての結果

ついでに、$K$の値を変えた場合の結果も見ておきましょう。

$K=1$

$K=2$

$K=4$

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な問題なので、値が大きいものが本当に「良い」のかどうかは、正直なところ判断がつかないです。 ↩

PRML第9章の混合ガウスモデルに対するEMアルゴリズムをPythonで実装

1 理論のおさらい

1.1 設定

1.2 モデル

1.3 アルゴリズム

2 数式からコードへ

2.1 概要

2.2 初期化: _init_params

2.3 Nmatの計算: _calc_nmat

2.4 E-step : _Estep

2.5 M-step : _Mstep

2.6 確率密度 calc_prob_density

2.7 対数尤度の計算: calc_log_likelihood

2.8 fit

2.9 predict_proba

2.10 predict