xiangze's sparse blog

機械学習、ベイズ統計、コンピュータビジョンと関連する数学について

Stan2.0では離散変数をparameterとできないことについて

Stan2.0で配列のindexとしてCategorical分布から生成したものは使えないためにモデルの記述がわかりにくくなってしまう場合があります。

Stanマニュアルの説明が簡潔すぎて自分にはすぐに理解できなかったためその補足です。

Stanで混合モデルを記述する場合にこの問題に直面します。

該当部分はberoberoさんが日本語訳されています(11. Mixture Modeling)。

混合モデルは一般的には観測変数y_i,その分布関数をA,  パラメータをm, indexの確率変数z_i、そのパラメータをθとして

 z_i \sim Categorical(\theta)

 y_i \sim A(m_{z_i})

と書くことができます*1

全体の確率密度関数はカテゴリの総数をKとすると各カテゴリkになる確率とその場合のAの値の積の総和になるので

 p(y_i)=\sum_{k=0}^{K} \theta_k A(m_{z_k})

と書くことができます。z_iの分布が式中に現れなくなったので、Stanの制約があってもモデルが記述できるということになります。

 

あるサンプリングした変数の集合に対するモデルの尤度はp(y_i)の積

 \prod_{i=0}^{N} p(y_i) =\prod_{i=0}^{N} \sum_{k=0}^{K} \theta_k A(m_{z_k} )

になりますが、計算上p(y_i) (=p_i)の値が非常に大きく、あるいは小さくなる可能性があるのでそのlogを取ったものを用いるのが常道らしいです。その場合log(p_i)は

 log(p_i)=\sum_{i=0}^{N} log( \sum_{k=0}^{K} \theta_k A(m_{z_k}))

 =\sum_{i=0}^{N}log( \sum_{k=0}^{K} exp ( log(\theta_k)+logA(m_{z_k})))

 =\sum_{i=0}^{N}logsumexp( log(\theta_k)+logA(m_{z_k}))

 

とかけます。log,sum,expの並びはlogsumexpという一つの関数で書いています*2

和の部分はstanではincremental_log_probという関数を使うことができ、最終的にStanマニュアルのpage 85のように

 

real ps[K];

for (n in 1:N) {

  for (k in 1:K) {

    ps[k] <- log(theta[k])+log(A(m_{z_k}));

    }

   increment_log_prob(log_sum_exp(ps));

}

 

とかけます*3

 LDA

 Latent Dirichlet Allocation (LDA)はグループ(document)毎にグループ分けされた観測変数(word)wがそれぞれ隠れ変数(トピック)zに対応した分布φzに従っているようなモデルです。

事前分布としてパラメータα、βのディリクレ分布

Dir(\vec{x}|\vec{\alpha})=\frac{1}{B(\vec{\alpha})}\prod_{i=0}^N x_i^{\alpha_i-1}

を用いた

\vec{\theta_m } \sim Dir(\alpha)

\vec{\phi_k } \sim Dir(\beta)

 z_ {m,n}\sim Categorical(\vec{\theta_m})

 w_{m,n} \sim Categorical(\vec{\phi_{z_{m,n}}})

が2重に出てくるモデルで、文書の分類などに用いられるそうです。

f:id:xiangze:20131219222903p:plain

φのindexに単語ごとに定義されたトピック変数zが使われていますが、これのCategorical分布をあらわに含まないように書き換えようとすると、以下のようになります。

Kをトピックの数として確率密度関数は変数間の依存関係にしたがって

 p(\theta,\phi|w,\alpha,\beta)\propto p(\theta|\alpha) p(\phi|\beta) p(w|\theta,\phi)

 =\prod_m^M p(\theta_m|\alpha) \prod_k^K p(\phi_k,\beta) \prod_m^M\prod_n^{M_n}p(w_{m,n}|\theta_m,\phi)

となります。wの分布関数は隠れ変数zを用いてさらに

 p(w_{m,n}|\theta ,\phi)=\sum_{z=1}^{K} p(z,w_{m,n} |\theta_m \phi)

 =\sum_{z=1}^{K} p(z,\theta_m) p(w_{m,n} | \phi_z)

と分解できます。これのlogをとると

 log(p(\theta,\phi|w,\alpha,\beta)) =\sum_m^M log(p(\theta_m|\alpha))+\sum_k^K log(p(\phi_k|\beta)+\sum_m^M\sum_n^{N_m} log(p(w_{m,n}|\phi_z)

 =\sum_m^M Dir(\theta_m|\alpha)+\sum_k^K Dir(\phi_k|\beta)+\sum_m^M\sum_n^{N_m}log(\sum_{z=1}^{K}Categorical(z|\theta_m)Categorical(w_{m,n}|\phi_z))

となります。和の部分 \sum_m^M\sum_n^{N_m}はドキュメントm

とそれに属する単語のindex n(<=N_M)に対するものです。

この第3項ではトピックzに関する総和がとられているのでCategorical分布は現れなくなり、文書mごとに存在するベクトルθmのトピックzに対応する成分θm[z],またword nごとのベクトルφzの n[w]成分が選択されることになります。

 \sum_m^M\sum_n^{N_m}log(\sum_{z=1}^{K}\theta_{m,z}\phi_{z,w_{m,n}})

 \sum_m^M\sum_n^{N_m}log(\sum_{z=1}^{K} exp (log(\theta_{m,z})+log(\phi_{z,w_{m,n}})))

となります。Stanでは要素数が変動するarray(ragged array)をサポートしていないこともあってwの要素ごとにドキュメントのindex mを与えてるようにしています。これをdoc[n]とすると普通の混合モデルの場合と同様にpage 128のようなStanのコード

for (n in 1:N) {

 real gamma[K];

  for (k in 1:K) 

    gamma[k]<-log(theta[doc[n],k])+log(phi[k,w[n]]);

  increment_log_prob(log_sum_exp(gamma));

  }

になります。

 

Reference

Stanのマニュアルの8章~12章の私的メモ

 

おまけ: HaskellでのMCMC,学習の実装

increment_log_probなどincrementalな記述はモデルの数式との差異が大きく、コードに落とし込みづらい、あるいStanのC++に由来した厳密な型の取り扱いをうまく使いたい、型推論を使いたいと感じる方にはHaskellベースでのHierarchical Bayes Compilerを用いるのがよいかもしれません。LDAも非常に簡単に記述することができます。

またパラメータ推定を学習ととらえるのであれば逐次的な学習とdatasetの追加が可換であることを利用して効率的に学習できるらしいHlearnというライブラリがあるそうです。

参考:

The categorical distribution’s algebraic structure  

Algebraic classiers: a generic approach to fast cross-validation, online training, and parallel training(pdf)

 あとで読みます...

*1:マニュアルではCategorical分布としており、他の専門書でもMultinodal distribution(多項分布)としている場合がありますが、 後者は前者の特別な場合であり、厳密には異なります

*2:scipyなどにもある関数です。参考:logsumexp (log sum exponential)とは

*3:マニュアルでは分布Aは正規分布(normal),log(A)はnormal_logです