pythonのscikit-learnで混合ガウスモデルを用いたクラスタリング

この記事では, pythonのscikit-learnで提供されている混合ガウスモデル(Gaussian Mixture Model, GMM)を用いたクラスタリングの実装について解説する.
目標としては, scikit-learnの公式ドキュメント[1]にある図をirisのデータセットに対して置き換える. そして図の下をPCA(Principal Component Analysis)を用いて, クラス分けの結果を可視化した図にすることである.

scikit-learnのオプションについて

始めに, 混合ガウス分布を利用するにあたってモデルの簡単な説明とscikit-learnのオプションの設定の意味を概説する.

混合ガウス分布モデルは, 次の式によって表現される.


$$ p(x) = \sum_{k=1}^{K} \pi_k N(x|\mu_k,\Sigma_k) $$

ここで, Kはクラスターの総数を示している. クラスの分け方としては, 各kに関してΣの中身が最大, つまりπ*Nが最大となるようなものにクラスを割り振ることとなる.
scikit learnで実装されている混合ガウス分布によるクラスタリングでは, π, μ, ΣをEM algorithmを用いて推定している. 理論については, [2][3]を参照のこと.

[4]みると, irisのデータセットに関してGMMを用いた結果が表示されている. x軸, y軸はそれぞれ, “sepal length (cm)”, “sepal width (cm)”である.

ここで, 4つのtype, “spherical” ,”diag”, “tied”, “full” の4つの設定がある. 公式ドキュメント[5]には次のように載っている.

“spherical”の場合は, 一つのクラスターに関して, 分散は一つしかない, つまりモデルは次のように表される.


$$ p(x) = \sum_{k=1}^{K} \pi_k N(x|\mu_k, \alpha_k^ I) $$

“diag”の場合は, 各クラスターに関して, 対角要素の分散が与えられる. つまり,


$$ p(x) = \sum_{k=1}^{K} \pi_k N(x|\mu_k,\Sigma_k) $$

において, データxの次元が3次元だとすると,


$$
\Sigma_k =
\begin{pmatrix}
\sigma_{k11} & 0 & 0 \\
0 & \sigma_{k22} & 0 \\
0 & 0 & \sigma_{k33}
\end{pmatrix}
$$

で与えられる.

“tied” は, クラスターに関わらず同じ分散共分散行列を持つことを述べている. 数式としては以下のように与えられる. 分散共分散行列から添字のkがなくなっていることに注意. 実際は, ガウス分布の制約から分散共分散行列は対称行列となる.


$$ p(x) = \sum_{k=1}^{K} \pi_k N(x|\mu_k,\Sigma) $$

“full”は制限なしのモデルである.

準備

実装に移っていく.
以下のモジュールをインポートしておく. 実行はjupyter notebook上を想定.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(context="paper" , style ="whitegrid",rc={"figure.facecolor":"white"})

from scipy import linalg
import itertools
# for Gaussian Mixture model for classification
from sklearn import mixture
# for PCA (principal component analysis)
from sklearn.decomposition import PCA

# for sample data
from sklearn import datasets

また, irisのデータセットを構造化されたデータと同じように加工しておく.

dataset = datasets.load_iris()
dic = {}
for k,v in zip(dataset.feature_names,dataset.data.T):
    dic[k] = v
df = pd.DataFrame(dic)
df.head().T

GMMの実装と可視化

実装内容は, ほぼ[1]参考. データの入力の仕方と変数の名前の付け方, fittingの際の初期値を変えての実行回数(defaultだと1回)を10回に設定しているところだけが相違点. 何回も初期値を変えて実行しているのは, EM algorithmにおいては最適解にいかない場合があるから. これは, GMM自体が設定不良問題になってしまっているから. 詳細は[2]のCh.9参照.

X = df.values
lowestBIC = np.infty
bic = []
nComponentsRange = range(1, 7)
cvTypes = ['spherical', 'tied', 'diag', 'full']
nInit = 10
for cvType in cvTypes:
    for nComponents in nComponentsRange:
        # Fit a Gaussian mixture with EM
        gmm = mixture.GaussianMixture(n_components=nComponents,
                                      covariance_type=cvType,n_init =nInit)
        gmm.fit(X)
        bic.append(gmm.bic(X))
        if bic[-1] < lowestBIC:
            lowestBIC = bic[-1]
            bestGmm = gmm
            
bic = np.array(bic)
colorIter = itertools.cycle(['navy', 'turquoise', 'cornflowerblue',
                              'darkorange'])

bars = []

# Plot the BIC scores
plt.figure(figsize=(8, 6),dpi=100)
ax = plt.subplot(111)
for i, (cvType, color) in enumerate(zip(cvTypes, colorIter)):
    xpos = np.array(nComponentsRange) + .2 * (i - 2)
    bars.append(plt.bar(xpos, bic[i * len(nComponentsRange):
                                  (i + 1) * len(nComponentsRange)],
                        width=.2, color=color))
plt.xticks(nComponentsRange)
plt.ylim([bic.min() * 1.01 - .01 * bic.max(), bic.max()])
plt.title('BIC score per model')
xpos = np.mod(bic.argmin(), len(nComponentsRange)) + .65 +
    .2 * np.floor(bic.argmin() / len(nComponentsRange))
plt.text(xpos, bic.min() * 0.97 + .03 * bic.max(), '*', fontsize=14)
ax.set_xlabel('Number of components')
ax.legend([b[0] for b in bars], cvTypes)

可視化に関しては, 以前の記事( pythonでk means clusteringの結果をPCAを用いて2次元と3次元で確認する ) に載っている関数を用いて行う. 関数部分のコードは省略した.

varNames = df.columns
# normalization
dfScaled = df.copy()[varNames]
dfScaled = (dfScaled - dfScaled.mean())/dfScaled.std()

# principal components analysis
pca = PCA(n_components=3)
trans= pca.fit_transform(dfScaled[varNames])

# insert results of GMM clustring 
dfScaled["clusterId"] = bestGmm.predict(X)

# these functions are defined in the previous article. 
summaryPCA(pca,varNames)
contriPCAPlot(pca,varNames)
clusterPCAPlot(trans,dfScaled)

結果の解釈について

iris datasetはもともと, "setosa", "versicolor", "virginica"の3種から構成されたデータセットであったが, BICに基づいた分類では次数が2の"full"モデルが最も良いと判断された. これはモデルの適合度(goodness of fit)を表す指標が, 現実を正確に表すわけではない例といえる.


以下の画像はscikit-learnが出している, 解析の際にどのalgorithmを使えば良いかに関してのcheat-sheet[6]である.

このcheet-sheetに従えば, クラスタリングがカテゴリーでクラスターの数が決まっていない場合は, "MeanShift"か"VBGMM"を用いることが推奨されている. しかし, どちらの方法もクラスターの数を定めなくても良い代わりに, hyper parameterの調整が必要になってくる. だが, そのhyper parameterを選ぶのは恣意的な作業で交差検定の作業が必要になってくる. 個人的には, それをやるくらいならば, 今回のように次数ごとにBICを計算して, 最もBICが良かったクラスターの数を採用するのが楽なように思える.

----------雑感(`・ω・´)----------
先のような議論は, reversible jump MCMCを使った時に感じたことである. reversible jump MCMCは, モデル間の推移も含めてパラメーターを推定しようという方法なのだが, モデルによっては十分に片方のパラメーターの探索が行われないという自体に陥ってしまう. それくらいならば, しっかりと一つ一つのモデルに対してパラメーターの探索をさせて, その上でモデルのgoodness of fitに関する指標で比較した方が良いと感じた.

[2]の文献のCh.9のEM algorithmの実装が今回のscikit-learnのGMMに該当して, Ch.10.2の変分法の実装が, scikit-learnのVBGMMの "weight_concentration_prior_type"を"dirichlet_distribution"にしたものに該当する. [7]参照のこと.

パターン認識と機械学習 ベイズ理論による統計的予測 上/C.M.ビショップ/元田浩/栗田多喜夫【3000円以上送料無料】

参考文献

[1]Gaussian Mixture Model Selection, https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_selection.html#sphx-glr-auto-examples-mixture-plot-gmm-selection-py
[2] Christopher M. Bishop, Pattern Recognition and Machine Learning
[3] pythonで混合正規分布実装, https://qiita.com/ta-ka/items/3e8b127620ac92a32864
[4] GMM covariances, https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html#sphx-glr-auto-examples-mixture-plot-gmm-covariances-py
[5] sklearn.mixture.GaussianMixture , https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture
[6] Choosing the right estimator, https://scikit-learn.org/stable/tutorial/machine_learning_map/index.html
[7] 2.1. Gaussian mixture models, https://scikit-learn.org/stable/modules/mixture.html

コメント

タイトルとURLをコピーしました