Stanï¼šéš ã‚Œãƒžãƒ«ã‚³ãƒ•ãƒ¢ãƒ‡ãƒ«â‘

はじめに

---

今回は隠れマルコフモデルをStanで実装する。
隠れマルコフモデル自体は以前に書いた。
kento1109.hatenablog.com

今回は教師ありモデルを考える。
教師あり「隠れ状態」が既知のモデル。
前回の例で考えると、「晴れ→雨」などの遷移状態が与えられているモデル。
その状態から、「遷移行列」と「出力行列」を推定する。

モデル

---

コードは以下のようにシンプル。

data {
  int<lower=1> K;  // カテゴリーの数
  int<lower=1> V;  // 単語(word)の数
  int<lower=0> T;  // 時点の数
  int<lower=1,upper=V> w[T]; // 単語(word)
  int<lower=1,upper=K> z[T]; // カテゴリー
  vector<lower=0>[K] alpha;  // 推移(transit)確率の事前確率
  vector<lower=0>[V] beta;   // 単語vを出力(emit)する確率の事前確率
}
parameters {
  simplex[K] theta[K];  // 推移(transit)確率
  simplex[V] phi[K];    // 単語vを出力(emit)する確率
}
model {
  for (k in 1:K)
    theta[k] ~ dirichlet(alpha);
  for (k in 1:K)
    phi[k] ~ dirichlet(beta);
  for (t in 1:T)
    w[t] ~ categorical(phi[z[t]]);
  for (t in 2:T)
    z[t] ~ categorical(theta[z[t - 1]]);
}

これは、
Stan モデリング言語: ユーザーガイド・リファレンスマニュアル
の「10.6. 隠れマルコフモデル」のまんま。

が遷移行列、が出力行列に該当する。

実行

---

データは以下のように作成した。

emission <- c(1,2,3,3,3,2,2,1,2)
hidden_state <- c(1,1,2,2,2,2,1,1,1)

stan.dat <- list(K=2, V=3, T=9, w=emission, z=hidden_state,
                 alpha=c(1.0,1.0), beta=c(1.0,1.0,1.0))

データを前回の例を用いて表にするとこんな感じ。

day	emission	hidden_state
1	walk	Sunny
2	shop	Sunny
3	clean	Rainy
4	clean	Rainy
5	clean	Rainy
6	shop	Rainy
7	shop	Sunny
8	walk	Sunny
9	shop	Sunny

これをに渡してあげる。

stan.fit <- stan("supervised_hmm.stan", data=stan.dat)

結果

---

以下の通り。収束判定は問題なさそう。

summary(stan.fit)
$summary
                  mean     se_mean        sd          2.5%          25%
theta[1,1]   0.6687044 0.002852309 0.1803958   0.285459555   0.54520998
theta[1,2]   0.3312956 0.002852309 0.1803958   0.050861661   0.18992343
theta[2,1]   0.3390158 0.002919518 0.1846466   0.052310405   0.19037473
theta[2,2]   0.6609842 0.002919518 0.1846466   0.272114157   0.53280306
phi[1,1]     0.3739677 0.002550008 0.1612767   0.102355908   0.24580265
phi[1,2]     0.4998130 0.002675016 0.1691828   0.179900530   0.37727482
phi[1,3]     0.1262193 0.001702130 0.1076522   0.003180911   0.04258686
phi[2,1]     0.1434164 0.001966619 0.1243799   0.004601771   0.04711737
phi[2,2]     0.2863562 0.002532753 0.1601854   0.044368247   0.15728610
phi[2,3]     0.5702274 0.002743551 0.1735174   0.223375958   0.44862691
                    50%         75%       97.5%    n_eff      Rhat
theta[1,1]   0.68911775   0.8100766   0.9491383 4000.000 0.9997148
theta[1,2]   0.31088225   0.4547900   0.7145404 4000.000 0.9997148
theta[2,1]   0.31960050   0.4671969   0.7278858 4000.000 0.9999138
theta[2,2]   0.68039950   0.8096253   0.9476896 4000.000 0.9999138
phi[1,1]     0.36392534   0.4872123   0.7034251 4000.000 0.9997530
phi[1,2]     0.50182133   0.6216414   0.8188262 4000.000 1.0001396
phi[1,3]     0.09887729   0.1803597   0.3941908 4000.000 0.9998056
phi[2,1]     0.10823715   0.2085361   0.4664891 4000.000 0.9999784
phi[2,2]     0.26591206   0.3916398   0.6373429 4000.000 0.9993131
phi[2,3]     0.57647391   0.6975200   0.8830145 4000.000 0.9997687

「遷移行列」と「出力行列」を表にしてみた。

• 遷移行列（）

	Sunny	Rainy
Sunny	0.67	0.33
Rainy	0.34	0.66

→約66％で次の日も同じ天気であることを示している。

• 出力行列（）

	walk	shop	clean
Sunny	0.37	0.50	0.13
Rainy	0.14	0.29	0.57

→晴れの日は「掃除」をする確率は低い、雨の日は「散歩」をする確率は低いことを示している。

その他

---

複数人をデータとする場合、以下のような行列を用意すればよい。

(emission  <- rbind(c(1,1,3,1,2,1,2,1,3,1),
                    c(2,2,3,2,2,1,1,2,2,1),
                    c(1,2,2,2,1,1,1,1,1,2)))
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    3    1    2    1    2    1    3     1
[2,]    2    2    3    2    2    1    1    2    2     1
[3,]    1    2    2    2    1    1    1    1    1     2

コードは以下のように修正する。

data {
  int<lower=1> K;  // カテゴリーの数
  int<lower=1> V;  // 単語(word)の数
  int<lower=0> T;  // 時点の数
  int<lower=1> N;  // 人数
  int<lower=1,upper=V> w[N,T]; // 単語(word)
  int<lower=1,upper=K> z[T]; // カテゴリー
  vector<lower=0>[K] alpha;  // 推移(transit)確率の事前確率
  vector<lower=0>[V] beta;   // 単語vを出力(emit)する確率の事前確率
}
parameters {
  simplex[K] theta[K];  // 推移(transit)確率
  simplex[V] phi[K];    // 単語vを出力(emit)する確率
}
model {
  for (k in 1:K)
    theta[k] ~ dirichlet(alpha);
  for (k in 1:K)
    phi[k] ~ dirichlet(beta);
  for (n in 1:N){
    for (t in 1:T)
      w[n,t] ~ categorical(phi[z[t]]);
  }
    for (t in 2:T)
      z[t] ~ categorical(theta[z[t - 1]]);
}

おわりに

---

今回は遷移行列、出力行列の推定を行った。
ただ、遷移行列、出力行列はパラメータに過ぎずこれ自体に大きな意味があるものではない。
これが既知の場合、次は潜在変数の推定に興味がある。

次は、そのあたりをまとめたい。