続・わかりやすいパターン認識ひとり読書会 1 - 渡邊哲司（Satoshi Watanabe）の日記

続パタの核心部である11章～12章（ディリクレ過程混合モデル）、13章（共クラスタリング）を見返しながらより深く理解したいという欲望の下、ひとり読書会と称してまとめていく予定です。

続パタ読書会 1:ディリクレ過程混合モデル
続パタ読書会 2:ディリクレ過程混合モデルのアルゴリズム
 続パタ読書会 3:実験
 続パタ読書会 4:実装
 続パタ読書会 5:事前確率P(s)の妥当性 1
続パタ読書会 6:事前確率P(s)の妥当性 2
続パタ読書会 7:ベル数 1
続パタ読書会 8:ベル数 2
続パタ読書会 9:P(s_k|x_k, s_-k, θ)の導出
 続パタ読書会 10:演習問題12.5の計算
 続パタ読書会 11:演習問題12.6の計算 1
続パタ読書会 12:演習問題12.6の計算 2

使用するテキスト「続・わかりやすいパターン認識」は第1版第3刷のものを使用します。テキストには詳細が書いてない部分も力が及ぶ範囲で埋めていきたいと思います。

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

背景

$d$ 次元特徴空間上に $n$ 個のパターン $\boldsymbol{x}=\{x_1, \ldots , x_n\}$ が分布しており、各パターンは $c$ 種の異なるクラス $\omega_1, \ldots , \omega_c$ に所属しているものとする。ただし観測できるのは特徴空間上の位置に相当するパターン $\boldsymbol{x}$ のみである。

という前提のもと

クラスタの数はいくつか

パターン $\boldsymbol{x}$ はそれぞれどのクラスタに所属するか

各クラスタに想定した確率密度関数の最適なパラメータは何か

を推定するのがいわゆるクラスタリングです。方法として

混合正規分布のパラメータ推定

k-means法

凸クラスタリング法

等があります。続パタでは9章（混合正規分布のパラメータ推定）、10章（k-means法、凸クラスタリング法）で扱われています。ところが前者二つではクラスタ数は固定されており「1.クラスタ数はいくつか」は推定できず、最後の凸クラスタリング法ではパラメータが一定かつクラスタ間で同一という制約があるため「3.各クラスタに想定した確率密度関数の最適なパラメータは何か」が推定できません。そこで1,2,3のすべての推定が可能な方法としてノンパラメトリックベイズモデルを使用することになります。

ではこれら3つとノンパラベイズでは何が異なるのでしょうか。上に挙げたクラスタリング法はすべて最尤推定によるものでクラスタ構造についてモデル化されていません。一方ノンパラベイズではベイズ推定に基づきクラスタ構造をも事前分布としてモデル化するところが異なる部分です。

ここでは表題の通りディリクレ過程混合モデル（DPM: Dirichlet Process Mixture Model）というクラスタリング法について考察します。

具体的に何を求めればよいか

$n$ 個のパターンを観測した場合を考えます。 $\boldsymbol{x}$ をパターン集合、 $\boldsymbol{s}$ を各パターンがどのクラスタに所属しているかを示す潜在変数の集合として

$\begin{eqnarray} \boldsymbol{x} &=& \{x_1, \ldots , x_n\}\\ \boldsymbol{s} &=& \{s_1, \ldots , s_n\} \end{eqnarray}$

とあらわすこととし、また、 $c$ 種のクラスタを各々 $\omega_1, \ldots , \omega_c$ とし、この各クラスタの確率モデルのパラメータ集合 $\boldsymbol{\theta}$ を

$\boldsymbol{\theta}=\{\theta_1, \ldots , \theta_c\}$

とあらわすこととします。このときパターン $\boldsymbol{x}$ のみが観測されたという条件下で最適な所属クラスタ $\boldsymbol{s}$ と各クラスタのパラメータ $\boldsymbol{\theta}$ を求めればよいことになります。つまり

$\DeclareMathOperator*{\argmax}{arg\,max} %argmax$
$(\boldsymbol{\hat{s}}, \boldsymbol{\hat{\theta}})=\displaystyle\argmax_{\boldsymbol{s}, \boldsymbol{\theta}}{\{p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})\}}$

となる $(\boldsymbol{\hat{s}}, \boldsymbol{\hat{\theta}})$ を求めればよいわけです。これはテキストでは「クラスタリング法１」とされているものですが、まずはこの方法について考えていきます。

二つの問題

上式の $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ はベイズの定理より

$p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})\propto p(\boldsymbol{x}\vert\boldsymbol{s}, \boldsymbol{\theta})P(\boldsymbol{s})p(\boldsymbol{\theta})$

となります。式中で離散分布は $P(\cdot)$ 、連続分布は $p(\cdot)$ と分けて書いています。
ここで二つ解決するべき問題が出てきます。

事前確率 $P(\boldsymbol{s})$ をどのように定めるか

事後分布 $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ を最大にする $\boldsymbol{s}, \boldsymbol{\theta}$ をどのように求めるか

1.の事前確率 $P(\boldsymbol{s})$ は結果として

$P(\boldsymbol{s})=P(s_1, \ldots , s_n)=P_E(n_1, \ldots , n_c)=\displaystyle\frac{\alpha^c\displaystyle\prod_{i=1}^{c}{(n_i-1)!}}{\text{AF}(\alpha, n)}$

となります。ここで $n_i$ はクラスタ $\omega_i$ に所属する $x_k$ の数であり $\sum_{i=1}^{c}{n_i}=n$ です。また $\text{AF}(\cdot)$ は上昇階乗（Acsending Factorial）で

$\text{AF}(\alpha, n)\overset{\mathrm{def}}{=}\alpha\cdot (\alpha+1)\cdot\cdots\cdot (\alpha +n-1)$

と定義されます。また $\alpha$ は集中度パラメータであり、 $P_E(\cdot)$ はイーウェンスの抽出公式（Ewens's sampling formula）として知られています。何故事前確率 $P(\boldsymbol{s})$ をこのように定めるのが妥当なのか、また集中度パラメータ $\alpha$ の意味は何なのかについては話が長いので別途考察します。

では次に2.の事後分布 $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ を $\boldsymbol{s}, \boldsymbol{\theta}$ に関してどのように最大化すればよいかという問題です。何故これが問題になるかというと $\boldsymbol{s}$ の組み合わせが膨大な事から直接 $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ を求めることが困難なのです。具体的には $\boldsymbol{s}$ の組み合わせはベル数 $B_n$ になります。ベル数ではサンプル数が53、つまり $n=53$ 程度の組み合わせの数ですら $10^{51}$ 以上あります。ベル数についての具体的な話は別途行うとして、ここではギブスサンプリングを用いて解決します。

ギブスサンプリング

ギブスサンプリングは分布に従ってランダムなパターンを生成するアルゴリズムです。このアルゴリズムを用いて $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ に従うランダムな $(\boldsymbol{s}, \boldsymbol{\theta})$ を十分な数だけ生成すればその中には $(\boldsymbol{\hat{s}}, \boldsymbol{\hat{\theta}})$ に近いものも含まれるだろうという事を利用します。ギブスサンプリングでは一度に生成する元を一つに絞り、それ以外を固定した分布から順番に各元を生成します。ここでは $(\boldsymbol{s}, \boldsymbol{\theta})$ を生成するので、

$\boldsymbol{\theta}$ を固定し $p(\boldsymbol{s}\vert\boldsymbol{\theta}, \boldsymbol{x})$ に従う $\boldsymbol{s}$ を生成
$\boldsymbol{s}$ を固定し $p(\boldsymbol{\theta}\vert\boldsymbol{s}, \boldsymbol{x})$ に従う $\boldsymbol{\theta}$ を生成

を繰り返すというわけです。のちに具体的に述べますが2.の $p(\boldsymbol{\theta}\vert\boldsymbol{s}, \boldsymbol{x})$ だけならば解析的に解けるようにすることが可能です。ところが1.の $p(\boldsymbol{s}\vert\boldsymbol{\theta}, \boldsymbol{x})$ は $\boldsymbol{s}$ に関する分布であるため、まだ組み合わせ爆発を抱えたままです。そこでこの分布に対してもさらにギブスサンプリングを用います。

$\boldsymbol{s}$ の推定

以降の説明のために以下を定義します。

$\boldsymbol{x}_{-k}\overset{\mathrm{def}}{=}\{x_1, \ldots , x_{k-1}, x_{k+1}, \ldots , x_n\}$
$\boldsymbol{s}_{-k}\overset{\mathrm{def}}{=}\{s_1, \ldots , s_{k-1}, s_{k+1}, \ldots , s_n\}$

つまり $\boldsymbol{x}_{-k}$ は $\boldsymbol{x}$ から $x_k$ を除いた集合、 $\boldsymbol{s}_{-k}$ は $\boldsymbol{s}$ から $s_k$ を除いた集合です。この表記の下

現在のパターン集合 $\boldsymbol{x}$ から $x_k$ を抜き取る
$P(s_k\vert x_k, \boldsymbol{s}_{-k}, \boldsymbol{\theta})$ より確率的に新しい $s_k$ を生成
1.2.を各パターンすべてにわたり一度ずつ行う

という手続きがここでのギブスサンプリングの方法です。このようにして $n$ 個のパターンすべてにわたり一通り所属クラスタ $s_k$ を更新してできた新しい $\boldsymbol{s}=\{s_1, \ldots , s_n\}$ が $P(\boldsymbol{s}\vert\boldsymbol{x}, \boldsymbol{\theta})$ に従うランダムな $\boldsymbol{s}$ です。では $P(s_k\vert x_k, \boldsymbol{s}_{-k}, \boldsymbol{\theta})$ はどのようにして求めればよいでしょうか。具体的な計算は後に回して結果のみを提示すると

$\begin{eqnarray} P(s_k\vert x_k, \boldsymbol{s}_{-k}, \boldsymbol{\theta}) &\propto& \begin{cases} \displaystyle\frac{n_i'}{\alpha +n-1}\cdot p(x_k\vert\theta_i) &(s_k=\omega_i)\\ \displaystyle\frac{\alpha}{\alpha +n-1}\cdot\int{p(x_k\vert\theta_{\text{new}})p(\theta_{\text{new}})d\theta_{\text{new}}} &(s_k=\omega_{\text{new}}) \end{cases} \end{eqnarray}$

となります。式中の $\alpha$ は前述した集中度パラメータ、 $n_i'$ は $\boldsymbol{x}_{-k}$ のうちクラスタ $\omega_i$ に属するパターンの数を表します。この式は正規化項を無視して（正規化されているものかのごとく）述べれば、上段の計算結果の確率で各既存クラスタ $\omega_i$ に、下段の計算結果の確率で新規クラスタ $\omega_{\text{new}}$ に $s_k$ が割り当てられることを意味しています。実際には $P(s_k\vert x_k, \boldsymbol{s}_{-k}, \boldsymbol{\theta})$ に比例しているだけの式であり、正規化されていませんので実装時にはすべての計算を一通りしたのちに正規化します。

$\boldsymbol{\theta}$ の推定

続いて $\boldsymbol{\theta}$ の推定に必要な $p(\boldsymbol{\theta}\vert\boldsymbol{x}, \boldsymbol{s})$ の計算方法です。これは各クラスタに所属するパターンを用いてクラスタ毎に計算できます。すなわち $p(\theta_i\vert\{x_k\vert x_k\in\omega_i\})$ を各クラスタ毎に求めればよいことになります。 $p(\boldsymbol{\theta}\vert\boldsymbol{x}, \boldsymbol{s})$ 中の条件 $\boldsymbol{s}$ は条件 $\{x_k\vert x_k\in\omega_i\}$ に吸収されていることに注意します。つまり $\{x_k\vert x_k\in\omega_i\}$ の各 $x_k$ が $\omega_i$ に属しているかどうかは $\boldsymbol{s}$ によって決定されているということです。結果として $p(\boldsymbol{\theta}\vert\boldsymbol{x}, \boldsymbol{s})$ から $\boldsymbol{\theta}$ を推定することと $p(\theta_i\vert\{x_k\vert x_k\in\omega_i\})$ から各 $\theta_i$ を推定することは等価であることが分かります。言い換えれば各 $\theta_i$ の推定は $\omega_i$ に属していないパターンには影響されないために各 $\theta_i$ を個別に求められるということです。具体的にはベイズの定理より

$\begin{eqnarray} p(\theta_i\vert\{x_k\vert x_k\in\omega_i\})=\displaystyle\frac{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})}{\displaystyle\int{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})d\theta_i}} \end{eqnarray}$

となります。

クラスタの確率モデル

ここまでの式において、クラスタの確率モデル $p(x_k\vert\theta_i)$ と、 $\theta_i$ の事前分布 $p(\theta_i)$ を保留してきました。これら二つの分布の決定は同時に行う必要があります。というのもある確率モデルをクラスタの確率モデルとして適用する際、パラメータの事前分布 $p(\theta_i)$ をその確率モデルの共役事前分布とすることで上述の式のうち煩雑な二つの式が解析的に解けるからです。具体的には

$\displaystyle\int{p(x_k\vert\theta_{\text{new}})p(\theta_{\text{new}})d\theta_{\text{new}}}$
$\displaystyle\frac{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})}{\displaystyle\int{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})d\theta_i}}$

の二つの式です。ではテキストである続パタに従って各クラスタの確率モデルを平均 $\mu$ 、共分散行列 $\boldsymbol{\Sigma}$ が共に未知の $d$ 次元正規分布とし $\mathcal{N}(x;\mu, \boldsymbol{\Sigma})$ とあらわすとします。これはパターン全体の分布を混合正規分布であると仮定することと同じです。 $\theta_i=(\mu_i, \boldsymbol{\Sigma}_i)$ とあらわすことにすると

$p(x_k\vert\theta_i)=\mathcal{N}(x_k;\mu_i, \boldsymbol{\Sigma}_i)$

となります。この平均 $\mu$ 、共分散行列 $\boldsymbol{\Sigma}$ が共に未知の $d$ 次元正規分布の場合、 $p(\theta_i)$ に設定するべき共役事前分布は下記の正規ウィシャート分布となります。

$p(\theta_i)=p(\mu_i, \boldsymbol{\Lambda}_i)=\mathcal{N}(\mu_i;\mu_0, (\beta\boldsymbol{\Lambda}_i)^{-1})\mathcal{W}(\boldsymbol{\Lambda}_i;\nu, \mathbf{S})$

ここで $\mu_0, \beta, \nu, \mathbf{S}$ はすべてハイパーパラメータです。 $\nu$ は自由度と呼ばれるスカラー、 $\mathbf{S}$ は $d\times d$ 正定値対象行列です。また $\boldsymbol{\Lambda}$ は精度行列と呼ばれ共分散行列 $\boldsymbol{\Sigma}$ の逆行列つまり $\boldsymbol{\Lambda}=\boldsymbol{\Sigma}^{-1}$ です。以降は計算の便宜上、主に共分散行列ではなく精度行列を用います。確率モデルが定まったので上述の二つの式を具体的に解くことが出来ます。やはり具体的な計算は読書会10と読書会11, 12に回して結果のみを提示すると

$\displaystyle\int{p(x_k\vert\theta_{\text{new}})p(\theta_{\text{new}})d\theta_{\text{new}}}=\left(\displaystyle\frac{\beta}{(1+\beta)\pi}\right)^{\frac{d}{2}}\cdot\displaystyle\frac{\vert \mathbf{S}_b\vert^{\frac{\nu+1}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1}{2}\right)}{\vert\mathbf{S}\vert^{\frac{\nu}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1-d}{2}\right)}$

ただし $\mathbf{S}_b$ は以下

$\mathbf{S}_b^{-1}=\mathbf{S}^{-1}+\displaystyle\frac{\beta}{1+\beta}(x_k-\mu_0)(x_k-\mu_0)^t$

となります。また、

$\displaystyle\frac{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})}{\displaystyle\int{p(\theta_i)\cdot\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i})d\theta_i}}=\mathcal{N}(\mu_i;\mu_c, \boldsymbol{\Lambda}_c^{-1})\mathcal{W}(\boldsymbol{\Lambda}_i;\nu_c, \mathbf{S}_q)$

ただし $\mu_c, \boldsymbol{\Lambda}_c, \nu_c, \mathbf{S}_q$ は以下

$\begin{eqnarray} \bar{x} &=& \displaystyle\frac{1}{n_i}\displaystyle\sum_{x_k\in\omega_i}{x_k}\\ \mu_c &=& \displaystyle\frac{n_i\bar{x}+\beta\mu_0}{n_i+\beta}\\ \boldsymbol{\Lambda}_c &=& (n_i+\beta)\boldsymbol{\Lambda}_i\\ \nu_c &=& \nu+n_i\\ \mathbf{S}_q^{-1} &=& \mathbf{S}^{-1}+\displaystyle\sum_{x_k\in\omega_i}{(x_k-\bar{x})(x_k-\bar{x})^t}+\displaystyle\frac{n_i\beta}{n_i+\beta}(\bar{x}-\mu_0)(\bar{x}-\mu_0)^t \end{eqnarray}$

となります。

今回はここまでにして次回は具体的なアルゴリズムを記述していきます。