混合ディリクレ分布でクラスタリング

クラスタリングには様々な手法がありますが、目的、データの分布などに合わせて適切なものを選択します(1).
今は強力なツール、ライブラリが揃っていますので大体はそれでまかなえてしまいますし、様々な高度なアルゴリズムをすぐに使うことができます(2).
しかし、アルゴリズムによっては提供されていないこともありますので、実装ができると幅が広がります.

今回紹介する「混合ディリクレ分布」は運悪く、私のよく使う sckit-learn にないため、頑張って実装してみようと思います(私自身、アルゴリズムを実装からするのは久々です..).

混合ディリクレ分布

今回は、「足したら 1 になるデータ」の集合に対する分布表現・クラスタリング手法を考えてみました.
このようなデータには例えば、当社の顧客価値観モデル「Societas」があります.
Societas では人の価値観を 12 の特性パターンの割合(所属確率)で表現します.具体的には例えばある人の価値観を、

[パターン 2-2 家庭的な真面目タイプ]:60%
[パターン 2-1 家族大好き悠々タイプ]:15%
[パターン 4-1 自分中心的なアクティブタイプ]:10%
[パターン 5-2 社交的な堅実ホームメーカータイプ]:10%
...

のように表現しますが、こういったデータに分布にあてはめる場合です.
例えば、上の例のソシエタスの割合からできる価値観イメージは「実用性を重視した、しっかり派」で、下図のようなパイチャートで表現したりします.
societas_cluster

さて、このようなデータの分布モデルとしてはディリクレ分布があり、「サイコロの出目のでやすさ」の確率分布と言われ、いびつなサイコロで確率分布を表現するイメージになります(=ある面が出やすかったり出にくかったりしますが、出目の確率の和は 1 ).

今回はこのディリクレ分布を複数重ねあわせてデータ分布を表現します.
と言っても、実際には混合ガウス分布モデルのガウス分布をディリクレ分布に置き換えるだけです(3).

アルゴリズム

まず、データの分布を複数のディリクレ分布の重ね合わせで定式化します. α はディリクレ分布パラメータ、π は混合比と呼ばれる潜在変数、k は重ね合わせる分布の数です.

p(x) = \sum_{k=1}^K \pi_k \ dirichlet(x|\alpha_k)

データの分布を表現できたらその尤度を最大化することで最適なパラメータを見つけることができますが、今回は混合ガウス分布と同様に EM アルゴリズムで見つけます.
異なるところは、M ステップでディガンマ関数の逆関数をニュートン法で解くところだけです.

  1. E ステップ
    各混合要素の負担率を計算します.

    r_{nk} = \frac {\pi_k \ p(x_n | \theta_k)} {\sum_{c=1}^K \ \pi_c\ p(x_n | \theta_c)}

  2. M ステップ
    現在の負担率を使って、パラメータ値を再計算します。
    Q 関数は混合ガウス分布と同様ですが、ガウス分布をディリクレ分布と置き換えます。

    Q(\theta, \theta^{old}) = \sum_{n=1}^N \ \sum_{k=1}^K r_{nk}\{ log\ \pi_k \ dirichlet(x_n | \alpha_k) \}

    これは凸関数ですのでパラメータで偏微分し、最大化することで新しいパラメータと置き換えますが、

    \alpha_{kd} = \psi^{-1} \bigl( \psi(\sum_{d=1}^D \alpha_{kd}) + \frac {\sum_{n=1}^N r_{nk} \log x_{nd}}{\sum_{n=1}^N r_{nk}} \bigr)

    この時、ディガンマ関数の逆関数をニュートン法で解きます(4).

    x^0 = \psi^{-1}(y) = \left\{\begin{aligned}    exp(y) + 1/2 &: if\ y \geq -2.22 \\    -1/(y-\psi(1)) &: if\ y < -2.22\end{aligned}\right. \\x^{t+1} = x^{t} - \frac{\psi(x) - y}{\psi^{

  3. E と M のステップを繰り返し、対数尤度が収束したときに得られる α、π が最適なパラメータ、潜在変数となります.

実際にやってみる

さて、この混合ディリクレ分布の EM アルゴリズムを実際に次の 3 次元データに適用してみます.
本当は Societas の分類数と同じ 12 次元でやってみたいのですが、各イテレーションでの分布を可視化したいですので 3 次元データにします.
データは足したら 1 ですので次のようなデータの集合になります.

本来は 3 次元ですが、足したら 1 という制約がついていますので、2 次元上にマッピングできて、下図のような三角形領域内に分布します.

test_data

今回は 3 つのディリクレ分布から生成した、こんな感じのデータをテスト分布として、前述のアルゴリズムでこの分布パラメータが推定できるかどうかを試してみます.

dmm_test_data_dist

初期値は適当に設定して、各イテレーションでの推定分布の推移を見てみます.

dmm_iteration

赤点が推定分布の中心点(ディリクレ分布の期待値)ですが、イテレーションを繰り返すにつれ、テストデータに近づいてゆきます.
複数の分布の境界領域は微妙ですが(この領域の取り合いに分布モデルの個性がでそうですね.)、なんとなく(汗)元データの生成分布を推定できていそうです.

イテレーション毎の対数尤度の推移もプロットします.100 イテレーションくらいで収束に近づきます.

dmm_loglikelihood

終わりに

今回は簡単な実装とデータでやってみましたが、なんとなく混合分布モデルのイメージがつかめたと思います。
実際にはデータサイズや運用など諸々総合的に判断して(5)、実績のある機械学習ライブラリのアルゴリズムを採用することがほとんどでしょうが、実装すると理解が深まり適切なアルゴリズムの選択やパラメータ探索など役立つことも多いのではと思います.

分布の混合と潜在的変数はトピックモデルへと続く潜在的意味解析の入り口ですので、機会があればこの道を辿ってみたいと思います.

参考ですが、私は Emacs を使っており iimage-mode でよく関数のコメントに数式や手描きのポンチ絵の画像を差し込んだりします.後々、見なおしたりする際に、なかなか便利です.

iimage_mode

<脚注>
1) 例えば K-means や 混合ガウスモデルは、データ変数がカテゴリカルな場合には適切ではなかったりします(ユークリッド距離が適切ではないため)
2) 利点の反面、「ツールのできること > 自分のできること」、となってしまう傾向にあります.
3) パターン認識と機械学習 下 (ベイズ理論による統計的予測) 9.2.2章
4) Estimating a Dirichlet distribution. 数イテレーションで 10-14 桁精度になるそうです.
5) 最適なクラスタ数の探索や、尤度関数の発散する特異性に対する対策、等々

Comments are closed, but you can leave a trackback: Trackback URL.