Stanで隠れマルコフモデル：Taglibro de H：SSブログ

	ブログをはじめるログイン

Stanで隠れマルコフモデル　[統計] [編集]

Stanで、隠れマルコフモデルのパラメータ推定をやってみたので、そのメモです。

データは以下のようなものとしました。

t-1のとき潜在状態1ならば、tのとき、0.8の確率で潜在状態1、0.2の確率で潜在状態2。
t-1のとき潜在状態2ならば、tのとき、0.5の確率で潜在状態1、0.5の確率で潜在状態2。
潜在状態1のとき観測値はNormal(2, 1)、潜在状態2のとき観測値はNormal(6, 1)にしたがう。
以上のような、長さNt=80の時系列が、独立にNs=20だけ観測されている。

以下のRコードでこのようなデータを生成します。

set.seed(1)
Ns <- 20
Nt <- 80
P <- matrix(c(0.8, 0.5, 0.2, 0.5), ncol = 2)
mu <- c(2, 6)
sigma <- c(1, 1)
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] <- rbinom(1, 1, P[z[i, t - 1], 2]) + 1
  y[i, ] <- rnorm(Nt, mu[z[i, ]], sigma[z[i, ]])
}

このようなデータになります。

Stanのコードです。muは潜在状態2のときの方が大きいとして、orderedとして宣言しています。この制約をつけないと、chainによってmuの値が入れ替わります。（追記: この例では、sigmaの値が両方で同じなのでこのコードで推定できましたが、ちがうようならやはり制約をつけないといけないようです。sigma[2] > sigma[1]なら、positive_orderedとしてsigmaを宣言すればよいようです。）また、尤度の計算は、Stan User's Guide and Reference Manual の7.6節にあるforwardアルゴリズムのコードを使用しています。

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

parameters {
  simplex[2] theta[2];  // transition probability
  ordered[2] mu;
  vector<lower=0,upper=100>[2] sigma;
}

model {
  real acc[2];
  real gamma[Nt, 2];
  
  // priors
  for (k in 1:2)
    theta[k] ~ dirichlet(alpha);
  mu ~ normal(0, 100);
  sigma ~ uniform(0, 100);

  // likelihood
  // code from the Stan User's Guide and Reference Manual
  // (section 7.6)
  for (i in 1:Ns) {
    for (k in 1:2)
      gamma[1, k] <- normal_log(y[i, 1], mu[k], sigma[k]);
    for (t in 2:Nt) {
      for (k in 1:2) {
        for (j in 1:2)
          acc[j] <- gamma[t - 1, j] + log(theta[j, k]) +
            normal_log(y[i, t], mu[k], sigma[k]);
        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_test1.stan"
data <- list(Ns = Ns, Nt = Nt, y = y, alpha = c(0.5, 0.5))

fit <- stan(model_file, data = data,
            iter = 2000, warmup = 1000, thin = 1)

データを生成したのと近い値の推定値が得られました。

Inference for Stan model: hmm_test1.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

               mean se_mean   sd     2.5%      25%      50%      75%    97.5%
theta[1,1]     0.79    0.00 0.01     0.76     0.78     0.79     0.80     0.82
theta[1,2]     0.21    0.00 0.01     0.18     0.20     0.21     0.22     0.24
theta[2,1]     0.52    0.00 0.03     0.47     0.50     0.52     0.54     0.57
theta[2,2]     0.48    0.00 0.03     0.43     0.46     0.48     0.50     0.53
mu[1]          1.95    0.00 0.04     1.88     1.92     1.95     1.97     2.02
mu[2]          5.91    0.00 0.07     5.77     5.87     5.91     5.96     6.04
sigma[1]       1.01    0.00 0.03     0.96     0.99     1.01     1.03     1.06
sigma[2]       1.08    0.00 0.05     0.98     1.04     1.08     1.12     1.18
lp__       -3107.41    0.04 1.72 -3111.62 -3108.35 -3107.08 -3106.15 -3105.00
           n_eff Rhat
theta[1,1]  2480    1
theta[1,2]  2480    1
theta[2,1]  2750    1
theta[2,2]  2750    1
mu[1]       2502    1
mu[2]       2645    1
sigma[1]    2661    1
sigma[2]    2493    1
lp__        1602    1

タグ：STAn