クラスタリングとGMM ― 複合システムモデルを同定する二つのアプローチ

やりたいこと

ある未知の代数的システム（与えられた代数方程式を満たす全ての変数の集合）$\mathcal{S}=\left\{\boldsymbol{x}|\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{0}\right\}$（$\boldsymbol{x}$はシステム変数）から標本群$\mathcal{X}=\{\boldsymbol{x}_{i}\}$（$i=1,\cdots,N$）が観測によって得られたとして，これらからシステムのモデル$\hat{\mathcal{S}}=\left\{\boldsymbol{x}|\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{0}\right\}$（$\boldsymbol{\theta}$はシステムパラメータ）を同定することを考えます．具体的には，$\boldsymbol{\theta}$を最小二乗法で求める問題になります．

\boldsymbol{\theta}=\mathop{\mathrm{arg~min}}_{\boldsymbol{\theta}}
\left\{
E(\boldsymbol{\theta};\mathcal{X})=\frac{1}{2}\sum_{i=1}^{N}\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{e}_{i}
\right\}

ただし，$\boldsymbol{e}_{i}\overset{\mathrm{def}}{=}\hat{\boldsymbol{f}}(\boldsymbol{x}_{i};\boldsymbol{\theta})$です．$E(\boldsymbol{\theta};\mathcal{X})$が「$\mathcal{X}$が与えられた下での$\boldsymbol{\theta}$の関数」であることに注意して下さい．$\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})$が$\boldsymbol{\theta}$に対し微分可能かつ$E(\boldsymbol{\theta};\mathcal{X})$が下に凸であれば，次の方程式を解くことにより$\boldsymbol{\theta}$が得られます．

\left(\frac{\partial E}{\partial\boldsymbol{\theta}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\sum_{i=1}^{N}\left(\frac{\partial\hat{\boldsymbol{f}}(\boldsymbol{x}_{i};\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right)^{\mathrm{T}}\hat{\boldsymbol{f}}(\boldsymbol{x}_{i};\boldsymbol{\theta})=\boldsymbol{0}

以下，$\boldsymbol{x}$，$\boldsymbol{f}(\boldsymbol{x})$，$\boldsymbol{\theta}$の次元がそれぞれ$n$，$m$，$p$であるとします．

【例1】 $\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{x}-\boldsymbol{\theta}$（$p=m=n$）の場合，

\left(\frac{\partial E}{\partial\boldsymbol{\theta}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\sum_{i=1}^{N}(\boldsymbol{x}_{i}-\boldsymbol{\theta})=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{\theta}=\frac{1}{N}\sum_{i=1}^{N}\boldsymbol{x}_{i}

つまり，$\boldsymbol{\theta}$は$\boldsymbol{x}_{i}$の平均です．

【例2】 $\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{x}-1$（$p=n, m=1$）の場合，

\left(\frac{\partial E}{\partial\boldsymbol{\theta}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\sum_{i=1}^{N}\boldsymbol{x}_{i}(\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{\theta}-1)=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{\theta}=\left(\sum_{i=1}^{N}\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{\mathrm{T}}\right)^{-1}\sum_{i=1}^{N}\boldsymbol{x}_{i}

ちなみにこの例の$\mathcal{S}$は，$\boldsymbol{\theta}$を法線方向ベクトルとする平面上の点群により構成されるシステムを表しています．右辺の計算は$\sum_{i=1}^{N}\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{\mathrm{T}}$が正則な場合にのみ可能で，このため少なくとも$N\geq p$でなければなりません．

次に，システム全体が$K$個の相異なる部分システム$\mathcal{S}_{k}$（$k=1,\cdots,K$）の直和$\mathcal{S}=\mathcal{S}_{1}+\cdots+\mathcal{S}_{K}$である場合を考えましょう．$\mathcal{S}_{k}$（$k=1,\cdots,K$）を，パラメータ$\boldsymbol{\theta}_{k}$を用いて$\hat{\mathcal{S}}_{k}=\left\{\boldsymbol{x}\left|\hat{\boldsymbol{f}}_{k}(\boldsymbol{x};\boldsymbol{\theta}_{k})=\boldsymbol{0}\right.\right\}$とモデル化します（$\boldsymbol{\theta}_{k}$の要素数は全て$p$であるとします）．$\mathcal{X}$から全ての$\boldsymbol{\theta}_{k}$を同定するためには，先にそれぞれの標本がどの部分システムに所属するか推定しなければなりません．しかしそれぞれの標本がどの部分システムに所属するか正しく推定するためには，全ての$\boldsymbol{\theta}_{k}$が同定されていなければなりません．この典型的なchicken-and-egg問題をどのように解くかが本記事の本題です．

クラスタリング（K-means法）によるモデル同定

K-means法

$K$-means法 は，全標本がそれぞれどの部分システムに所属するか，陽に推定する方法です．

• J. B. MacQueen,"Some Methods for Classification and Analysis of Multivariate Observations,"
  in Proceedings of the fifth Berkeley Symposium on Mathematical Statistics and Probability 1, University of California Press, 1967.
• R. O. Duda and P. E. Hart, "Pattern Classification and Scene Analysis," John Wiley & Sons, 1973.

アルゴリズムは次の通りです．

1. 観測された全標本を$K$個のグループに適当に分割し，各部分システム$\mathcal{S}_{k}$（$k=1,\cdots,K$）に所属する標本群（クラスタ）の初期候補$\mathcal{C}_{k}$（$k=1,\cdots,K$）とする．
2. 各$\mathcal{C}_{k}$に所属する標本群を用いて$\boldsymbol{\theta}_{k}$をそれぞれ同定する．
3. 各$\boldsymbol{x}_{i}$に対し$\boldsymbol{e}_{ik}=\hat{\boldsymbol{f}}_{k}(\boldsymbol{x}_{i};\boldsymbol{\theta}_{k})$（$k=1,\cdots,K$）を求め，最も小さい$\|\boldsymbol{e}_{ik}\|$を与える$\mathcal{C}_{k}$にその標本を改めて所属させる．
4. 3において所属変更された標本が無ければ終了．さもなければ2に戻る．

ただし，オリジナルの$K$-means法における平均および平均から各標本までの距離計算を，それぞれ$\boldsymbol{\theta}_{k}$の同定と$\boldsymbol{e}_{ik}$の計算に読み替えています．可視化すればアルゴリズムのプロセスは一目瞭然なのですが，横着してそれは他記事に委ねます．

良い結果を得るための初期クラスタ候補配置

$K$-means法の結果は初期クラスタ候補の配置に依存します．良い結果（クラスタ間の分離度が高い結果）を得るためには，初期クラスタ候補をお互いになるべく離れた所に配置すべきです．このような発想に基づく方法として，KKZ法 と呼ばれる方法が知られています．

• I. Katsavounidis, C. C. J. Kuo and Z. Zhang, "A New Initialization Technique for Generalized Lloyd Iteration," IEEE Signal Processing Letters, Vol. 1, No. 10, pp. 144--146, 1994.

具体的には次のように初期クラスタを作成します．

1. $\mathcal{X}$の元のうち，互いに最も離れた二つを$\boldsymbol{x}^{*}_{1}$，$\boldsymbol{x}^{*}_{2}$とする．
2. $k=3$とする．
3. $\mathcal{X}$から$\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{k-1}$を除いた$\mathcal{X}_{k}\overset{\mathrm{def}}{=}$ $\mathcal{X}\setminus\left\{\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{k-1}\right\}$の元の中から，$\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{k-1}$のうち最小距離を与える点への距離が最大となるものを$\boldsymbol{x}^{*}_{k}$とする．
4. $k<K$ならば，$k\leftarrow k+1$として3へ戻る．
5. 残りの全ての標本$\forall\boldsymbol{x}_{i}\in\mathcal{X}_{K+1}=\mathcal{X}\setminus\left\{\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{K}\right\}$を，$\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{K}$のうち最も近いものと同じグループに所属させ，初期クラスタ候補$\mathcal{C}_{k}$（$k=1,\cdots,K$）を確定する．

なお，KKZというのは上記引用文献の著者3人の頭文字をつなげたものです．

この方法は，観測された標本の中に外れ点（アウトライア）が混ざっている場合に，それに引っ張られて結果が悪化するという弱点があります．そこで，上記手順のうち3を

1. $\mathcal{X}$から$\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{k-1}$を除いた$\mathcal{X}_{k}\overset{\mathrm{def}}{=}\mathcal{X}\setminus\left\{\boldsymbol{x}^{*}_{1},\cdots,\boldsymbol{x}^{*}_{k-1}\right\}$の元それぞれについて$d_{i}=\min_{j=1,\cdots,k-1}d(\boldsymbol{x},\boldsymbol{x}^{*}_{j})$を求め，$\mathcal{X}_{k}$の元から$\boldsymbol{x}_{i}$が選ばれる確率を$\displaystyle\frac{d_{i}}{\sum_{j=1}^{N-k+1}d_{j}}$として，乱数的に一つ選択し$\boldsymbol{x}^{*}_{k}$とする．

と置き換えることで，確率的に外れ点の影響から逃れる方法も提案されています．これは $K$-means$++$法 と呼ばれます．

• D. Arthur and S. Vassilivitskii, "k-means++: the advantages of careful seeding," in Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithm, pp. 1027--1035, 2007.

いずれにしても確率的に外れ点を選択してしまう可能性は免れないので，$K$-means法を何度か試行して最も良い結果を採択する，というやり方が実用上は必須です．

結果の良し悪しの評価方法としては，次に説明する シルエット分析 がよく知られています．

• P. J. Rousseeuw, "Silhouettes: a Graphical Aid to the Interpretation and Validation of Cluster Analysis," Computational and Applied Mathematics, Vol. 20, pp. 53--65, 1987.

まず，全標本に対して次の$s_{i}$（シルエット）を計算します．

s_{i}=\frac{b_{i}-a_{i}}{\mathrm{max}\left\{a_{i},b_{i}\right\}}

ただし，

\displaylines{
a_{i}=D(\boldsymbol{x}_{i};\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}})
\\
b_{i}=\min_{\mathcal{C}\neq\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}}}D(\boldsymbol{x}_{i};\mathcal{C})
\\
D(\boldsymbol{x};\mathcal{C})\overset{\mathrm{def}}{=}
\sum_{\boldsymbol{x}_{j}\in\mathcal{C}}
\frac{d(\boldsymbol{x}_{j},\boldsymbol{x})}{\left|\mathcal{C}\right|}
}

で，$\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}}$は$\boldsymbol{x}_{i}$の所属するクラスタ，$\left|\mathcal{C}\right|$はクラスタ$\mathcal{C}$の要素数を表します．原典では$a_{i}$の計算において$\boldsymbol{x}_{i}$自身を$\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}}$から除いているのですが，その場合は$\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}}=\left\{\boldsymbol{x}_{i}\right\}$となったときに例外処理が必要となるので，本記事では敢えてそのようなことをしない定義を採用します．

$s_{i}$は$-1\sim 1$の値をとり得ます．定性的には，クラスタ$\mathcal{C}_{\left\{\boldsymbol{x}_{i}\right\}}$の凝集度が高いと$a_{i}$が小さく，クラスタ間の分離度が高いと$b_{i}$が大きくなるので，$\boldsymbol{x}_{i}$がうまくクラスタリングされていれば$s_{i}$は$1$に近く，そうでなければ$-1$に近くなります．平均シルエット$\bar{s}$

\bar{s}=\frac{1}{N}\sum_{i}^{N}s_{i}

が大きいほど，良好な結果が得られていると言えます．

クラスタ数はどのように選ぶ？

$K$が未知の場合は，幾つか異なる$K$について$\hat{\mathcal{S}}$を得た後，最も良い結果を与えた$K$を選択するのが良いでしょう．この評価にもシルエット分析が応用できます．

• R. C. de Amorim and C. Hennig, "Recovering the number of clusters in data sets with noise features using feature rescaling factors," Information Sciences, Vol. 324, pp. 126--145, 2015.

前節で説明した事柄から単純に考えれば，平均シルエット$\bar{s}$が最大となる$K$を選べば良いように思えますが，$K$を大きくして外れ点を全て独立したクラスタとすることで$\bar{s}$を大きくすることも，理論上可能です．したがって，各クラスタの持つ標本数に大きく隔たりが無いかどうかも考慮に入れなければなりません．実用上は，全てのクラスタの最大シルエット$\max_{\boldsymbol{x}_{i}\in\mathcal{C}_{k}}s_{i}$が平均シルエット$\bar{s}$より大きいものの中で，最小クラスタ数と最大クラスタ数の比$\displaystyle\frac{\max_{k=1,\cdots,K}\left|\mathcal{C}_{k}\right|}{\min_{k=1,\cdots,K}\left|\mathcal{C}_{k}\right|}$が最小となるような$K$を採用する，という規範が採用されます．

X-means法：クラスタを再帰的に自動分割するアルゴリズム

評価値に基づいて再帰的にクラスタを自動分割していく $X$-means法 という方法も提案されています．

• D. Pelleg and A. Moore, "$X$-means: Extending $K$-means with Efficient Estimation of the Number of Clusters," ICML-2000.

アルゴリズムは次の通りです（オリジナルのものを少々変更しています）．

1. まず$K=1$と仮定し，全標本から$\hat{\mathcal{S}}$を同定する．
2. モデルの同定精度を評価する（同定精度に対し単調減少な）指標値$I(\hat{\mathcal{S}})$を計算する．
3. 次に，システムを2つのサブシステムの直和$\hat{\mathcal{S}}^{\prime}=\hat{\mathcal{S}}_{\mathrm{A}}+\hat{\mathcal{S}}_{\mathrm{B}}$としてモデル化し，同定する．
4. 同様に$I(\hat{\mathcal{S}}^{\prime})$を計算する．
5. $I(\hat{\mathcal{S}}^{\prime})<I(\hat{\mathcal{S}})$ならば，$\hat{\mathcal{S}}_{\mathrm{A}}$，$\hat{\mathcal{S}}_{\mathrm{B}}$へのシステム分割によるモデル同定精度向上の効果があったと判断して$\hat{\mathcal{S}}$を$\hat{\mathcal{S}}^{\prime}$で置き換える．
6. さらに$\hat{\mathcal{S}}_{\mathrm{A}}$と$\hat{\mathcal{S}}_{\mathrm{B}}$のそれぞれについて，再帰的に上記2$\sim$5を繰り返す．

モデル同定精度を評価する指標$I(\hat{\mathcal{S}})$が重要です．オリジナルの$X$-means法で用いられているのは，ベイズ情報量基準 （Baysian Information Criterion, BIC ）です．

• G. Schwarz, "Estimating the Dimension of a Model," Annual Statistics, Vol. 6, No. 2, pp.461--464, 1978.

これは次式により定義されます．

I(\hat{\mathcal{S}})\equiv -2\log L(\hat{\mathcal{S}};\mathcal{X})+q\log N

第1項の$L(\hat{\mathcal{S}};\mathcal{X})$は「標本群$\mathcal{X}$を出力するシステムのモデルとしての$\hat{\mathcal{S}}$の尤もらしさ」を意味するもので， 尤度（ゆうど） （likelihood）と呼ばれます．次の数式で定義されます．

L(\hat{\mathcal{S}};\mathcal{X})
\overset{\mathrm{def}}{=}P(\mathcal{X}|\forall\boldsymbol{x}_{i}\in\hat{\mathcal{S}})
=\prod_{i=1}^{N}P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}})

ただし$P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}})$は 確率密度関数 と呼ばれ，「モデル$\hat{\mathcal{S}}$で近似されるシステムから$\boldsymbol{x}$が観測される確率」を意味します．対数をとっている理由および2が掛けられている理由は後ほど分かります．
第2項の$q$はシステムパラメータの数です．ある標本が観測される確率が常に1より小さいことを考えれば，同一のモデルに対し$N$が多いほど尤度は下がります．また$q$が多いほどより多くの標本に適合するモデルとなり得ますが，それは同時にモデルの過剰適合（標本が観測される確率のみ高め，他のモデル要素が観測される確率を下げてしまい，モデルの汎化能力を妨げること）を引き起こします．BICの第2項は，$N$と$q$に対しペナルティをかけている，と定性的に解釈できます．

ここで，$P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}})$が標準正規分布（ガウス分布）

P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}})
=\frac{1}{\sqrt{(2\pi)^{m}|\boldsymbol{\Sigma}|}}\exp\left\{-\frac{1}{2}
\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})^{\mathrm{T}}
\boldsymbol{\Sigma}^{-1}
\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})
\right\}

で近似できると仮定します．ただし$\boldsymbol{\Sigma}$は次式で得られる誤差の 分散共分散行列 （variance-covariance matrix）です．

\boldsymbol{\Sigma}=\sum_{i=1}^{N}\frac{\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{\mathrm{T}}}{N}

このとき，$q=p$であることに注意すれば，BICは次のようになります．

\displaylines{
I(\hat{\mathcal{S}})
=-2\log\prod_{i=1}^{N}P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}})+p\log N
\\
=-2\sum_{i=1}^{N} \log
\frac{1}{\sqrt{(2\pi)^{m}|\boldsymbol{\Sigma}|}}\exp\left(-\frac{1}{2}\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{e}_{i}\right)
+p\log N
\\
=-2\sum_{i=1}^{N}\left\{
-\frac{1}{2}\log (2\pi)^{m}|\boldsymbol{\Sigma}|
-\frac{1}{2}\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{e}_{i}
 \right\}
+p\log N
\\
= N ( m\log 2\pi +\log |\boldsymbol{\Sigma}| )
+\sum_{i=1}^{N}\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{e}_{i}
+p\log N
}

対数をとることで，乗算されていた確率密度関数を加法的に扱えるようになること，および指数関数をキャンセルできることがお分かりになるかと思います．ちなみに$\log L(\hat{\mathcal{S}};{\boldsymbol{x}_{i}})$は 対数尤度 （log-likelihood）と呼ばれます．
また，2を掛けていたことで，計算の過程で各項に現れる$\frac{1}{2}$もキャンセルされています．

原典ではさらに次の近似を用いています．

\boldsymbol{\Sigma}\simeq\sigma^{2}\boldsymbol{1}

ただし

\sigma^{2}=\sum_{i=1}^{N}\frac{\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{e}_{i}}{N}

です．このときのBICは，

I(\hat{\mathcal{S}})= N ( m\log 2\pi +\log \sigma^{2} + 1 ) +p\log N

となり，計算がさらに簡単化されます．

次に，$\hat{\mathcal{S}}=\hat{\mathcal{S}}_{\mathrm{A}}+\hat{\mathcal{S}}_{\mathrm{B}}$である場合を考えましょう．ただし，$\hat{\mathcal{S}}_{\mathrm{A}}=\left\{\boldsymbol{x}\left|\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{A}})=\boldsymbol{0}\right.\right\}$，
$\hat{\mathcal{S}}_{\mathrm{B}}=\left\{\boldsymbol{x}\left|\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{B}})=\boldsymbol{0}\right.\right\}$とモデル化されているとします．このとき尤度は次のようになります．

L(\hat{\mathcal{S}};\mathcal{X})
=\prod_{i=1}^{N}\left\{
P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}}_{\mathrm{A}})P(\hat{\mathcal{S}}_{\mathrm{A}})
+P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}}_{\mathrm{B}})P(\hat{\mathcal{S}}_{\mathrm{B}})\right\}

ただし，$P(\hat{\mathcal{S}}_{\mathrm{A}})$は「$\mathcal{S}$が$\boldsymbol{x}$を出力するに当たってサブシステム$\mathcal{S}_{\mathrm{A}}$が関与する確率」です．少々分かりにくいですが，平たく言えば「確率$P(\mathcal{S}_{\mathrm{A}})$でサブシステム$\mathcal{S}_{\mathrm{A}}$が働き，さらに確率$P(\boldsymbol{x}|\boldsymbol{x}\in\mathcal{S}_{\mathrm{A}})$で$\boldsymbol{x}$を出力する」ということです．$P(\hat{\mathcal{S}}_{\mathrm{B}})$も同様です．これらは$P(\hat{\mathcal{S}}_{\mathrm{A}})+P(\hat{\mathcal{S}}_{\mathrm{B}})=1$を満たします．このままBICを計算することは出来ないのですが，原典では次の仮定をおいています．

• $P(\hat{\mathcal{S}}_{\mathrm{A}})=|\mathcal{C}_{\mathrm{A}}|/N$，$P(\hat{\mathcal{S}}_{\mathrm{B}})=|\mathcal{C}_{\mathrm{B}}|/N$である．ただし$\mathcal{C}_{\mathrm{A}}$，$\mathcal{C}_{\mathrm{B}}$はそれぞれ$\hat{\mathcal{S}}_{\mathrm{A}}$，$\hat{\mathcal{S}}_{\mathrm{B}}$に対応するクラスタ．
• $P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}}_{\mathrm{A}})$，$P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}}_{\mathrm{B}})$が次のように近似できる．

\displaylines{
P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}}_{\mathrm{A}})
=\frac{1}{\sqrt{(2\pi)^{m}\sigma^{2}}}\exp\left\{-\frac{1}{2\sigma^{2}}
\hat{f}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{A}})^{\mathrm{T}}\hat{f}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{A}})
\right\}
\\
P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}}_{\mathrm{B}})
=\frac{1}{\sqrt{(2\pi)^{m}\sigma^{2}}}\exp\left\{-\frac{1}{2\sigma^{2}}
\hat{f}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{B}})^{\mathrm{T}}\hat{f}(\boldsymbol{x};\boldsymbol{\theta}_{\mathrm{B}})
\right\}
}

これらの下で，次のように（強引に）BICを導くことが出来ます（システムパラメータの数が$2p$に増えていることにご注意下さい）．

I(\hat{\mathcal{S}})
=
N(m\log2\pi+\log\sigma^{2}-2\log 2)
+2p\log N
+\frac{1}{\sigma^{2}}\left(
\sum_{i=1}^{N_{\mathrm{A}}}\boldsymbol{e}_{i,\mathrm{A}}^{\mathrm{T}}\boldsymbol{e}_{i,\mathrm{A}}
+\sum_{i=1}^{N_{\mathrm{B}}}\boldsymbol{e}_{i,\mathrm{B}}^{\mathrm{T}}\boldsymbol{e}_{i,\mathrm{B}}
\right)

ただし，$\boldsymbol{e}_{i,\mathrm{A}}=\hat{\boldsymbol{f}}(\boldsymbol{x}_{i};\boldsymbol{\theta}_{\mathrm{A}})$，$\boldsymbol{e}_{i,\mathrm{B}}=\hat{\boldsymbol{f}}(\boldsymbol{x}_{i};\boldsymbol{\theta}_{\mathrm{B}})$とそれぞれおきました．

クラスタ分割によって全体の分散を低減する主旨の方法でありながら，システム分割後も分散が維持されるという近似を置いているのは，筆者にとっては，不可解です．後に計算例を示しますが，$X$-meansの性能はそれほど高くありません．恐らくこの強引な仮定が原因ではないかと想像します．

混合ガウス分布（EMアルゴリズム）による同定

EMアルゴリズム

全標本を陽にクラスタに分けるのとは異なるアプローチを考えましょう．$K$個のサブシステム$\mathcal{S}_{k}$（$k=1,\cdots,K$）の直和であるシステム$\mathcal{S}$から標本$\boldsymbol{x}$が観測される確率は，次式のように表せます．

P(\boldsymbol{x}|\boldsymbol{x}\in\mathcal{S})
=\sum_{k=1}^{K}P(\boldsymbol{x}|\boldsymbol{x}\in\mathcal{S}_{k})P(\mathcal{S}_{k})

ただし，$P(\mathcal{S}_{k})$は「$\mathcal{S}$が$\boldsymbol{x}$を出力するに当たってサブシステム$\mathcal{S}_{k}$が関与する確率」です．$P(\mathcal{S}_{k})$は次式を満たします．

\sum_{k=1}^{K}P(\mathcal{S}_{k})=1
\tag{1}

したがって，観測された標本群$\mathcal{X}$に対するモデル$\hat{\mathcal{S}}$の尤度は次式で表せます．

L(\hat{\mathcal{S}};\mathcal{X})
\overset{\mathrm{def}}{=}
\prod_{i=1}^{N}P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}})
=\prod_{i=1}^{N}\sum_{k=1}^{K}P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}}_{k})P(\hat{\mathcal{S}}_{k})

簡単のため，以降は$P(\boldsymbol{x}_{i}|\boldsymbol{x}_{i}\in\hat{\mathcal{S}}_{k})=p_{ik}$，$P(\hat{\mathcal{S}}_{k})=\pi_{k}$とそれぞれ表記することにします．対数尤度は次のようになります．

\log L(\hat{\mathcal{S}};\mathcal{X})
=\log\prod_{i=1}^{N}\sum_{k=1}^{K}\pi_{k}p_{ik}
=\sum_{i=1}^{N}\log\sum_{k=1}^{K}\pi_{k}p_{ik}

$\sum_{k=1}^{K}\pi_{k}=1$の制約の下で上記の対数尤度を最大にする$\hat{\mathcal{S}}_{k}$が，$\mathcal{S}_{k}$のモデルとして尤もらしいと言えます．このような考え方に基づいて$\hat{\mathcal{S}}_{k}$を求める方法を 最尤推定法 （maximum likelihood estimation method）と呼びます．具体的にはラグランジュ乗数$\lambda$を用いて

J=\log L(\hat{\mathcal{S}};\mathcal{X})+\lambda\left(\sum_{k=1}^{K}\pi_{k}-1\right)

を定義し，

\left(\frac{\partial J}{\partial\boldsymbol{\theta}_{k}}\right)^{\mathrm{T}}=\boldsymbol{0}
\qquad\mbox{かつ}\qquad
\frac{\partial J}{\partial\boldsymbol{\Sigma}_{k}}=\boldsymbol{O}
\qquad\mbox{かつ}\qquad
\frac{\partial J}{\partial\pi_{k}}=0

を満たす$\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$を求めれば良いでしょう．$P(\boldsymbol{x}|\boldsymbol{x}\in\hat{\mathcal{S}}_{k})$を標準正規分布で近似し，上式を実際に計算すれば，

\displaylines{
\begin{align}
\left(\frac{\partial J}{\partial\boldsymbol{\theta}_{k}}\right)^{\mathrm{T}}
=\sum_{i=1}^{N}\frac{\pi_{k}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}\left(\frac{\partial p_{ik}}{\partial\boldsymbol{\theta}_{k}}\right)^{\mathrm{T}}
 =-\sum_{i=1}^{N}\frac{\pi_{k}p_{ik}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}
\left(\frac{\partial\boldsymbol{e}_{ik}}{\partial\boldsymbol{\theta}_{k}}\right)^{\mathrm{T}}
\boldsymbol{\Sigma}_{k}^{-1}\boldsymbol{e}_{ik}
 =\boldsymbol{0}
\tag{2}
\\
\frac{\partial J}{\partial\boldsymbol{\Sigma}_{k}}
=\sum_{i=1}^{N}\frac{\pi_{k}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}
\frac{\partial p_{ik}}{\partial\boldsymbol{\Sigma}_{k}}
 =-\frac{1}{2}\sum_{i=1}^{N}\frac{\pi_{k}p_{ik}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}
\boldsymbol{\Sigma}_{k}^{-1}\left(
\boldsymbol{1}-\boldsymbol{\Sigma}_{k}^{-1}\boldsymbol{e}_{ik}\boldsymbol{e}_{ik}^{\mathrm{T}}
\right)
 =\boldsymbol{O}
\tag{3}
\\
\frac{\partial J}{\partial\pi_{k}}
=\sum_{i=1}^{N}\frac{p_{ik}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}+\lambda=0
\tag{4}
\end{align}
}

となります．なお，任意の行列$\boldsymbol{A}$およびベクトル$\boldsymbol{x}$に関する公式

\boldsymbol{x}^{\mathrm{T}}\boldsymbol{A}\boldsymbol{x}\equiv{\mathrm{Tr}}\boldsymbol{A}\boldsymbol{x}\boldsymbol{x}^{\mathrm{T}},
\qquad
\frac{\partial|\boldsymbol{A}|}{\partial\boldsymbol{A}}\equiv|\boldsymbol{A}|(\boldsymbol{A}^{-1})^{\mathrm{T}},
\qquad
\frac{\partial\left({\mathrm{Tr}}\boldsymbol{A}^{-1}\right)}{\partial\boldsymbol{A}}\equiv-(\boldsymbol{A}^{-2})^{\mathrm{T}}

および$\boldsymbol{\Sigma}_{k}$が対称行列である事実を用いました．

上の三式に共通して出てくる

\gamma_{ik}\overset{\mathrm{def}}{=}\frac{\pi_{k}p_{ik}}{\sum_{j=1}^{K}\pi_{j}p_{ij}}
\tag{5}

は，「$\boldsymbol{x}_{i}$が観測されたときに，それを出力したのが$\mathcal{S}_{k}$である事後確率」を意味し， 負担率 （contribution ratio）と呼ばれます．これは次式を満たします．

\sum_{k=1}^{K}\gamma_{ik}=1
\tag{6}

これを用いれば式（2），式（3）はそれぞれ

\displaylines{
\begin{align}
\sum_{i=1}^{N}\gamma_{ik}\boldsymbol{e}_{ik}^{\mathrm{T}}\boldsymbol{\Sigma}_{k}^{-1}\frac{\partial\boldsymbol{e}_{ik}}{\partial\boldsymbol{\theta}_{k}}=\boldsymbol{0}
\tag{7}
\\
\sum_{i=1}^{N}\gamma_{ik}\left(
\boldsymbol{1}-\boldsymbol{\Sigma}_{k}^{-1}\boldsymbol{e}_{ik}\boldsymbol{e}_{ik}^{\mathrm{T}}
\right)
=\boldsymbol{O}
\tag{8}
\end{align}
}

となります．

【例1】 $\hat{\boldsymbol{f}}_{k}(\boldsymbol{x};\boldsymbol{\theta}_{k})=\boldsymbol{x}-\boldsymbol{\theta}_{k}$の場合，式（7）は，

\sum_{i=1}^{N}\gamma_{ik}(\boldsymbol{x}_{i}-\boldsymbol{\theta}_{k})^{\mathrm{T}}\boldsymbol{\Sigma}_{k}^{-1}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{\theta}_{k}=\frac{1}{N_{k}}\sum_{i=1}^{N}\gamma_{ik}\boldsymbol{x}_{i}
\tag{9}

（ただし$N_{k}\overset{\mathrm{def}}{=}\sum_{i=1}^{N}\gamma_{ik}$）となります．

【例2】 $\hat{\boldsymbol{f}}_{k}(\boldsymbol{x};\boldsymbol{\theta}_{k})=\boldsymbol{\theta}_{k}^{\mathrm{T}}\boldsymbol{x}-1$の場合，式（7）は，

\sum_{i=1}^{N}\gamma_{ik}(\boldsymbol{\theta}_{k}^{\mathrm{T}}\boldsymbol{x}_{i}-1)\boldsymbol{x}_{i}^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{\theta}_{k}=\left(\sum_{i=1}^{N}\gamma_{ik}\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{x}_{i}\right)^{-1}\sum_{i=1}^{N}\gamma_{ik}\boldsymbol{x}_{i}
\tag{10}

となります．

式（8）は

\boldsymbol{\Sigma}_{k}=\frac{1}{N_{k}}\sum_{i=1}^{N}\gamma_{ik}\boldsymbol{e}_{ik}\boldsymbol{e}_{ik}^{\mathrm{T}}
\tag{11}

また，式（4）は

\frac{\partial J}{\partial\pi_{k}}
=\sum_{i=1}^{N}\frac{\gamma_{ik}}{\pi_{k}}+\lambda=0
\quad\Leftrightarrow\quad
\sum_{i=1}^{N}\gamma_{ik}+\lambda\pi_{k}=0

となります．これと式（1），（6）より，

\sum_{k=1}^{K}\sum_{i=1}^{N}\gamma_{ik}+\sum_{k=1}^{K}\lambda\pi_{k}=0
\quad\Leftrightarrow\quad
\lambda=-N

となるので，結局

\pi_{k}=\frac{N_{k}}{N}
\tag{12}

を得ます．

式（7）（具体的には式（9）や（10））と式（11），（12）から$\boldsymbol{\theta}_{k}$と$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$が求まったかのように思えますが，実際にはそうではありません．これらの式は$\gamma_{ik}$を含み，定義式（5）から明らかに，$\gamma_{ik}$は$\boldsymbol{\theta}_{k}$や$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$を含んでいます．すなわちこれらは$\boldsymbol{\theta}_{k}$と$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$についての非線形連立方程式となっているわけです．直接的に$\boldsymbol{\theta}_{k}$と$\pi_{k}$を求めることはできないので，$\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$の暫定値から$\gamma_{ik}$を推定し（Expectationステップ），それを用いて最尤推定法に基づき$\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$を更新する（Maximizationステップ），という二つのステップを交互に繰り返す方法が考えられます．これは EMアルゴリズム と呼ばれます．

• A. P. Dempster, N. M. Laird, D. B. Rubin, "Maximum Likelihood from Incomplete Data via the EM Algorithm," Journal of the Royal Statistical Society, Series B (Methodological), Vol. 39, No. 1, pp.1--38, 1977.

アルゴリズムは次のように具体化されます．

1. $\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$の適当な初期推定値を定める．
2. 現在の$\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$の推定値を用いて，式（5）より$\gamma_{ik}$を推定する．
3. 現在の$\gamma_{ik}$を用いて，式（7）（11）（12）より$\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$を更新する．
4. 次のいずれかの状況となったら終了．(i) $\boldsymbol{\theta}_{k}$，$\boldsymbol{\Sigma}_{k}$，$\pi_{k}$，$\gamma_{ik}$いずれかの変化が$\varepsilon$以下となった．(ii) 対数尤度の変化が$\delta$以下となった．(iii) 反復回数が一定数を越えた．さもなければ2に戻る．

初期推定値の選び方が計算結果に影響を及ぼしますが，一般的にEMアルゴリズムよりも高速な$K$-means法を用いて得たクラスタからこれらを定めるのが良いとされています．

このようにして得られたモデルは，$\boldsymbol{x}$を出力する確率が標準正規分布（ガウス分布）の重み付き和となるので， 混合ガウス分布モデル （Gaussian Mixture Model, GMM）と呼ばれます．

ガウス分布数はどのように選ぶ？

$K$-means法のときと同様に，$K$が未知の場合は，幾つか異なる$K$についてEMアルゴリズムを実行した後，最も良い結果を与えた$K$を選択することになります．モデル$\hat{\mathcal{S}}$の良し悪しを評価する指標は幾つか提案されていますが，BICがよく使われます．混合ガウス分布においては次のようになります．

I(\hat{\mathcal{S}})=-2\sum_{i=1}^{N}\log\sum_{k=1}^{K}\pi_{k}p_{ik}+Kp\log N

モデルパラメータ数が$Kp$であることにご注意下さい．これについても後ほど計算の実例を示します．

実際の計算例によるK-means法，X-means法，EMアルゴリズムの結果比較

$K$-means法，$X$-means法，EMアルゴリズムをZMに実装し，試してみました．用いたのは次の3種のデータセットです．

(A) 乱数的に選んだ6つの座標それぞれのまわりに一様乱数的に1000点配置，
(B) (A)と同じやり方で点群を配置した後にアウトライア10点を追加配置，
(C) 乱数的に選んだ6つの座標それぞれのまわりに標準正規分布に基づいて乱数的に1000点配置．

まず$K$-means法（$K$-means++法）を試しました．いずれの試行においても，同一の$K$に対し3回の結果を比較して最良のものを採用しました．さらにシルエット分析に基づいて最適な$K$を選択しました．

(A)のケースでは$K=5$，(B)(C)では$K=6$となりました．それほど違和感の無い結果となっています．

次に$X$-meansを試しました．

(B)のケースではクラスタ分割が起こりませんでした．(A)(C)も十分分割されているとは言えない結果になりました．

最後にEMアルゴリズムの結果です．データセットに得られた確率密度関数を重ねて，立体的に可視化します．

ただし$K$は全て6としています．参考までに，$K$を$4\sim 9$と変えて得られたモデルのBICをプロットしたものが次のグラフです．

紫線がBICです．緑線はAIC（赤池情報量基準）と呼ばれる別の指標です（詳細は割愛します）．

• H. Akaike, "Information theory and an extension of the maximum likelihood principle," in Proceedings of the 2nd International Symposium on Information Theory, Akadimiai Kiado, Budapest, pp. 267--281, 1973.

どちらの値も，(A)(B)のケースでは生成クラスタ数$K=6$に対し大きく低下しているものの，最小値とはなっておらず，それより大きい$K$に対して減少していっています．(C)のケースでは減少傾向にあまり特徴がありません．いずれにしても，これを見て$K$を自動的に決定するのは難しそうです．

おまけで，$\hat{\boldsymbol{f}}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{x}-1$（$k=1,2,3$，$n=3$）とし$\boldsymbol{\theta}_{1}$，$\boldsymbol{\theta}_{2}$，$\boldsymbol{\theta}_{3}$を乱数的に作成したときの$K$-means法およびEMアルゴリズムの結果を示します．

ちょっと分かりにくいですが，どちらの方法でも平面をうまく推定できていることが見てとれるかと思います．

まとめ

複合代数的システムのモデル同定は典型的なchicken-and-egg問題であり，標本の所属する部分システムの推定とモデルパラメータの同定を同時に行う必要があります．これには$K$-means法とEMアルゴリズムが応用できます．

使用者観点での両者の違いは，標本をクラスタリングするかしないかです．性能差は，本記事で示したような単純なケースでは大きく無さそうです．部分システムの重複が大きかったり，アウトライアが多かったりするケースではEMアルゴリズムが有利になるだろうと予想は出来ますが，未検証です．

部分システム数$K$が未知の場合，$K$-means法ならばシルエット分析によって自動推定することができます．
$X$-means法は残念ながら信頼できる方法では無さそうです．
EMアルゴリズムにおいては，BICやAICは参考にはなるものの，直接的な自動推定の規範にはなりません．他にうまい自動推定のやり方があるかどうかは十分調査できていませんが，恐らくは$K$-means法で初期推定値を得る際に，シルエット分析でクラスタ数を推定し，それをそのままEMアルゴリズムのクラスタ数として用いるのが良いのでしょう．

ZMでは、$K$-means法と$X$-means法はzm_mca_cluster.cで、EMアルゴリズムはzm_mca_gmm.cでそれぞれ実装しています．シルエット分析に基づくクラスタ数自動決定は、app/zm_kmeansにて実装しています．