ベイズ推論により混合分布のパラメータ推論をやってみる　〜混合数の推論２〜

【概要】

混合分布（混合モデル）はモデルを潜在変数でスイッチする構造を持ったモデルであり、実用的な観点でも面白いです
弊ブログでは数回にわたって、混合分布を使って遊んでみています（これが5記事目）
第4弾、第5弾では、これまで既知としていた混合数（クラスタ数）を含めて推論してみます（記事が長くなったので分割しました）
まだまだ勉強中なので、間違いなど指摘いただけると助かります

【目次】

はじめに
混合モデルとディリクレ過程
ディリクレ過程混合モデルの実装
全体のnotebook
まとめ
参考文献

はじめに

機械学習や統計の問題では、手元にあるデータを解釈して応用しようとしますね。この時、明に暗になんらかの「確率モデル」を仮定しているはずです。

確率モデルの中でも、混合分布（混合モデル）は、複数の確率モデルの組み合わせとして定義されており、複雑なデータ構造を表現できます。応用としても、クラス分類や多クラスの回帰など面白い実用例があります。

当ブログではこれまで複数回にわたって混合モデルを使って遊んできました。第5弾となる本記事では、これまで既知としていた混合数（クラスタ数）もデータに合わせて推論してみたいと思います。ディリクレ過程混合モデル(DPMM; Dirichlet Process Mixture Model)と呼ばれるモデルです。ほんとは一つの記事でまとめたかったのですが、長くなってしまったので、第4弾と合わせて読んでもらえたらと思います

【ベイズ推論により混合分布のパラメータ推論をやってみる記事一覧】

第1弾： learning-with-machine.hatenablog.com
第2弾： learning-with-machine.hatenablog.com
第3弾： learning-with-machine.hatenablog.com
第4弾： learning-with-machine.hatenablog.com
第5弾：当記事

【トップに戻る】

混合モデルとディリクレ過程

混合モデルとは、複数の確率モデルを組み合わせた構造をしており、潜在変数 $s$ でそれらの確率モデルをスイッチする構造を持った確率モデルです。

f:id:hippy-hikky:20200308143035p:plain

前回の記事で説明した通り、上図の $\pi$ が潜在変数 $s$ の種類数と生成され方を表現したものです。この $\pi$ の次元数などでクラスタの数を固定しているのが通常の混合モデルです。

混合数（クラスタ数）を固定すると、現実の問題では柔軟さが足りない場合があります。そのため、無限の混合モデルを考えるのがディリクレ過程混合モデルでした（詳しくは前回記事や参考資料2, 3などを確認してください）。

ディリクレ過程を実現するために二つのアプローチがあり、対応する二つの代表的な確率過程があります。

中華料理店過程(CRP; Chinese Restaurant Process)
- $\pi$ を周辺化除去して、 $s_n$ について無限クラスの分割方法に関する事前分布を考える
棒折り過程(SBP; Stick Breaking Process)
- 無限次元の $\pi$ を表現する事前分布を考える

前回の記事では、この二つの確率過程を確認するところまで実装しました。本記事では実際に混合モデルの推論に適用してみたいと思います。

【トップに戻る】

ディリクレ過程混合モデルの実装

ディリクレ過程混合モデル(DPMM)を実装してみます。

DPMMによるデータの生成プロセス

ディリクレ過程の実現例として、CRPとSBPを紹介しました。実際にこれらを利用して混合モデルのデータ生成過程（モデル）を構成すると次の通りとなります。

CRPを利用した混合モデルのデータ生成プロセス

f:id:hippy-hikky:20200826162819p:plain:h100

CRPはデータの分割に関する事前分布( $P(S)$ )となるのでした。そのため、潜在変数（クラスタ番号）の生成に使われます(1)。それぞれのコンポーネントのパラメータ（例えば平均と分散）を $\theta$ とすると、 $\theta$ はなんらかの基底分布 $G_0(\theta)$ から生成されます(3)。そして、それぞれのデータ $x_k$ は $s_k=i, \theta_i$ に基づいて生成されます(3)。

通常の混合モデルにおいてsの生成がCRPからされると見れば良さそうです。

SBPを利用した混合モデルのデータ生成プロセス

次にSBPを利用した場合のデータ生成プロセスです。

f:id:hippy-hikky:20200826164253p:plain:h120

SBPの場合、無限次元の混合比を生成する確率過程ですが、計算上、上限の値(K)を定めています(1)。コンポーネント毎のパラメータについてはCRPと同様に生成します(2)。上記の式では、潜在変数 $s$ を周辺化除去しています。そのため、(3)では混合比とK個の確率モデルの混合モデルからxが生成されています。（SBPの場合、sを積分消去しないといけないのかについてはちゃんと調べていないです。参考資料[1]の記述に従っています。）

SBPでは、無限次元（計算上は有限次元）の混合比が生成されるとみなせば良さそうです。

実装方針

DPMMの実装については、文献2では、CRPを利用してギブスサンプリングによってパラメータの推論をしています。実際、混合モデルのグラフィカルモデルにおいて、確率変数 $s$ の事前分布として利用するCRPはわかりやすいと思います。クラスタリングタスクにおいても、潜在変数sが推論できればそのままクラスタ番号として利用できますので。

しかしここでは、MCMC アルゴリズムによって事後分布からのサンプルを取得して推論する方法を取りたいと思います。

MCMC アルゴリズムを利用する場合、混合モデルにおける離散確率変数である $s$ を積分消去することで計算効率と計算の安定性が得られるとされています（この記事やこちらなど）。そのため、混合比 $\pi$ の事前分布となるSBPが相性が良いです(参考資料1より)。

ということで、SBPを利用してDPMMを構成し、MCMC アルゴリズムを利用してサンプルを取得してみようと思います（参考資料1を参考にしてください）。なお本記事では、SBPのパラメータである $\alpha$ は固定の値として、推論はしないものとします。また、実装はPPLの一つであるPyMC3を利用します(PyMC3については参考資料7など)。

問題設定

本記事では、1次元のガウス混合モデルを考えてみます。

f:id:hippy-hikky:20200826175813p:plain — サンプルデータのヒストグラム．3つのコンポーネントからなる１次元ガウス混合モデル．実線はKDEによる密度関数．

上図のデータに対して、コンポーネント数は不明として、ガウス混合モデルのパラメータを推論します。

実装結果

実装については、添付のnotebookを参照してください。サンプルデータの生成から推論結果の可視化まで全てのコードを載せています。重要なのは以下の部分です。

これはPyMC3において確率モデルを定義する部分です。PyMC3ではこのように先に定義した確率モデルをほぼそのまま列挙していくだけでモデルの定義ができます。この中でpm.Mixtureを使っていますが、これはPyMC3で混合モデルを定義するAPIです。クラスタの割り当てを示す確率変数 $s$ は周辺化するため、パラメータとして指定していないことに注意です。

with pm.Model() as model:
    vi = pm.Beta('vi', 1., alpha, shape=K)
    pi = pm.Deterministic('pi', SBP(vi))
    
    comp_dist = []
    for k in range(K):
        mu = pm.Normal('mu%i'%k, mu=0, sd=10)
        tau = pm.Gamma('tau%i'%k, alpha=1., beta=0.1)
        comp_dist.append(pm.Normal.dist(mu=mu, tau=tau))
    obs = pm.Mixture('obs', pi, comp_dist, observed=x)

結果は次の通りです。まずは、混合比の事後分布を確認してみます。

f:id:hippy-hikky:20200826180902p:plain — 混合比の事後分布の平均．横軸はコンポーネントを示す．

この図から、設定通り3つ程度のコンポーネントに意味があり、以降の割合は小さいことがわかります。

次に、割合が上位のコンポーネントのパラメータのMCMCサンプル系列を確認します。

f:id:hippy-hikky:20200826181207p:plain — 混合比が5%以上のコンポーネントのパラメータのトレース．上が平均で下は精度(分散の逆数)．

この結果から、上位3つのコンポーネントは設定値に近い値で収束していることがわかります。4番目のコンポーネントは収束せずデータ全体を彷徨っていますね。

最後に、この事後分布を利用してクラスタリングをした結果を示します。クラスタリングについては、sを周辺化除去しているため、各コンポーネントの尤度を計算し、尤度の割合からクラスタの推定をしています（添付のnotebook参照）。

f:id:hippy-hikky:20200826181615p:plain — クラスタリング結果．黒の実線はデータ全体のKDEプロット．破線はコンポーネント毎の確率密度．色は推定クラスタを示す．

図から、各データが3つのクラスタにクラスタリングされていることがわかります。

【トップに戻る】

全体のnotebook

【トップに戻る】

まとめ

今回は、混合モデルにおける混合数の推論まで行うために、ディリクレ過程混合モデルを実装してみました。実装と言っても、推論の実態はPPL（PyMC3）に任せっきりですが。、、

DPMMを実装してみたのですが、私自身まだ数学的に理解しているとは言えないです。なので、本記事は表面的な理解のもとで、PPLを使って実装してみたという内容です。詳しい方、間違いなどありましたら指摘してくれると助かります。

SBPを使ってDPMMを実装したのですが、SBPではコンポーネント数の上限を定めます。ディリクレ過程ではない混合モデルでも、大きい混合数で推論することは可能で、混合比で使うコンポーネントを選択することもできると思います。この違いについてまだ理解ができてないので、詳しい方、教えていただけると助かります。（DPMMではデータが割り当てられないコンポーネントがあることが許されるみたいなことが違うのかな？？）

混合モデルの推論についてのシリーズは一旦これで終了とするつもりです。ディリクレ過程についてはまだまだ理解が浅いので、引き続き勉強を続けていきます。CRPを利用してギブスサンプリングでのDPMMの推論なども実装して記事を書いていこうと思います。

【トップに戻る】

参考文献

[1] Dirichlet process mixtures for density estimation, PyMC3 Examples, 2020/08/20アクセス
[2] 石井, 上田, 続・わかりやすいパターン認識, オーム社, 2014(初版)
- これは実際わかりやすかった。扱ってる内容がハードなので難しいんだけど、それでもわかりやすかったです。これは必読です。