Stan：隠れマルコフモデル② - 機械学習・自然言語処理の勉強メモ

はじめに

前回は教師ありHMMについて勉強した。
kento1109.hatenablog.com

今回は教師なしHMMについて勉強する。
今回のデータはセンサーデータのようなものを想定したものを用いる。

例えば、以下のような時系列データがあったとする。
f:id:kento1109:20180623112239p:plain:w500

$t=100$ 時点と $t=200$ 時点で何かしらの状態変化が起きていることは明らかであろう。
これは、あるセンサーが $t=100$ 時点で異常状態になり、 $t=200$ 時点で正常状態に戻ったことを示すデータなどが想定される。
つまり、下図の赤線で観測されない状態（正常・異常）が変化しており、その結果観測データが変わったと考える。
f:id:kento1109:20180623113157p:plain:w500
こう考えると、HMMでモデリング可能だということが分かる。

データ

今回のデータは上記のようなセンサーデータを用いる。
データは以下のようにして作成した。

data <- rpois(n=100, lambda=50)
data <- c(data,rpois(n=10, lambda=100))
data <- c(data,rpois(n=30, lambda=50))
data <- c(data,rpois(n=10, lambda=100))
data <- c(data,rpois(n=50, lambda=50))

赤線の範囲内が異常状態を想定して作成したデータの範囲である。
f:id:kento1109:20180623122543p:plain:w500

今回は観測されるセンサーデータは、非負の整数を仮定してポアソン分布から生成した。
正常状態は $\lambda=50$ 、異常状態は $\lambda=100$ で生成している。

時系列で考えないのであれば、以前に書いたポアソン混合モデルを用いることが出来る。
kento1109.hatenablog.com

今回の場合、各時点でのポアソン混合モデルでも問題ないかもしれないが、一般にはセンサーデータのようなものは現在の状態が直前の状態に依存していると考えるほうが自然である。

今回は

直前が正常状態なら、現在も正常状態である可能性が高い
直前が異常状態なら、現在も異常状態である可能性が高い
ただし、異常状態が続く可能性は正常状態のそれよりも可能性は低い

ことを想定してデータを生成した。
このあたりも推論できるかも確認したい。

結果①

まずは、以下のコードの実行例

stan.dat <- list(K=2,T=length(data),x=data,alpha=c(1.0,1.0))
stan.fit <- stan("hmm_poison.stan", data=stan.dat)

※結論から言うと、このコードではうまく推論が出来ない

data {
  int<lower=1> K;  // num categories
  int<lower=0> T;  // num supervised items
  int<lower=0> x[T]; // timeseries data
  vector<lower=0>[K] alpha;  // transit prior
}
parameters {
  simplex[K] theta[K];  // transit probs
  real lambda[K];       // poisson means
}
model {
  for (k in 1:K) 
    theta[k] ~ dirichlet(alpha);
  { 
    // forward algorithm computes log p(u|...)
    real acc[K];
    real gamma[T,K];
    for (k in 1:K)
      gamma[1,k] = poisson_lpmf(x[1] | lambda[k]);
    for (t in 2:T) {
      for (k in 1:K) {
        for (j in 1:K)
          acc[j] = gamma[t-1,j] + log(theta[j,k]) + poisson_lpmf(x[t] | lambda[k]);
        gamma[t,k] = log_sum_exp(acc);
      }
    }
     target += log_sum_exp(gamma[T]);
  }
}

結果は以下のとおりである。

              mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff
theta[1,1]    0.92    0.04  0.08    0.73    0.88    0.96    0.99    1.00     3
theta[1,2]    0.08    0.04  0.08    0.00    0.01    0.04    0.12    0.27     3
theta[2,1]    0.08    0.04  0.08    0.00    0.01    0.04    0.12    0.27     3
theta[2,2]    0.92    0.04  0.08    0.73    0.88    0.96    0.99    1.00     3
lambda[1]    75.61   18.37 26.02   48.81   49.64   72.23  101.58  105.51     2
lambda[2]    75.59   18.36 26.01   48.78   49.62   72.96  101.58  105.28     2
lp__       -697.50    0.03  1.43 -700.98 -698.23 -697.20 -696.43 -695.69  1977
            Rhat
theta[1,1]  1.63
theta[1,2]  1.63
theta[2,1]  1.57
theta[2,2]  1.57
lambda[1]  16.44
lambda[2]  17.07
lp__        1.00

$\lambda$ がほとんど同じ値であり、Rhatも高い値になっている。
これは、ラベルスイッチングの問題が起きているからである。
詳しくは「Stan モデリング言語: ユーザーガイド・リファレンスマニュアルの『混合分布モデルでのラベルスイッチング』」参照

結果②

では、どうすればよいか。
対策の一つとして紹介されているのが、「パラメータの順序の制約」を追加することである。
$\lambda$ をorder型で指定して、順序に制約を与える。

data {
  int<lower=1> K;  // num categories
  int<lower=0> T;  // num supervised items
  int<lower=0> x[T]; // timeseries data
  vector<lower=0>[K] alpha;  // transit prior
}
parameters {
  simplex[K] theta[K];  // transit probs
  //real lambda[K];       // poisson means
  ordered[K] lambda;
}
model {
  for (k in 1:K) 
    theta[k] ~ dirichlet(alpha);
  { 
    // forward algorithm computes log p(u|...)
    real acc[K];
    real gamma[T,K];
    for (k in 1:K)
      gamma[1,k] = poisson_lpmf(x[1] | lambda[k]);
    for (t in 2:T) {
      for (k in 1:K) {
        for (j in 1:K)
          acc[j] = gamma[t-1,j] + log(theta[j,k]) + poisson_lpmf(x[t] | lambda[k]);
        gamma[t,k] = log_sum_exp(acc);
      }
    }
     target += log_sum_exp(gamma[T]);
  }
}

これの実行結果は以下のとおりである。

              mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
theta[1,1]    0.98    0.00 0.01    0.96    0.98    0.99    0.99    1.00  4000    1
theta[1,2]    0.02    0.00 0.01    0.00    0.01    0.01    0.02    0.04  4000    1
theta[2,1]    0.13    0.00 0.07    0.03    0.08    0.12    0.18    0.30  4000    1
theta[2,2]    0.87    0.00 0.07    0.70    0.82    0.88    0.92    0.97  4000    1
lambda[1]    49.65    0.01 0.52   48.65   49.28   49.64   50.01   50.66  4000    1
lambda[2]   101.70    0.04 2.26   97.38  100.15  101.67  103.21  106.12  4000    1
lp__       -693.51    0.03 1.41 -697.10 -694.22 -693.20 -692.47 -691.73  1940    1