混合ガウス分布とEMアルゴリズムを解説して実装してみる

こんにちは、ぐぐりら(@guglilac)です。

前回の記事でk-means法を解説して実装までしてみたので、勢いに乗って今回はEMアルゴリズムについて解説してみようと思います。実装もします。


 (なぜクラスタリングばかりやっているのか。)

混合分布には月並みですが、多次元正規分布を使ったものを用いてみます。

混合ガウス分布とは

混合ガウス分布の前に、混合分布について説明します。

混合分布とは、複数の確率分布が混合係数という重みを用いて重み付け平均をとって作られる確率分布のことです。

混合される複数の分布の種類によってバリエーションがありますが、今回扱う混合ガウス分布はガウス分布(すなわち正規分布。通常多次元バージョンを扱います。)

数式で表すとこんな感じ。

\[P(X;\theta)=\sum_{i=1}^{K}\pi_iP(X;\mu_i,\Sigma_i)\]

$pi_i$が$i$番目の確率分布に対応する混合率になります。
ガウス分布を用いる場合は、パラメータ$\theta$は$(\pi,\mu,\Sigma)$の事を指します。

「この混合分布から生成したデータ列を用いて各パラメータを推定したい」というのが今回の気持ちです。

自然な発想としては、得られたデータから対数尤度関数を計算して、これを最大化するようなパラメータをみつけてあげればよいのですが、混合分布の対数尤度関数は偏微分が解析的に計算できません。

そこで、EMアルゴリズムの出番です。

EMアルゴリズムとは

対数尤度関数の偏微分がうまく計算できないので、少しずつパラメータを修正して対数尤度が最大になるパラメータを探索するという手法を用います。

いきなりですが、手順を書いちゃおうと思います。
よくわかんなくてもこれやっとけば大丈夫。なんでこれがうまくいくかは導出のところで書こうと思います。

1.適当な値を各パラメータ$(\pi,\mu,\Sigma)$に設定する。(初期化)
2.更新式に従って各パラメータを更新する。
$j$ステップ目のパラメータから$j+1$ステップ更新式は
\[\pi_i^{j+1}=\frac{1}{n}\sum_{t=1}^n\gamma_i^j(t)\]
\[\mu_i^{j+1}=\frac{\sum_{t=1}^n\gamma_i^j(t)x_t}{\sum_{t=1}^n\gamma_i^j(t)}\]
\[\Sigma_i^{j+1}=\frac{\sum_{t=1}^n\gamma_i^j(t)x_tx_t^T}{\sum_{t=1}^n\gamma_i^j(t)}-\mu_i^{j+1}(\mu_i^{j+1})^T\]
ただし、
\begin{eqnarray*}\gamma_i^j(t)&=&P(i|x_t;\mu_i^{j},\Sigma_i^{j})\\&=&\frac{\pi_i^jP(x_t;\mu_i^{j},\Sigma_i^{j})}{\sum_{l=1}^{j}\pi_k^jP(x_t;\mu_i^{j},\Sigma_i^{j})}\end{eqnarray*}
3.新しいパラメータとデータ列で対数尤度を計算し増加が止まるまで2を繰り返す。
ちょっとした注意ですが、三本目の更新式の右辺に出てくる$\mu_i^{j+1}$は新しく更新したばっかりの$\mu_i^{j+1}$であることに注意。なお、Tは転置を表しています。

$\gamma_i^j(t)$は$j$ステップ目においてデータ$x_t$が$i$番目の確率分布から現れたデータである確率を表しています。この割合を用いてより尤度が高くなるような各パラメータを推定して更新します。

EMアルゴリズムのEMはEステップとMステップを繰り返すことに由来します。

EステップはExpectation
MステップはMaximization

を表しています。

上のアルゴリズム中の処理でいうと、

Eステップは$\gamma_i^j(t)$を求める部分
Mステップは各パラメータの更新部分

に該当します。

更新式の導出

では、更新式の導出をしましょう。
細かい計算は省略します。きつかったら数学をやり直しましょう。 (たまに厳しい)

1ステップ前のパラメータ$\theta'$が与えられた時に尤度関数を最大化するパラメータ$\theta$は、次のように定義されるQ関数を最大化する$\theta$と同じです。

\[Q(\theta|\theta')=\sum_{z^n}P(z^n|x^n;\theta')\log P(x^n,z^n;\theta)\]

$z^n$は$n$個のデータがそれぞれどの分布から出たデータなのかを表す潜在変数です。
右辺はその全ての組み合わせについて和をとっている形になっています。

この式を変形すると次のようになります。
\begin{eqnarray*}Q(\theta|\theta')&=&\sum_{z^n}P(z^n|x^n;\theta')\log P(x^n,z^n;\theta)\\&=&\sum_{t=1}^n\sum_{i=1}^k\gamma_i'(t)\log \pi_iP_i(x_t;\mu_i,\Sigma_i)\end{eqnarray*}

1ステップ前のパラメータ$\theta'$を定数とみなして、(あと当然データ列&x_t&も定数として)Q関数を最大化するパラメータ$\theta$を求めれば良いです。

ただ今回はパラメータ$\theta$のうち、混合係数$\pi$はすべて足すと1になるという制約があるので、単純に偏微分して0に等しい!とやればいいわけではありません。

制約つき最適化問題ですね。ラグランジュの未定乗数法を用います。

\[f(\theta)=Q(\theta|\theta')+\lambda\sum_{i=1}^k\pi_i\]
として、
$\frac{\partial f}{\partial \theta}=0$と$\sum_{i=1}^k\pi_i=1$を連立してあげればいいでしょう。

更新式はそれぞれ、$\frac{\partial f}{\partial \pi_i}=0$,$\frac{\partial f}{\partial \mu_i}=0$,$\frac{\partial f}{\partial \Sigma_i}=0$から得られます。

2本目はスカラーをベクトルで微分、3本目はスカラーを行列で微分しています。
このような演算は機械学習界隈では頻出なので慣れておくといいと思います。

この辺の記事が参考になるかも。

「ベクトルで微分・行列で微分」公式まとめ

実装

実際に混合ガウス分布に従うデータを生成して、得られたデータ列をEMアルゴリズムによってクラスタリングしてみます。

あ、EMアルゴリズムでクラスタリングするっていうのは、EMアルゴリズムによって推定した複数の分布の中からそのデータが一番出やすい分布に分類するという感じです。

$\gamma_i^j(t)$は$j$ステップ目においてデータ$x_t$が$i$番目の確率分布から現れたデータである確率を表しているので、この値が最大になるインデックス$i$に分類してしまえば良さそうです。

Python3で実装しました。結果の後ろにソースコードも載せてあるので参考にしてみてください。

前回のk-means法の時と似たようなデータをクラスタリングします。


これを4つのグループにクラスタリングします。うまくいくといいな〜〜

結果

イテレーションを5回行うごとにその時点でのパラメータを用いてクラスタリングしたものを載せていきます。

黒丸は推定した各グループの中心を表しています。等高線は各分布のガウス分布を表しています。2次元の確率分布も等高線を使えばわかりやすいですね。
 5回終了時。割と初期化した平均点が正解に近いところだったため、最初から比較的うまくクラスタリングされています。

初期化をランダムにしている部分があるので、同じコードを試しても結果は異なります。今回は運が良いですね。
 10回終了時。少し落ち着いてきました。
 15回終了時。だいたいもう終わってます。

20回で収束しました!綺麗にクラスタリングされています。わーい

対数尤度関数の値はこんな感じ。上昇して行って最後には増加が止まるので収束していることがわかります。
iter:1 likelihood:-3842.720847712288
iter:2 likelihood:-3649.8295196239096
iter:3 likelihood:-3483.322755557905
iter:4 likelihood:-3345.8369806950527
iter:5 likelihood:-3271.919667255609
iter:6 likelihood:-3242.2223846206166
iter:7 likelihood:-3228.6747589443034
iter:8 likelihood:-3221.5026106295704
iter:9 likelihood:-3217.1428785467497
iter:10 likelihood:-3213.7649598658277
iter:11 likelihood:-3209.998035905365
iter:12 likelihood:-3203.6673419243098
iter:13 likelihood:-3191.6960964562663
iter:14 likelihood:-3176.4622164239267
iter:15 likelihood:-3169.8894090988674
iter:16 likelihood:-3169.68998823923
iter:17 likelihood:-3169.6892371688186
iter:18 likelihood:-3169.6891464002156
iter:19 likelihood:-3169.6891259683107
iter:20 likelihood:-3169.6891213337417
混合係数の推定値はこんな感じになりました。


ピンク
0.099975030.19976180.299940960.40032222

生成元の混合係数は次のように設定したので、ほぼ正しく推定できてることがわかります。


ピンク
0.10.20.30.4


ソースコード

今回使ったプログラムです。
自分でも試してみたい!という方は参考にしてみてください。
データ生成部分は載せないので、自分で作るか元々あるデータをクラスタリングしてみてください!
import numpy as np
from genarateData import genarateData, plot
import math
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt

Colors = ["red", "green", "m", "blue"]
N = 1000
iter_num = 100
Means = [[2, 2], [-2, -2], [2, -2], [-2, 2]]
Covs = [[[0.5, 0], [0, 0.25]], [[0.3, 0], [0, 0.6]],
        [[0.5, 0], [0, 0.25]], [[0.3, 0], [0, 0.6]]]
Pi = [0.1, 0.2, 0.3, 0.4]
K = len(Pi)
d = len(Means[0])


def gaussian(x, y, mu, sigma):
    x = np.matrix([x, y])
    mu = np.matrix(mu)
    sig = np.matrix(sigma)
    assert sig.shape == (d, d), sig.shape
    if np.linalg.det(sig) * (2 * np.pi)**sig.ndim < 0:
        print("error!")
    a = np.sqrt(np.linalg.det(sig) * (2 * np.pi)**sig.ndim)
    b = np.linalg.det(-0.5 * (x - mu) * sig.I * (x - mu).T)
    return np.exp(b) / a


def likelihood(X, Y, pi, mu, sigma):
    L = sum([np.log(sum([pi[i] * gaussian(X[t], Y[t], mu[i], sigma[i]) for i in range(K)]))
             for t in range(N)])
    return L


def initialize(X, Y, classes):
    pi = [1 / classes] * classes
    mu = [np.random.uniform(min(X), max(X), d) for _ in range(classes)]
    sigma = [np.eye(d) for _ in range(classes)]
    return pi, mu, sigma


def expectation(X, Y, pi, mu, sigma):
    r = []
    for t in range(len(X)):
        denom = sum([pi[l] * gaussian(X[t], Y[t], mu[l], sigma[l])
                     for l in range(K)])
        numer = np.array([pi[i] * gaussian(X[t], Y[t], mu[i], sigma[i])
                          for i in range(K)])
        r.append(numer / denom)
    assert np.array(r).shape == (N, K)
    return r


def maximization(X, Y, r):
    N_k = np.sum(r, axis=0)
    assert len(N_k) == K
    pi = N_k / N
    assert len(pi) == K
    mu = []
    sigma = []
    for i in range(K):
        mu_i = np.zeros(d)
        for t in range(N):
            mu_i += r[t][i] * np.array([X[t], Y[t]])
        mu_i /= N_k[i]
        mu.append(mu_i)
        sigma_i = np.zeros((d, d))
        for t in range(N):
            temp = np.array([X[t], Y[t]]) - mu[i]
            sigma_i += r[t][i] * \
                np.matrix(temp).reshape(d, 1) * np.matrix(temp).reshape(1, d)
        sigma_i /= N_k[i]
        sigma.append(sigma_i)
    assert np.array(mu).shape == (K, d), mu.shape
    assert np.array(sigma).shape == (K, d, d), sigma.shape
    return pi, mu, sigma


def show_param(pi, mu, sigma):
    print("pi:{}".format(pi))
    print("mu:")
    for i in range(K):
        print(mu[i])
    print("sigma:")
    for i in range(K):
        print(sigma[i])


def EM_plot(X, Y, r, iter, mu, sigma):
    plt.title = "iter{}".format(iter)
    # クラスタリング
    for t in range(N):
        i = np.argmax([r[t][i] for i in range(K)])
        plt.scatter(X[t], Y[t],
                    c=Colors[i], marker="o", s=5)
    # 平均点と等高線の表示
    xlist = np.linspace(min(X), max(X), 100)
    ylist = np.linspace(min(Y), max(Y), 100)
    x, y = np.meshgrid(xlist, ylist)
    for i in range(K):
        plt.scatter(mu[i][0], mu[i][1], c="k", marker="o")
        z = bivariate_normal(x, y, np.sqrt(sigma[i][0][0]), np.sqrt(
            sigma[i][1][1]), mu[i][0], mu[i][1], sigma[i][0][1])
        cs = plt.contour(x, y, z, 3, colors='k', linewidths=1)

    plt.show()


def EM(X, Y, classes, iter_num):
    pi, mu, sigma = initialize(X, Y, classes)
    for j in range(iter_num):
        prev = likelihood(X, Y, pi, mu, sigma)
        r = expectation(X, Y, pi, mu, sigma)
        pi, mu, sigma = maximization(X, Y, r)
        now = likelihood(X, Y, pi, mu, sigma)
        print("iter:{} likelihood:{}".format(j + 1, now))
        if now - prev < 0.00001:
            print("converged!")
            show_param(pi, mu, sigma)
            EM_plot(X, Y, r, j + 1, mu, sigma)
            return pi, mu, sigma
        if (j + 1) % 5 == 0:
            EM_plot(X, Y, r,  j + 1, mu, sigma)
    print("not converged")
    return pi, mu, sigma


if __name__ == "__main__":
    X, Y = genarateData(Means, Covs, Pi, N)
    plot(X, Y, Pi)
    EM(X, Y, K, iter_num)

まとめ

ふう。長かった。
解説して実装して、、ってやると力つくかんじするけど、やっぱ大変だ笑

前回のk-means法と同様にクラスタリングの手法の解説だったので、合わせて読むと知識が整理されるかもしれません。よかったら読んでみてね。


実装はそんなに難しいところはないと思うけれど、初期化した値によってはうまく計算が進まないこともありました。(オーバーフローや、共分散行列が半正定値じゃ無くなっちゃうとか)

結局、平均ベクトルのみランダムに初期化して、混合係数や共分散行列はそれっぽい値で初期化しました。
こうすれば計算が落ちることはなくなったので「これでいっか笑」となりました。

何か質問等ありましたらコメント下さい!

どうもありがとうございました。




コメント