用RSTAN拟合poisson-HMM-JAGS模型

沃尔特·西葫芦在他的书中 Hidden Markov Models for Time Series An Introduction Using R ,在第8章第129页,使用R2OpenBUGS调整泊松HMM,然后我展示了代码。我对使用rstan调整这个相同的模型很感兴趣,但是由于我是新使用这个软件包的,所以我不清楚它的语法和任何建议。

数据

dat <- read.table("http://www.hmms-for-time-series.de/second/data/earthquakes.txt")

RJAGS

library(R2jags)
library(rjags)

HMM <- function(){
  for(i in 1:m){
    delta[i] <- 1/m
    v[i] <- 1}
  s[1] ~ dcat(delta[])
  for(i in 2:100){
    s[i] ~ dcat(Gamma[s[i-1],])}
  states[1] ~ dcat(Gamma[s[100],])
  x[1]~dpois(lambda[states[1]])
  for(i in 2:n){
    states[i]~dcat(Gamma[states[i-1],])
    x[i]~dpois(lambda[states[i]])}
  for(i in 1:m){
    tau[i]~dgamma(1,0.08)
    Gamma[i,1:m]~ddirch(v[])}
  lambda[1]<-tau[1]
  for(i in 2:m){
    lambda[i]<-lambda[i-1]+tau[i]}}

x = dat[,2]
n = dim(dat)[1]
m = 2

mod = jags(data = list("x", "n", "m" ), inits = NULL, parameters.to.save = c("lambda","Gamma"), 
           model.file = HMM, n.iter = 10000, n.chains = 1)

输出

mod
Inference for Bugs model at "C:/Users/USER/AppData/Local/Temp/RtmpOkrM6m/model36c8429c5442.txt", fit using jags,
 1 chains, each with 10000 iterations (first 5000 discarded), n.thin = 5
 n.sims = 1000 iterations saved
           mu.vect sd.vect    2.5%     25%     50%     75%   97.5%
Gamma[1,1]   0.908   0.044   0.805   0.884   0.915   0.940   0.971
Gamma[2,1]   0.155   0.071   0.045   0.105   0.144   0.195   0.325
Gamma[1,2]   0.092   0.044   0.029   0.060   0.085   0.116   0.195
Gamma[2,2]   0.845   0.071   0.675   0.805   0.856   0.895   0.955
lambda[1]   15.367   0.763  13.766  14.877  15.400  15.894  16.752
lambda[2]   26.001   1.321  23.418  25.171  25.956  26.843  28.717
deviance   645.351   8.697 630.338 639.359 644.512 650.598 665.405

DIC info (using the rule, pD = var(deviance)/2)
pD = 37.8 and DIC = 683.2
DIC is an estimate of expected predictive error (lower deviance is better).

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

HMM <- '
data{
  int<lower=0> n;    // number of observations (length)
  int<lower=0> x[n]; // observations  
  int<lower=1> m;    // number of hidden states
 }

parameters{
  simplex[m] Gamma[n]; // t.p.m  
  vector[m] lambda;  // mean of poisson ordered
 }

model{
  vector[m] delta[m];
  vector[m] v[m];
  vector[100] s[100];
  vector[n] states[n];
  vector[m] tau;

  for(i in 1:m){
    delta[i] = 1/m;
    v[i] = 1;}
  s[1] ~ categorical(delta[]);

  for(i in 2:100){
    s[i] ~ categorical(Gamma[s[i-1],]);}
  states[1] ~ categorical(Gamma[s[100],]);
  x[1] ~ poisson(lambda[states[1]]);

  for(i in 2:n){
    states[i] ~ categorical(Gamma[states[i-1],]);
  x[i] ~ poisson(lambda[states[i]])};

  for(i in 1:m){
    tau[i] ~ gamma(1,0.08);
    Gamma[i,1:m] ~ dirichlet(v[]);}
  lambda[1] = tau[1];

  for(i in 2:m){
    lambda[i] = lambda[i-1] + tau[i]};}'

data <- list(n = dim(dat)[1], x = dat[,2], m = 2)
system.time(mod2 <- stan(model_code = HMM, data = data, chains = 1, iter = 1000, thin = 4))
mod2

但是,运行stan模型时会发生错误。