Stanで隠れマルコフモデル (2) うまくいかなかった例

今回は3状態で、出力も離散値としました。が、単純なモデルではうまくいきませんでした。

潜在状態は{1, 2, 3}の3状態をとり、観測値も{1, 2, 3}の3状態です。遷移行列は、出力行列はとしました。

以下のようなRコードで、仮想データを生成しました。Nsが時系列の数、Ntが各時系列の長さ、zが潜在状態、yが観測値です。

set.seed(1)
Ns <- 20
Nt <- 40
Mt <- matrix(c(0.7, 0.3, 0.0, 0.2, 0.6, 0.0, 0.1, 0.1, 1.0), ncol = 3)
Me <- matrix(c(0.8, 0.1, 0.1, 0.1, 0.8, 0.1, 0.1, 0.1, 0.8), ncol = 3)

y <- z <- matrix(0, nrow = Ns, ncol = Nt)
for (i in 1:Ns) {
  z[i, 1] <- 1
  for (t in 2:Nt)
    z[i, t] <- grep(1, rmultinom(1, 1, Mt[z[i, t - 1], ]))
  for (t in 1:Nt)      
    y[i, t] <- grep(1, rmultinom(1, 1, Me[z[i, t], ]))
}

Stanのコードは以下のとおりです。前回とおなじく、Stan User's Guide and Reference Manual のforwardアルゴリズムをつかっています。

data {
  int<lower=1> Ns;              // number of series
  int<lower=1> Nt;              // number of occasions
  int y[Ns, Nt];                // observed values
  vector<lower=0>[3] alpha;     // transition prior
  vector<lower=0>[3] beta;      // emission prior
}

parameters {
  simplex[3] theta[3];
  simplex[3] phi[3];
}

model {
  real acc[3];
  real gamma[Nt, 3];
  
  // priors
  for (k in 1:3)
    theta[k] ~ dirichlet(alpha);
  for (k in 1:3)
    phi[k] ~ dirichlet(beta);
  
  // likelihood
  for (i in 1:Ns) {
    for (k in 1:3)
      gamma[1, k] <- log(phi[k, y[i, 1]]);
    for (t in 2:Nt) {
      for (k in 1:3) {
        for (j in 1:3)
          acc[j] <- gamma[t - 1, j] + log(theta[j, k]) +
                    log(phi[k, y[i, t]]);
        gamma[t, k] <- log_sum_exp(acc);
      }
    }
    increment_log_prob(log_sum_exp(gamma[Nt]));
  }
}

RStanを呼び出すコードです。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

model_file <- "hmm_test2.stan"
data <- list(Ns = Ns, Nt = Nt, y = y,
             alpha = rep(1/3, 3), beta = rep(1/3, 3))
fit <- stan(model_file, data = data,
            iter = 7000, warmup = 2000, thin = 5)

しかし結果は収束しませんでした。マニュアルにも“The resulting posteriors are typically extremely multimodal.”とあるのですが、そのような雰囲気です。

> traceplot(fit, pars = "theta")

> traceplot(fit, pars = "phi")

おそらく、情報をくわえるか、制約をつけるかしないといけないのでしょう。

タグ：STAn