Dirichlet process mixture model (以下 DPMM)。"ノンパラメトリック(無限次元)ベイジアンクラスタリング" と呼ばれたりもする。クラスタ数を事前に指定する必要がない(モデルとハイパーパラメタによって決まる)。
DPMMは比較的有名なモデルのため、和文洋文含めて数多くの解説記事・論文が存在する。しかし自分のレベルに合った一連の説明・記法が欲しかったのでまとめた。実態は参考文献に挙げた論文や本、資料の換言や貼り合わせでしかない。前提知識はベイズ推論に対する軽い(真下にある最初の概観があーそなんだくらいの(？))理解だろうか。
DPMMをベイズ推論の枠組みで言うと
ベイズ推論は、観測変数(データ) $D$ および潜在変数(パラメタ) $\theta$ からなる(自分で作った)生成モデル $p(D,\theta)$ から事後分布 $p(\theta|D)$ を求めること。事後分布の最大値を与えるパラメタが MAP 推定値。
$$
p(\theta|D)=\frac{p(D,\theta)}{p(D)}=\frac{p(D|\theta)p(\theta)}{p(D)}
$$
DPMMでは、
- 生成モデル = [データ分割の確率分布] $\times$ [データ分割が与えられたときのデータの生成確率] となる
- 前者がディリクレ過程(DP; ただし実際の推論計算ではそれと等価な中華料理店過程が使われる)
- 後者が混合分布(MM)
- 事後分布は解析的に与えられない(これはディリクレ過程のせい)ので、MCMC で求める
ディリクレ過程 ― データの分割に関する生成モデル
元々1973年に考案されていたが、1990年代に MCMC による効率的な事後分布推定方法が可能となってから一気に実用性が上がった(cf. 下の「参考文献」)。
DPMMからの視点だけで言ってしまうと、ディリクレ過程は「クラスタ数、および各クラスタが持つパラメタの事前分布の確率分布(= 生成モデル)」となる。特に重要なのは前者のクラスタ数である。
ディリクレ「分布」の生成モデル
ディリクレ分布は、正の実数値 $\alpha_k\ (k=1,\cdots,K)$ が与えられた時に、$\pi_k\geq 0,\sum_k \pi_k=1$ なる確率変数の同時分布を与える。
$$
{\rm Dir}(\pi_1,\cdots,\pi_{K}|\alpha_1,\cdots,\alpha_{K})=C\prod_{k=1}^{K}\pi_{k}^{\alpha_{k}-1}
$$
($C$ は正規化項)。これは直感的に言うと、「各離散カテゴリの頻度 $\alpha_k$ が(パラメタとして)与えられた時の、各カテゴリの生起確率 $\pi_k$」についての分布を与える。実際にディリクレ分布はカテゴリ分布および多項分布の共役事前分布となっていて、例えばサイコロの $K$ 個の出目のそれぞれの観測頻度 $n_k$ が
$$
{\rm Multi}(n_1,\cdots,n_K|\pi_1,\cdots,\pi_K)=\prod_{k=1}^{K}\pi_k^{n_k}
$$
に従うとした場合には、事後分布を求める(=ベイズ推論する)ことは事前情報 $\alpha_k$ をもとに観測頻度 $n_k$ を使って生起確率 $\pi_k$ を更新することに相当する。
$$
p(\pi_1,\cdots,\pi_K|n_1,\cdots,n_K,\alpha_1,\cdots,\alpha_K)={\rm Dir}(\pi_1,\cdots,\pi_{K}|\alpha_1+n_1,\cdots,\alpha_{K}+n_K)
$$
上の記法と対応づけると、ディリクレ–多項分布の生成モデルは以下のようになる。
p(\theta)=p(\pi_1,\cdots,\pi_{K})={\rm Dir}(\pi_1,\cdots,\pi_{K}|\alpha_1,\cdots,\alpha_{K}) \\
p(D|\theta)=p(n_1,\cdots,n_K|\pi_1,\cdots,\pi_{K})={\rm Multi}(n_1,\cdots,n_K|\pi_1,\cdots,\pi_K)
ディリクレ「分布」混合モデルの生成モデル
データの観測列を ${x_i}$、$i$ 番目のデータが属する分布(クラスタ)のIDを $z_i$ 、$k\ (=1,\cdots,K)$ 番目の分布(クラスタ)が持つパラメタを $\theta_k$ とすると、混合分布の生成モデルは以下のようになる。
p(\pi_1,\cdots,\pi_{K})={\rm Dir}(\pi_1,\cdots,\pi_{K}|\alpha_1,\cdots,\alpha_{K})
\\
p(z_i=k|\pi_1,\cdots,\pi_{K})=\pi_k
\\
p(x_i|z_i=k,\theta_1,\cdots,\theta_K)=p(x_i|\theta_k)
ディリクレ過程混合モデルの生成モデル
上の生成モデルで、(1つの)クラスタパラメタ$\theta$ が従う分布 $p(\theta)$ の事前分布を $G_0$ とする。このとき、ディリクレ過程は「$\theta$ が従う(離散)確率分布 $G$ の確率分布」となる。具体的には、$G$ がディリクレ過程に従うとは、$\theta$ が取りうる値からなる標本空間の任意の分割 $(T_1,\cdots,T_n)$ に対して、
$$
p(G(T_1),\cdots,G(T_n))={\rm Dir}(G(T_1),\cdots,G(T_n)|\alpha G_0(T_1),\cdots,\alpha G_0(T_{n}))
$$
となることをいう。ただし、$G=G(\theta)$ は $\theta$ の生成分布を表すのに対し、$G(T_i)$ は $G$ から生成された $\theta$ が $T_i$ に属する確率を表すことに注意。$G_0$ は基底関数と呼ばれ、$G$ は $E[G(T_i)]=G_0(T_i)$ となるような離散分布となる。そして、この離散分布 $G(\theta)$ の非零値の個数がクラスタ数(の事前分布)となる。$\alpha$ は集中度と呼ばれ、$\alpha$ の値が大きいほど $G(T_i)$ は平均値 $G_0(T_i)$ に集中する(そしてクラスタ数が多くなる; cf. 後述のCRP)。
(http://chasen.org/~daiti-m/paper/ieice10npbayes.pdf より引用)
ディリクレ過程混合モデルの生成モデルは以下のようになる。この定式化ではデータ1つ1つに対してクラスタパラメタを割り当てていることになるが、このままでは効率的な計算ができない。
p(G(\theta))={\rm DP}(\alpha,G_0(\theta))
\\
p(\theta_i|G)=G(\theta)
\\
p(x_i|\theta_i)=p(x_i|\theta_i)
そこで、上の「ディリクレ分布による混合モデル」における生成モデルと同じく assignment variable $z_i$ を使った定式化をするために、ディリクレ過程の代わりに、それと等価な中華料理店過程(Chinese Restaurant Process; CRP)を用いる。
(http://proc-cpuinfo.fixstars.com/2017/06/clustering-by-dirichlet-process-mixture-model-1/ より引用)
CRP を使った DPMM の生成モデルは以下のようになる。ただし、$n_k$ は $z_1,\cdots,z_{i-1}$ のうち値が $k$ であるものの個数(= $z_{i-1}$ まででのクラスタ $k$ のサイズ)、 $\tilde{K}$ は $z_{i-1}$ まででのクラスタ数を表す。
p(z_i=k|z_1,\cdots,z_{i-1})={\rm CRP}(\alpha)=
\begin{cases}
\frac{n_k}{i-1+\alpha}\ \ \ (k=1,\cdots,\tilde{K}) \\
\frac{\alpha}{i-1+\alpha}\ \ \ (k=\tilde{K}+1)
\end{cases}
\\
p(\theta_k)=G_0(\theta)
\\
p(x_i|z_i=k,\theta_1,\cdots,\theta_K)=p(x_i|\theta_k)
クラスタ数 $K$ は理論上無限個が可能だが、実際は高々データ数以下だけを考えればよく、さらに実用上は CRP で生成されたクラスタの数をそのまま使えば良い。
そして、CRP が持つ重要な性質が交換可能性である。これは ${z_i}$ の同時確率
\begin{eqnarray}
p(z_1,\cdots,z_N)&=&p(z_1)p(z_2|z_1)\cdots p(z_N|z_1,\cdots,z_{N-1}) \\
&=& \frac{\alpha^K\prod_{k=1}^{K}(n_k-1)!}{\alpha(\alpha+1)\cdots(\alpha+N-1)} \\
&=&{\rm Ewens}(n_1,\cdots,n_K;\alpha, N)
\end{eqnarray}
が観測データ ${x_i}$ の並び順によらないことを意味する。これにより以下で述べる Gibbs sampling が可能となる。一方で CRP だと変分推論ができない(変分推論には棒折り過程などが使われる)。
DPMMのベイズ推定
観測データの数を $N$ とする。結局、CRPによるDPMMの事後確率は以下のようになる。
$$
p(z_i,\theta_k)\propto {\rm Ewens}(n_1,\cdots,n_K;\alpha, N)\prod_{k=1}^{K}G_0(\theta_k)\prod_{{x_i|z_i=k}}p(x_i|\theta_k)
$$
これは $G_0$ のせいで明らかに解析的に解けないので、MCMC 等に頼ることになる。$G_0(\theta_k)$ が $p(x_i|\theta_k)$ の共役事前分布である場合には、最もシンプルには (collapsed) Gibbs sampling による $z_i$ の更新
$$
z_i\sim p(z_i|z_{-i},x,\theta)
$$
が使える。これを計算すると
\begin{eqnarray}
p(z_i=k|z_{-i},x,\theta)&=&p(z_i=k|x_i,z_{-i},x_{-i},\theta) \\
&=&\frac{p(z_i=k,x_i|z_{-i},x_{-i},\theta)}{p(x_i|z_{-i},x_{-i},\theta)} \\
&\propto&p(z_i=k,x_i|z_{-i},x_{-i},\theta) \\
&=&p(z_i=k|z_{-i})p(x_i|z_i=k,\theta)
\end{eqnarray}
前半部分はCRPとその交換可能性($z_i$ を最後に持ってくる)により
p(z_i=k|z_{-i})=
\begin{cases}
\frac{n_{-i,k}}{N-1+\alpha}\ \ \ ({\rm if}\ k=z_{i^′}\ {\rm for\ some}\ i^′\neq i) \\
\frac{\alpha}{N-1+\alpha}\ \ \ ({\rm if}\ k\neq z_{i^′}\ {\rm for\ all}\ i^′\neq i)
\end{cases}
後半部分は、$x_i$ を既存クラスタに割り当てる場合は既存のクラスタパラメタをそのまま使い、新しいクラスタを作る場合は可能なパラメタで積分消去したものを考えるとこの積分は $G_0$ が共役事前分布なので解析的に解くことができる。
p(x_i|z_i=k,\theta)=
\begin{cases}
p(x_i|\theta_k)\ \ \ ({\rm if}\ k=z_{i^′}\ {\rm for\ some}\ i^′\neq i) \\
\int G_0(\theta_k)p(x_i|\theta_k)d\theta_k\ \ \ ({\rm if}\ k\neq z_{i^′}\ {\rm for\ all}\ i^′\neq i)
\end{cases}
$z$ が更新されると $\theta$ も(共役事前分布を使って)同時に更新される。$G_0$ が共役事前分布でない場合のサンプリングについては [Neal, 2000] などを参照。
ただし、Gibbs sampling だと(DPMMでは特に？)局所最適解から抜け出す確率が極めて低くなりやすい ([Neal, 2004] etc.)。これを解消するために考案されたのが split-merge sampling である[Neal, 2004; Dahl, 2004]。一度に更新する単位を1つのデータではなくクラスタとすることで、局所解に陥りにくくする効果がある。split-merge sampling は Metropolis-Hastings 法の一種である。
具体的には、各更新ステップにおいてクラスタの結合もしくは分割の提案が行われ、それが以下の acceptance ratio に基づいて受理もしくは棄却される。ここでは全 assignment variable の値(これはクラスタリングにほかならない) $z$ が1つの「状態」として扱われている。
a(z\rightarrow z^*)=\min\left[1, \frac{q(z^*\rightarrow z)}{q(z\rightarrow z^*)}\frac{p(z^*)}{p(z)}\frac{L(z^*|x)}{L(z|x)}\right],\ \ \ z^*\in\{z^{split},z^{merge}\}
ただし、$q$ は提案分布、$p$ は状態の事前分布、$L$ は状態の尤度(事後分布 $p(z|x)\propto p(z)L(z|x)$ が目標分布)。それぞれの具体的な中身は [Neal, 2004] に明記されているので省略するが、$G_0$ が共役事前分布であれば $a(z\rightarrow z^*)$ も解析的に計算できる。
また、2つのクラスタが与えられた場合の結合のやり方は1つしかないが、逆に1つのクラスタが与えられた場合の2分割は一般に何通りも考えられる。ランダムな分割では非常に効率が悪いので、これを避けるために例えば resticted Gibbs sampling が行われる。すなわち、起点となった2点をそれぞれ新しいクラスタの初期元として、分割対象のクラスタ内部だけで Gibbs sampling を行い、妥当な分割を提案する[Neal, 2004]。
発展
- 他のサンプリング手法(particle Gibbs sampling、NUTS sampling)
- MCMC 以外の推論(変分ベイズ、A$^*$ search、beam search)
- DP を使った混合モデル以外のモデル(LDA、階層DPMM、無限HMM)
- DP 以外の確率過程(ベータ過程、Pitman-Yor 過程)によるノンパラメトリックベイズ
参考文献
- 持橋, 最近のベイズ理論の進展と応用 (III) ノンパラメトリックベイズ, http://chasen.org/~daiti-m/paper/ieice10npbayes.pdf
- Teh, Y. W. Dirichlet Process. https://www.stats.ox.ac.uk/~teh/research/npbayes/Teh2010a.pdf
- Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics.
- Jain, S. & Neal, R. M. (2004). A Split-Merge Markov Chain Monte Carlo Procedure For the Dirichlet Process Mixture Model. Journal of Computational and Graphical Statistics.
- Dahl, D. B. (2003). An improved merge-split sampler for conjugate Dirichlet process mixture models. Technical Report.
書籍がわりと充実しているのでそちらから入るのが良いかもしれない。
- Murphy, K. P. Machine Learning: A Probabilistic Perspective (Adaptive Computation and Machine Learning) の 25.2節
- 石井・上田, 続・わかりやすいパターン認識―教師なし学習入門続 の11, 12章
- 佐藤, ノンパラメトリックベイズ 点過程と統計的機械学習の数理 (機械学習プロフェッショナルシリーズ)
DPMMの生物学での応用例は、
- Prabhakaran, S., et al. (2014). HIV Haplotype Inference Using a Propagating Dirichlet Process Mixture Model. IEEE/ACM Trans. Comput. Biol. Bioinform.
- Xing, E. P., Sohn, K. A., Jordan, M. I., Teh, Y. W. (2006). Bayesian multi-population haplotype inference via a hierarchical dirichlet process mixture. ICML.
- Duan, T., Pinto, J. & Xie, X. (2018). Parallelized Inference for Single Cell Transcriptomic Clustering with Split Merge Sampling on DPMM Model. BioRxiv.