[ベイズ推論による機械学習入門]多次元ガウス混合モデルにおける推論(ギブスサンプリング編)

須山敦志著の『ベイズ推論による機械学習入門』の第4章を読みました。4.4章 ガウス混合モデルにおける推論 をpythonで実装してみました。全コードは私のgitにあります。

記号の使い方

小文字の太字は縦ベクトル、大文字の太字は行列は表します。

多次元ガウス分布

$D$次元ガウス分布を、平均パラメータベクトル$\boldsymbol{\mu}\in \mathbb{R}^D$, 共分散行列${\bf \Sigma}\in \mathbb{R}^{D \times D}$を用いて
$$
{\mathcal N}({\bf x}|\boldsymbol{\mu},{\bf \Sigma})=\frac{1}{\sqrt{(2\pi)^D|{\bf \Sigma}|}} {\rm exp} \left\{ -\frac{1}{2}({\bf x}-\boldsymbol{\mu})^{\rm T}{\bf \Sigma^{-1}({\bf x}- \boldsymbol \mu)} \right\}
$$
で定義します。簡略のために精度行列${\bf \Lambda}={\bf \Sigma}^{-1}$を用います。

ウィシャード分布

$D\times D$の正定値行列${\bf \Lambda}$を生成するウィシャード分布を$D\times D$の正定値行列${\bf W}$を用いて次で表す。

$$
{\mathcal W}({\bf \Lambda}|\nu,{\bf W})= N |{\bf \Lambda}|^{\frac{\nu-D-1}{2}}{\rm exp}\left( -\frac{1}{2}{\rm Tr}({\bf W}^{-1}{\bf \Lambda}) \right)
$$

$N$は規格化定数で、パラメータ$\nu$は$\nu>D-1$を満たす。

カテゴリ分布

1 of K ベクトル${\bf s}$を生成するカテゴリ分布を次で表す。
$$
\rm{Cat}({\bf s}|\boldsymbol \pi) = \prod_{k=1}^K \pi_k^{s_k}
$$

ここで$\boldsymbol \pi= (\pi_1,\cdots , \pi_K)^{\rm T}$は分布を決めるパラメータで、$\pi_k \in (0,1)$かつ$\sum_{k=1}^K \pi_k = 1$を満たす。

ディレクレ分布

ベータ分布を多次元に拡張したディレクレ分布を次で表す。

$$
\rm{Dir}({\boldsymbol \pi}|\boldsymbol \alpha ) = N \prod_{k=1}^K \pi_k^{\alpha_k-1}
$$

ここで$N$は規格化定数で$K$次元のパラメータ$\boldsymbol \alpha$のそれぞれの要素は正の実数値です。ディレクレ分布は$\pi_k \in (0,1)$かつ$\sum_{k=1}^K \pi_k = 1$を満たす$\boldsymbol \pi$を生成する。

理論

詳細は『ベイズ推論による機械学習入門』の4.4章を読んでください。ここでは、前提と結論のみを書くに留めます。

前提

$K$個のガウス分布

\begin{align}
p({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k) = {\mathcal N}({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k^{-1})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,k=1,2,\cdots,K
\end{align}

から$N$個のデータ${\bf X}=\{{\bf x_1},\cdots,{\bf x}_N\}$を得る。${\bf x}_n(n=1,2,\cdots,N)$は$D$次元のベクトル。$\boldsymbol{\mu}_k$,${\bf \Lambda}_k$は$k$番目のガウス分布の平均パラメータベクトルと精度行列であり、$\boldsymbol{\mu}_k\in \mathbb{R}^D$, ${\bf \Lambda}_k\in \mathbb{R}^{D \times D}$である。

ここで、各${\bf x}_n$がどのガウス分布によって生成されたのかを表現する$K$次元のベクトル${\bf s}_n$を導入する。${\bf x}_n$が$k$番目のガウス分布由来であるとき、${\bf s}_n$の第$k$成分$s_{n,k}$が$1$で、他の成分は$0$になる。ここで、$\boldsymbol{\mu}=\{\boldsymbol{\mu}_1,\cdots,\boldsymbol{\mu}_K\}$,$\boldsymbol{\Lambda}=\{\boldsymbol{\Lambda}_1,\cdots,\boldsymbol{\Lambda}_K\}$とすれば、さらに

\begin{align}
p({\bf x_n}|{\bf s}_n,{\bf \mu}_k,{\bf \Lambda}_k) = \prod_{k=1}^K {\mathcal N}({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k^{-1})^{s_{n,k}}
\end{align}

以下のアルゴリズムの目的は、${\bf X}$から$\boldsymbol{\mu}$,${\bf \Lambda}$,${\bf S}$を求めることである。ただし${\bf S}=\{{\bf s}_1,\cdots,{\bf s}_N\}$である。

この観測モデルに対応する事前分布として、ガウス分布の共役事前分布であるガウス・ウィシャード分布を採用する。すなわち、ハイパーパラメータ${\bf m} \in \mathbb{R}^D,\beta \in \mathbb{R}^+,{\bf W} \in \mathbb{R}^{D\times D},\nu > D-1$を用いて

\begin{align}
p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) = {\mathcal N}(\boldsymbol{\mu}_k|\boldsymbol{m},(\beta\boldsymbol{\Lambda}_k)^{-1}){\mathcal W}({\bf \Lambda}_k|\mu,{\bf W})
\end{align}

と書ける。なお、それぞれのクラスタの混合比率$\boldsymbol{\pi}=(\pi_1,\cdots,\pi_K)^{\rm T}$より${\bf s}_n$が決定する。

\begin{align}
p({\bf s}_n|\boldsymbol \pi) = {\rm Cat}({\bf s}_n|\boldsymbol \pi)
\end{align}

$p(\boldsymbol \pi)$はカテゴリ分布の共役事前分布であるディレクレ分布を用いて

\begin{align}
p(\boldsymbol \pi) = {\rm Dir}(\boldsymbol \pi|\boldsymbol \alpha)
\end{align}

と表現される。$\boldsymbol \alpha$はディレクレ分布のハイパーパラメータである。以上の議論をまとめて、モデル全体の設計図は以下の通りである。

\begin{align}
p({\bf X},{\bf S},\boldsymbol \mu,{\bf \Lambda},\boldsymbol \pi) &= p({\bf X}|{\bf S},\boldsymbol \mu,{\bf \Lambda})p({\bf S}|\boldsymbol \pi)p(\boldsymbol{\mu},\boldsymbol{\Lambda})p(\boldsymbol \pi) \\
&= \left\{ \prod_{n=1}^N p({\bf x_n}|{\bf s}_n,{\bf \mu}_k,{\bf \Lambda}_k)p({\bf s}_n|\boldsymbol \pi)\right\} \left\{\prod_{k=1}^K p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k)\right\} p(\boldsymbol \pi)
\end{align}

結論

以下のようなギブスサンプリングアルゴリズムが導出される。

\begin{align}
{\bf s}_n &\sim {\rm Cat}({\bf s}_n|\boldsymbol \eta_n) \\
{\bf \Lambda}_k &\sim {\mathcal W}({\bf \Lambda}_k | \hat{\nu}_k,\hat{{\bf W}}_k) \\
\boldsymbol \mu_k &\sim {\mathcal N}(\boldsymbol \mu_k | \hat{\bf m}_k,(\beta{\bf \Lambda}_k)^{-1}) \\
\boldsymbol \pi &\sim {\rm Dir}(\boldsymbol \pi | \hat{\boldsymbol \alpha}) 

\end{align}

ここで

\begin{align}
\eta_{n,k} &\propto {\rm exp}\left(-\frac{1}{2}  ({\bf x}_n-\boldsymbol{\mu}_k)^{\rm T} {\bf \Lambda} ({\bf x}_n-\boldsymbol{\mu}_k) +\frac{1}{2} {\rm log}|{\bf \Lambda}_k| + {\rm log} \pi_k \right) \\
\hat{\alpha_k} &= \sum_{n=1}^N s_{n,k} +\alpha_k \\
\hat{\beta}_k &= \sum_{n=1}^N s_{n,k} + \beta \\
\hat{\nu}_k &= \sum_{n=1}^N s_{n,k} + \nu \\
\hat{{\bf m}}_k &= \frac{\sum_{n=1}^Ns_{n,k}{\bf x}_n+\beta{\bf m}}{\hat{\beta}_k} \\
\hat{{\bf W}}^{-1}_k &= \sum_{n=1}^N s_{n,k} {\bf x}_n{\bf x}_n^{\rm T} + \beta {\bf m} {\bf m}^{\rm T}-\beta_k {\bf m}_k {\bf m}^{\rm T}_k + {\bf W}^{-1}
\end{align}

である。

実装

まずは2次元ガウス混合モデルのギブスサンプリングを実装する。

number_of_cluster$(=K)$種類の2次元ガウス分布からnumber_of_data$(=N)$個のデータを得る関数datasetを用意する。show_initialはTrueでnumber_of_data個のデータを2dplotする。正解となる平均ベクトルexpectation_vectors$(=\{\boldsymbol \mu_1,\cdots,\boldsymbol \mu_K\})$と共分散行列convariance_matrixs$(=\{{\bf \Lambda}^{-1}_1,\cdots,{\bf \Lambda}^{-1}_K\})$を設定する。デフォルトでnumber_of_cluster=3としているが、number_of_clusterを変化させる場合は、適宜、設定する。

def dataset(number_of_data,number_of_cluster,show_initial):

    expectation_vectors = [[] for i in range(number_of_cluster)]

    expectation_vectors[0] = [+12,+6]
    expectation_vectors[1] = [0,-6]
    expectation_vectors[2] = [-13,+7]

    convariance_matrixs = [[[]] for i in range(number_of_cluster)]

    convariance_matrixs[0] = [[3,1],[1,12]]
    convariance_matrixs[1] = [[11,2],[2,9]]
    convariance_matrixs[2] = [[5,0],[0,18]]

    values = [[] for i in range(number_of_cluster)]

    for k in range(number_of_cluster):
        values[k] = np.random.multivariate_normal(expectation_vectors[k],convariance_matrixs[k],int(number_of_data/number_of_cluster))

    data = np.vstack(values)

    if show_initial:
        plt.scatter(data[:,0],data[:,1])
        plt.xlim([-20,20])
        plt.ylim([-20,20])
        plt.xlabel("x",fontsize=16)
        plt.ylabel("y",fontsize=16)
        plt.show()

    return data

s_vector(=${\bf s}_n$)をカテゴリ分布からサンプリングする。eta[n][k]は$\eta_{n,k}$を意味する。

pi_vectorは${\boldsymbol\pi}$、predict_precision_matrixsとpredict_expectation_vectorsは予測された正解${\bf \Lambda},{\boldsymbol \mu}$を意味する。

def sampling_svectors(number_of_data, number_of_cluster, pi_vector, predict_precision_matrixs, predict_expectation_vectors):

    eta = np.zeros((number_of_data,number_of_cluster))

    for n in range(0,number_of_data):
        for k in range(0,number_of_cluster):
            eta[n][k] = math.exp(-0.5*(x[n] - predict_expectation_vectors[k] ) @ predict_precision_matrixs[k] @ (x[n] - predict_expectation_vectors[k] ) + 0.5*math.log(np.linalg.det(predict_precision_matrixs[k])) + math.log(pi_vector[k]))

    for n in range(0,number_of_data):
        normalized_const = 0
        for k in range(0,number_of_cluster):
            normalized_const = normalized_const + eta[n][k]

        s_vectors[n] = np.random.multinomial(1,eta[n]/normalized_const)

    return s_vectors

次にpi_vector($={\boldsymbol \pi}$)をディレクレ分布からサンプリングする。ここでsum2は$\sum_{n=1}^Ns_{n,k}$を意味する。

def sampling_pi_vector(number_of_data,number_of_cluster,s_vectors):

    alpha_vector = np.ones(number_of_cluster)   
    new_alpha_vector = np.zeros(number_of_cluster)


    for k in range(0,number_of_cluster):

        sum2 = 0

        for n in range(0,number_of_data):
            sum2 = sum2 + s_vectors[n][k]

        new_alpha_vector[k] = sum2 + alpha_vector[k]

    pi_vector = np.random.dirichlet(new_alpha_vector)

    return pi_vector

次に、predict_precision_matrixsとpredict_expectation_vectorsをウィシャード分布とガウス分布からサンプリングする。ただし、実際にウィシャード分布からサンプリングするのは大変そうだったので、ウィシャード分布の期待値$<{\bf \Lambda}>=\nu {\bf W}$を採用した。

def sampling_predict_precision_matrixs_and_predict_expectation_vectors(number_of_data,number_of_cluster,s_vectors,):

    #ガウスウィシャード分布の超パラメータ
    beta = 1.0
    v = 2 # > D-1
    m_vector = np.ones(2)
    W_matrix = np.identity(2)

    new_beta = np.zeros(number_of_cluster)
    new_m_vector = np.zeros((number_of_cluster,2))


    for k in range(0,number_of_cluster):
        sum1 = 0
        sum2 = 0

        for n in range(0,number_of_data):
            sum1 = sum1 + s_vectors[n][k]*x[n]
            sum2 = sum2 + s_vectors[n][k]

        new_beta[k] = sum2 + beta 
        new_m_vector[k] = (sum1 + beta*m_vector)/new_beta[k]

    for k in range(0,number_of_cluster):
        sum1 = 0
        sum2 = 0

        for n in range(0,number_of_data):
            sum1 = sum1 + s_vectors[n][k]*(x[n].reshape(2,1) @ x[n].reshape(1,2))
            sum2 = sum2 + s_vectors[n][k]

        new_W_inverse = sum1 + beta*(m_vector.reshape(2,1) @ m_vector.reshape(1,2)) - new_beta[k]*(new_m_vector[k].reshape(2,1) @ new_m_vector[k].reshape(1,2)) + np.linalg.inv(W_matrix)
        new_v = sum2 + v

        predict_precision_matrixs[k] = new_v * np.linalg.inv(new_W_inverse)
        predict_expectation_vectors[k] = np.random.multivariate_normal(new_m_vector[k],np.linalg.inv(new_beta[k]*predict_precision_matrixs[k]))

    return predict_precision_matrixs,predict_expectation_vectors

二次元ガウス分布を返す変数を用意しておく。

def multivariate_Gaussian_distribution(x,y,precision_matrix,expectation_vector):

    cordinate = np.array([x,y])
    convariance_matrix = np.linalg.inv(precision_matrix)

    return 1.0/((2.0*math.pi)**2*np.linalg.det(convariance_matrix))*math.exp(-0.5*(cordinate-expectation_vector)@np.linalg.inv(convariance_matrix)@(cordinate-expectation_vector))

求めた$\bf s$を用いてクラスタリングし、散布図をカラーリングして、等高線をひく関数を用意する。

def show_answer(number_of_data,number_of_cluster,x,s_vectors,predict_precision_matrixs,predict_expectation_vectors,show_contour):

    groups = [[] for i in range(number_of_cluster)] #空の配列をnumber_of_cluster個用意する。

    for n in range(0,number_of_data):
        for k in range(0,number_of_cluster):
            if s_vectors[n][k] == 1:
                groups[k].append(x[n].tolist())

    for k in range(0,number_of_cluster):
        groups[k] = np.array(groups[k])
        plt.scatter(groups[k][:,0],groups[k][:,1])

    if show_contour:

        N = 60
        cordinate_x = np.linspace(-20,20,N)
        cordinate_y = np.linspace(-20,20,N)
        cordinate_X,cordinate_Y = meshgrid(cordinate_x,cordinate_y)

        value_of_z = np.zeros((N,N))

        for k in range(0,number_of_cluster):
            for i in range(0,N):
                for j in range(0,N):
                    value_of_z[j][i] = multivariate_Gaussian_distribution(cordinate_x[i],cordinate_y[j],predict_precision_matrixs[k],predict_expectation_vectors[k])

            plt.contour(cordinate_X,cordinate_Y,value_of_z) 

    plt.ylim([-20,20])
    plt.ylim([-20,20])
    plt.xlabel("x",fontsize=16)
    plt.ylabel("y",fontsize=16)
    plt.show()

最後にmain関数を用意する。

if __name__  == "__main__":
    number_of_cluster = 3 #クラスタの数
    trial_times = 80 #サンプリングを行う回数
    number_of_data = 600 #入力データ数

    assert number_of_data % number_of_cluster == 0, "入力データ数をクラス数の倍数にしてください。"

    show_initial = True #分類前の２次元データ分布を見るか否か。
    show_contour = True #分類後に正解等高線を表示するか否か。

    #初期化
    predict_expectation_vectors = np.zeros((number_of_cluster,2))
    predict_precision_matrixs = np.zeros((number_of_cluster,2,2))
    for i in range(0,number_of_cluster):
        predict_precision_matrixs[i] = np.identity(2)

    pi_vector = np.ones(number_of_cluster)
    s_vectors = np.zeros((number_of_data,number_of_cluster))

    #データを渡す
    x = dataset(number_of_data,number_of_cluster,show_initial)

    #ギブスサンプリングを実行
    for i in range(0,trial_times):

        #sをサンプリングする。
        s_vectors = sampling_svectors(number_of_data, number_of_cluster, pi_vector, predict_precision_matrixs, predict_expectation_vectors)

        #predict_precision_matrixsとpredict_expectation_vectorsをサンプリングする。
        predict_precision_matrixs, predict_expectation_vectors = sampling_predict_precision_matrixs_and_predict_expectation_vectors(number_of_data, number_of_cluster, s_vectors)

        #pi_vectorをサンプリングする。
        pi_vector = sampling_pi_vector(number_of_data, number_of_cluster, s_vectors)

    #答えをプロット
    show_answer(number_of_data,number_of_cluster,x,s_vectors,predict_precision_matrixs,predict_expectation_vectors,show_contour)

    #最終結果を表示
    print("Trial回数:",trial_times,"\n")
    for k in range(0,number_of_cluster):
        print("グループ"+str(k)+"の平均\n",predict_expectation_vectors[k])
        print("グループ"+str(k)+"の共分散行列\n",np.linalg.inv(predict_precision_matrixs[k]),"\n")

結果

2次元ガウス混合モデル

1枚目の画像の600点を3つのクラスタに分類した後の様子が2枚目である。80回で正しくクラスタリングされている様子が分かる。

グループ0の平均
[-13.03527306 6.85378566]
グループ0の共分散行列
[[ 5.55089831 -0.06026977]
[-0.06026977 20.82277496]]

グループ1の平均
[-0.46999732 -6.09191948]
グループ1の共分散行列
[[9.62733973 1.28850163]
[1.28850163 9.60439715]]

グループ2の平均
[12.0153539 5.74785767]
グループ2の共分散行列
[[3.78214067 1.07161246]
[1.07161246 9.74919818]]

3次元ガウス混合モデル

3次元版でもやってみた！（ソースコードはこちら）