ディリクレ分布から中華料理店過程の分割ルールを導出する

公開:2020/11/08
更新:2020/11/23
8 min読了の目安(約7800字TECH技術記事
\gdef\ds{\displaystyle} \gdef\b#1{\mathbf{#1}} \gdef\bs#1{\boldsymbol{#1}}

まえがき

ディリクレ分布から中華料理店仮定(CRP)の分割ルールを導出します。CRPの詳細は参考にさせて頂いた『続・わかりやすい パターン認識 -教師なし学習入門-』を参照してください。

前準備

  • サンプルの数をnnとします。
  • クラスタの数をcc(ただしn<<cn << c)とし、各クラスタをωj  (j=1,,c)\omega_j \; (j = 1, \ldots, c)で表現します。
  • 各サンプルへのクラスタの割り当てをsks_kとして表現します。sks_kは1つのクラスタのみに属するとします。
  • sks_kωi\omega_iが割り当てられる確率をπi\pi_iとします。j=1cπj=1\sum_{j = 1}^c \pi_j = 1です。
  • π=[π1,,πc],α=[α1,,αc]\bs{\pi} = \left[\pi_1, \ldots, \pi_c \right]^\top, \bs{\alpha} = \left[\alpha_1, \ldots, \alpha_c \right]^\topという表記を用います。
  • π\bs{\pi}の確率分布として、以下のようにディリクレ分布を仮定します。
p(π)=Dir(πα)=Γ(j=1cαj)j=1cΓ(αj)j=1cπjαj1 p(\bs{\pi}) = \mathrm{Dir} ( \bs{\pi} | \bs{\alpha} ) = \frac{\Gamma\left(\sum_{j = 1}^c \alpha_j \right)}{\prod_{j = 1}^c \Gamma (\alpha_j)} \prod_{j = 1}^c \pi_j^{\alpha_j - 1}
  • π\bs{\pi}が与えられたとき、sks_kはカテゴリ分布に従うと考えられます。
p(sk=ωiπ)=Cat(sk=ωiπ)=πi p(s_k = \omega_i | \bs{\pi}) = \mathrm{Cat} (s_k = \omega_i | \bs{\pi}) = \pi_i

導出

事前に、今回想定しているグラフィカルモデルを示します。α\bs{\alpha}からπ\bs{\pi}が確率的に生成され、π\bs{\pi}からs1,,sns_1, \ldots, s_nが生成されるという構造になっています。

graphical-model
想定するグラフィカルモデル

CRPの分割ルールを導出するため、s1,,sn1s_1, \ldots, s_{n - 1}が与えられたとき、sns_nがあるクラスタωi\omega_iに割り当てられる確率p(sn=ωis1,,sn1)p(s_n = \omega_i | s_1, \ldots, s_{n - 1})を求めます。α\bs{\alpha}は与えられていますが、π\bs{\pi}は未知です。ここで、1つのパラメータα\alphaを用いて各クラスタのαj=α/c\alpha_j = \alpha / cとして、全て等しいとします。このとき、p(π)p(\bs{\pi})

p(π)=Dir(πα)=Γ(α)Γ(α/c)cj=1cπjα/c1 p(\bs{\pi}) = \mathrm{Dir} ( \bs{\pi} | \alpha ) = \frac{\Gamma(\alpha)}{\Gamma (\alpha / c)^c} \prod_{j = 1}^c \pi_j^{\alpha / c - 1}

となっています。s1,,sn1s_1, \ldots, s_{n - 1}の中でωj\omega_jに割り当てられた個数をnjn_jとすると、

p(s1,,sn1π)=k=1n1p(skπ)=j=1cπjnj p(s_1, \ldots, s_{n - 1} | \bs{\pi}) = \prod_{k = 1}^{n - 1} p(s_k | \bs{\pi}) = \prod_{j = 1}^c \pi_j^{n_j}

です。これを用いてp(sn=ωis1,,sn1)p(s_n = \omega_i | s_1, \ldots, s_{n - 1})を求めます。

p(sn=ωis1,,sn1)=p(sn=ωi,s1,,sn1)p(s1,,sn1) \begin{aligned} p(s_n = \omega_i | s_1, \ldots, s_{n - 1}) = \frac{p(s_n = \omega_i, s_1, \ldots, s_{n - 1})}{p(s_1, \ldots, s_{n - 1})} \end{aligned}

計算が煩雑になるため、まずは分母p(s1,,sn1)p(s_1, \ldots, s_{n - 1})を求めます。π\bs{\pi}についての周辺化を考慮することで

p(s1,,sn1)=p(s1,,sn1,π)dπ=p(s1,,sn1π)p(π)dπ={j=1cπjnj}Γ(α)Γ(α/c)cj=1cπjα/c1dπ=Γ(α)Γ(α/c)cj=1cπjnj+α/c1dπ=Γ(α)Γ(α/c)cj=1cΓ(nj+α/c)Γ(j=1c(nj+α/c))(補足を参照してください)=Γ(α)Γ(α/c)cj=1cΓ(nj+α/c)Γ(n1+α)(j=1cnj=n1,j=1cα/c=α)=Γ(α)Γ(n1+α)j=1cΓ(nj+α/c)Γ(α/c) \begin{aligned} & p(s_1, \ldots, s_{n - 1}) \\ &= \int p(s_1, \ldots, s_{n - 1}, \bs{\pi}) d \bs{\pi} = \int p(s_1, \ldots, s_{n - 1} | \bs{\pi}) p(\bs{\pi}) d \bs{\pi} \\ &= \int \left\{ \prod_{j = 1}^c \pi_j^{n_j} \right\} \frac{\Gamma(\alpha)}{\Gamma(\alpha / c)^c} \prod_{j = 1}^c \pi_j^{\alpha / c - 1} d \bs{\pi} = \int \frac{\Gamma(\alpha)}{\Gamma(\alpha / c)^c} \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1} d \bs{\pi} \\ &= \frac{\Gamma(\alpha)}{\Gamma(\alpha / c)^c} \frac{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)}{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))} \quad \left(\because \mathrm{補足を参照してください} \right) \\ &= \frac{\Gamma(\alpha)}{\Gamma(\alpha / c)^c} \frac{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)}{\Gamma(n - 1 + \alpha)} \quad \left(\because \sum_{j = 1}^c n_j = n - 1, \sum_{j = 1}^c \alpha / c = \alpha \right) \\ &= \frac{\Gamma(\alpha)}{\Gamma(n - 1 + \alpha)} \prod_{j = 1}^c \frac{\Gamma(n_j + \alpha / c)}{\Gamma(\alpha / c)} \end{aligned}

が得られます。分子p(sn=ωi,s1,,sn1)p(s_n = \omega_i, s_1, \ldots, s_{n - 1})ではs1,,sn1s_1, \ldots, s_{n - 1}に加えてsn=ωis_n = \omega_iとなる確率を求めることになります。これは分母p(s1,,sn1)p(s_1, \ldots, s_{n - 1})を求めるときと同様の手順によって式変形を行うことにより、

p(sn=ωi,s1,,sn1)=p(sn=ωi,s1,,sn1,π)dπ=p(sn=ωi,s1,,sn1π)p(π)dπ=p(sn=ωiπ){k=1n1p(skπ)}p(π)dπ=Γ(α)Γ(n+α)Γ(ni+1+α/c)Γ(α/c){j=1jiΓ(nj+α/c)Γ(α/c)} \begin{aligned} & p(s_n = \omega_i, s_1, \ldots, s_{n - 1}) \\ &= \int p(s_n = \omega_i, s_1, \ldots, s_{n - 1}, \bs{\pi}) d \bs{\pi} = \int p(s_n = \omega_i, s_1, \ldots, s_{n - 1} | \bs{\pi}) p(\bs{\pi}) d \bs{\pi} \\ &= \int p(s_n = \omega_i | \bs{\pi}) \left\{ \prod_{k = 1}^{n - 1} p(s_k | \bs{\pi}) \right\} p(\bs{\pi}) d \bs{\pi} \\ &= \frac{\Gamma(\alpha)}{\Gamma(n + \alpha)} \frac{\Gamma(n_i + 1 + \alpha / c)}{\Gamma(\alpha / c)} \left\{ \prod_{\substack{j = 1 \\ j \neq i}} \frac{\Gamma(n_j + \alpha / c)}{\Gamma(\alpha / c)} \right\} \end{aligned}

となります。したがって、求めたいp(sn=ωis1,,sn1)p(s_n = \omega_i | s_1, \ldots, s_{n - 1})

p(sn=ωis1,,sn1)=p(sn=ωi,s1,,sn1)p(s1,,sn1)=Γ(α)Γ(n+α)Γ(ni+1+α/c)Γ(α/c){j=1jiΓ(nj+α/c)Γ(α/c)}Γ(n1+α)Γ(α)j=1cΓ(α/c)Γ(nj+α/c)=Γ(α)Γ(n+α)Γ(ni+1+α/c)Γ(α/c){j=1jiΓ(nj+α/c)Γ(α/c)}Γ(n1+α)Γ(α)j=1cΓ(α/c)Γ(nj+α/c)=Γ(n1+α)Γ(n+α)Γ(ni+1+α/c)Γ(ni+α/c)=1n1+α(ni+α/c)=ni+α/cn1+α \begin{aligned} & p(s_n = \omega_i | s_1, \ldots, s_{n - 1}) \\ &= \frac{p(s_n = \omega_i, s_1, \ldots, s_{n - 1})}{p(s_1, \ldots, s_{n - 1})} \\ &= \frac{\Gamma(\alpha)}{\Gamma(n + \alpha)} \frac{\Gamma(n_i + 1 + \alpha / c)}{\Gamma(\alpha / c)} \left\{ \prod_{\substack{j = 1 \\ j \neq i}} \frac{\Gamma(n_j + \alpha / c)}{\Gamma(\alpha / c)} \right\} \cdot \frac{\Gamma(n - 1 + \alpha)}{\Gamma(\alpha)} \prod_{j = 1}^c \frac{\Gamma(\alpha / c)}{\Gamma(n_j + \alpha / c)} \\ &= \frac{\bcancel{\Gamma(\alpha)}}{\Gamma(n + \alpha)} \frac{\Gamma(n_i + 1 + \alpha / c)}{\bcancel{\Gamma(\alpha / c)}} \left\{ \prod_{\substack{j = 1 \\ j \neq i}} \frac{\Gamma(n_j + \alpha / c)}{\bcancel{\Gamma(\alpha / c)}} \right\} \cdot \frac{\Gamma(n - 1 + \alpha)}{\bcancel{\Gamma(\alpha)}} \prod_{j = 1}^c \frac{\bcancel{\Gamma(\alpha / c)}}{\Gamma(n_j + \alpha / c)} \\ &= \frac{\Gamma(n - 1 + \alpha)}{\Gamma(n + \alpha)} \frac{\Gamma(n_i + 1 + \alpha / c)}{\Gamma(n_i + \alpha / c)} \\ &= \frac{1}{n - 1 + \alpha} \cdot (n_i + \alpha / c) = \frac{n_i + \alpha / c}{n - 1 + \alpha} \end{aligned}

と計算することができます。

cc個の全てのクラスタを含む集合をΩ\Omegaとします。ここで、n<<cn << cと仮定していることから、ni=0n_i = 0、即ち1つも属するサンプルが存在しないクラスタが存在します。Ω\Omegaを1つも属するサンプルが存在しないクラスタの集合Ω0={ωini=0,i=1,,c}\Omega_0 = \left\{\omega_i | n_i = 0, i = 1, \ldots, c \right\}およびそれ以外のクラスタの集合Ω1={ωini0,i=1,,c}\Omega_1 = \left\{\omega_i | n_i \neq 0, i = 1, \ldots, c \right\}に分割します。つまりΩ=Ω0Ω1,=Ω0Ω1,c=Ω=Ω1+Ω2\Omega = \Omega_0 \cup \Omega_1, \emptyset = \Omega_0 \cap \Omega_1, c = |\Omega| = |\Omega_1| + |\Omega_2|です。
このとき、

p(snΩ0s1,,sn1)=ωiΩ0p(sn=ωis1,,sn1)=ωiΩ0ni+α/cn1+α=ωiΩ0α/cn1+α=cΩ1cαn1+α \begin{aligned} & p(s_n \in \Omega_0 | s_1, \ldots, s_{n - 1}) = \sum_{\omega_i \in \Omega_0} p(s_n = \omega_i | s_1, \ldots, s_{n - 1}) = \sum_{\omega_i \in \Omega_0} \frac{n_i + \alpha / c}{n - 1 + \alpha} \\ &= \sum_{\omega_i \in \Omega_0} \frac{\alpha / c}{n - 1 + \alpha} = \frac{c - |\Omega_1|}{c} \frac{\alpha}{n - 1 + \alpha} \end{aligned}

であり、cc \rightarrow \inftyの極限では

p(snΩ0s1,,sn1)=(1Ω1c)αn1+ααn1+α(c) \begin{aligned} & p(s_n \in \Omega_0 | s_1, \ldots, s_{n - 1}) = \left(1 - \frac{|\Omega_1|}{c} \right) \frac{\alpha}{n - 1 + \alpha} \rightarrow \frac{\alpha}{n - 1 + \alpha} \quad (c \rightarrow \infty) \end{aligned}

です。一方、Ω1\Omega_1に属するクラスタに対しては

p(sn=ωi Ω1s1,,sn1)=ni+α/cn1+αnin1+α(c) p(s_n = \omega_i \in \ \Omega_1 | s_1, \ldots, s_{n - 1}) = \frac{n_i + \alpha / c}{n - 1 + \alpha} \rightarrow \frac{n_i}{n - 1 + \alpha} \quad (c \rightarrow \infty)

です。以上の計算をまとめると

{p(snΩ0s1,,sn1)=αn1+αp(sn=ωi Ω1s1,,sn1)=nin1+α \begin{cases} \ds p(s_n \in \Omega_0 | s_1, \ldots, s_{n - 1}) = \frac{\alpha}{n - 1 + \alpha} \\ \\ \ds p(s_n = \omega_i \in \ \Omega_1 | s_1, \ldots, s_{n - 1}) = \frac{n_i}{n - 1 + \alpha} \end{cases}

となります。sns_nni0n_i \neq 0の既存のクラスタωiΩ1\omega_i \in \Omega_1のいずれかに属する確率とsns_nni=0n_i = 0の新規のクラスタωiΩ0\omega_i \in \Omega_0に属する確率はni:αn_i : \alphaの比になっています。したがって、これはCRPの分割ルールに対応していることが分かります。

補足

j=1cπjnj+α/c1dπ=j=1cΓ(nj+α/c)Γ(j=1c(nj+α/c)) \int \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1} d \bs{\pi} = \frac{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)}{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))}

を示します。
α~=[n1+α/c,,nc+α/c]\bs{\tilde{\alpha}} = \left[n_1 + \alpha / c, \ldots, n_c + \alpha / c \right]^\topとすると、α~\bs{\tilde{\alpha}}をパラメータとするディリクレ分布は

Dir(πα~)=Γ(j=1c(nj+α/c))j=1cΓ(nj+α/c)j=1cπjnj+α/c1 \mathrm{Dir} (\bs{\pi} | \bs{\tilde{\alpha}}) = \frac{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))}{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)} \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1}

と表せます。ディリクレ分布は確率分布なので、π\bs{\pi}で積分すると1になります。

Dir(πα~)dπ=Γ(j=1c(nj+α/c))j=1cΓ(nj+α/c)j=1cπjnj+α/c1dπ=Γ(j=1c(nj+α/c))j=1cΓ(nj+α/c)j=1cπjnj+α/c1dπ=1 \begin{aligned} \int \mathrm{Dir} (\bs{\pi} | \bs{\tilde{\alpha}}) d \bs{\pi} &= \int \frac{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))}{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)} \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1} d \bs{\pi} \\ &= \frac{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))}{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)} \int \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1} d \bs{\pi} \\ &= 1 \end{aligned}

分数の中にはπ\bs{\pi}を含んでいないため、積分の対象にはなっていません。よって、この結果から

j=1cπjnj+α/c1dπ=j=1cΓ(nj+α/c)Γ(j=1c(nj+α/c)) \int \prod_{j = 1}^c \pi_{j}^{n_j + \alpha / c - 1} d \bs{\pi} = \frac{\prod_{j = 1}^c \Gamma(n_j + \alpha / c)}{\Gamma(\sum_{j = 1}^c (n_j + \alpha / c))}

が得られます。

参考文献