機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan:離散パラメータの扱い

隠れ変数の扱い


前回、ポアソン混合モデルを勉強した。

kento1109.hatenablog.com

前回は触れなかったが、ポアソン混合モデルのグラフィカルモデルは一般的に書きのように書ける。

f:id:kento1109:20180428102058p:plain

観測データx_nは、k番目クラスタのパラメータ\lambda_kポアソン分布によって生成される。
どのクラスタが選択されるかは、隠れ変数s_nによって決定される。
この隠れ変数s_nはカテゴリ分布から生成される1 of kベクトルである。つまり、k番目の次元が「1」で、それ以外が「0」を取るようなベクトルである。
つまり、s_nは離散変数である。しかし、Stanでは離散変数を扱えないので、s_nをパラメータに定義することが出来ない。
Stanで実装する場合、このような隠れ変数は「周辺化除去(marginalization )」することで対応する。

marginalization



周辺化除去の簡単な例を用いて確認する。
離散変数の場合、除去したい変数を総和することでその変数を除去する。

以下は箱(A or B)の中にある玉の色(赤 or 白)から玉を取り出した時の確率を表したもの。
f:id:kento1109:20180428105110p:plain:w300
絵にするとこんな感じ。
f:id:kento1109:20180428105615p:plain

  • 同時確率

「箱がA&玉が赤色」の確率は以下のように書ける。
 \begin{eqnarray}
p(x=A,y=R)=\frac{3}{10}
 \end{eqnarray}

  • 周辺確率

「(箱に関係なく)玉が赤色」の確率は以下のように書ける。
 \begin{eqnarray}
p(y=R)=\sum_{x'}{p(x',y)}=\frac{3}{10}+\frac{1}{10}=\frac{4}{10}
 \end{eqnarray}

  • 条件付き確率

「玉が赤色だった時、それが箱Aから取り出された」の確率は以下のように書ける。
 \begin{eqnarray}
p(x=A|y=R)=\frac{p(x,y)}{\sum_{x'}{p(x',y)}}=\frac{\frac{3}{10}}{\frac{3}{10}+\frac{1}{10}}=\frac{3}{4}
 \end{eqnarray}

周辺確率でやったことは、「求めたい変数以外の変数の全ての取り得る確率の足し合わせ」である。
例えば、赤玉の個数は、「箱Aに入っている赤玉+箱Bに入っている赤玉」である。
箱がn個あった場合も同じであり、「箱1の赤玉+箱2の赤玉+...+箱nの赤玉」の個数となる。
確率も同じように全て足し込むことで計算できる。

では、ポアソン混合モデルに戻って考える。
観測データx_nの生成確率は、
\begin{eqnarray}
p(x_n|s_n,\lambda)=\prod_{k=1}^K{{\rm Poi}(x_n|\lambda_k)^{s_{n,k}}}
 \end{eqnarray}
で表せた。
今回は隠れ変数s_nが取り得る値を全て足し込むことで、変数s_nを除去する。
隠れ変数s_nは「n行k列の行列」であるため、この2つの変数に対して足し込めばよい。
今回、n=1(s_1)の場合のみ確認する。
s_1が取り得る値はk個存在する。
例えば、k=1の場合、観測データx_1の生成確率は、
\begin{eqnarray}
p(x_1|\lambda_1, \pi_1)={\rm Poi}(x_1|\lambda_1)\times p(\pi_1)
 \end{eqnarray}
と表すことが出来る。
k=2の場合、観測データx_1の生成確率は、
\begin{eqnarray}
p(x_1|\lambda_2, \pi_2)={\rm Poi}(x_1|\lambda_2)\times p(\pi_2)
 \end{eqnarray}
と表すことが出来る。
これをk個分足し込むので、観測データx_1の生成確率は、
\begin{eqnarray}
p(x_1|\lambda_k, \pi_k)=\sum_{k=1}^K{{\rm Poi}(x_1|\lambda_k)\times p(\pi_k)}
 \end{eqnarray}
で書くことが出来る。
これが観測データx_1の生成確率である。
独立した観測データx_nの生成確率は全てを掛けることで求められる。
\begin{eqnarray}
p(x_n|\lambda_k, \pi_k)=\prod_{n=1}^N\left\{{\sum_{k=1}^K{{\rm Poi}(x_n|\lambda_k)\times p(\pi_k)}}\right\}
 \end{eqnarray}
また、対数尤度は以下のとおりである。
\begin{eqnarray}
\log p(x_n|\lambda_k, \pi_k)&=&\log \prod_{n=1}^N\left\{{\sum_{k=1}^K{{\rm Poi}(x_n|\lambda_k)\times p(\pi_k)}}\right\}\\
\log p(x_n|\lambda_k, \pi_k)&=&\sum_{n=1}^N \log \left\{{\sum_{k=1}^K{{\rm Poi}(x_n|\lambda_k)\times p(\pi_k)}}\right\}
 \end{eqnarray}

log_sum_exp



最後に和の対数を計算したい場合に用いるテクニックである「log_sum_exp」を利用する。
log_sum_expを用いると和の対数を以下のように変換できる。
\log \sum_{i=1}^N {(x_i)}={\rm log\_sum\_exp}(\log x_i)
今回のモデルで「log_sum_exp」を用いると最終的にモデル式は以下のように書ける。
\begin{eqnarray}
\log p(x_n|\lambda_k, \pi_k)=\sum_{n=1}^N {\rm log\_sum\_exp} \left\{{\rm \log Poi}(x_n|\lambda_k)\times \log p(\pi_k)\right\}
 \end{eqnarray}

この処理をStanで書くと以下のように書ける。

model {
  real lp[K];
  for(n in 1:N){
    for (k in 1:K){
      lp[k] = log(pi[k]) + poisson_lpmf(x[n] | lambda[k]);
    }
    target += log_sum_exp(lp);
  }
}

これで、隠れ変数を直接扱わずに済んだ。
「周辺化ギブスサンプリング」も特定の変数を周辺化して除去することで高速化を実現している。
あれは連続変数でも問題無かった(解析的に積分計算ができる必要がある)が、考え方は似ているのかもしれないと思った。

※今回の内容は下記を参考とさせて頂いた。
statmodeling.hatenablog.com