はじめに
今回は隠れマルコフモデルを
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]]);
}
おわりに
今回は遷移行列、出力行列の推定を行った。
ただ、遷移行列、出力行列はパラメータに過ぎずこれ自体に大きな意味があるものではない。
これが既知の場合、次は潜在変数の推定に興味がある。
次は、そのあたりをまとめたい。