训练HMM模型后如何输出rjags上的隐藏状态序列？

Question

我对 HMM 和 JAGS 完全陌生。我的问题是如何根据训练好的模型结构和参数获得隐藏状态序列。 具体来说，我使用rjags训练了一个隐马尔可夫模型，并保存了模型参数的样本结果，其中包含59个参数，长度为1000（30000 iters，thining = 30）。 现在，我有了以 JAGS 的方式正确编写的模型结构，以及更新的参数（30000 个老化、30000 个样本迭代和 30 个细化）。问题是，怎样才能得到对应的隐藏状态呢？ 我从事机器学习有一段时间了，所以我的想法是首先加载结构和参数，然后输入可观察状态的测试序列，然后我就会得到我想要的。但似乎有所不同。 到目前为止，我没有找到这个问题的任何样本，即使在 JAGS 的文档中也是如此。 有人知道这件事吗？多谢。 我的训练代码如下：

# -------------------  model ------------------------#

load('./data/try_3_gtx_3d.RData')

library(rjags)
library(mcmcplots)

data_1 <- list(Y = three_dim_smc_data-1, 
             M = dim(team_id_unique_ordered)[1], 
             N = one_dim_sample_count, 
             T = two_dim_sample_length, 
             gp = three_dim_gp_data, 
             sm = three_dim_sm_data)   # Y should be 0 and 1


inits_HMM <- list( beta0.p = runif(2, -10, 10))

parameters <- c("beta0", "PAD", "PDA", "bAD", "bDA", "sigmaAD", "sigmaDA", 
                "sigmab", "betagp", 
                "b", "delta0",  'betasm')

result <- jags.model("model/longitudinal_model_gp_sm_try1_gtx.txt", data_1, inits_HMM, n.chains = 3, n.adapt = 0)

update(result, n.iter = 30000)

rsamps_relience_pg_sm_try2 <- coda.samples(result, parameters, n.iter = 30000, 
                                           thin = 30)

save(rsamps_relience_pg_sm_try2, file = './results/rsamps_resilience_pg_sm_try3.RData')

mcmcplot(rsamps_relience_pg_sm_try2)

summary(rsamps_relience_pg_sm_try1)$statistics[-(c(3 : 317)),]
summary(rsamps_relience_pg_sm_try1)$quantiles[-(c(3 : 317)),]

如果需要其他材料，请告知我。非常感谢！

我正在尝试根据模型结构和训练有素的参数输出与给定的测试可观察状态序列相对应的隐藏状态序列。

Answer 1

一种方法是根据您拥有的参数估计，使用维特比算法来获取新数据最可能的状态序列。有关维特比算法和状态解码的一般问题的背景，请参阅 Rabiner (1989)、隐马尔可夫模型教程（第 III.B 节）概述的“HMM 的三个基本问题”中的问题 2 ）。最常见的是，这是使用 HMM 参数的点估计来完成的（例如，您的情况下的后验平均值）。

如果您不想自己实现维特比算法，可以使用许多现有的 R 包。例如，下面的代码显示了 depmixS4 中的工作流程，其中将 2 状态高斯 HMM 拟合到包中包含的日志返回时间序列。步骤是：

1. 在训练数据上拟合模型
2. 为测试数据创建新模型
3. 使用

   setpars()

从训练模型中复制估计参数
4. 使用

   viterbi()

查找最可能的状态序列

# I'm using a financial data set provided in depmixS4
library(depmixS4)
data(sp500)

# Split into training and test sets
train <- sp500[1:500,]
test <- sp500[-(1:500),]

# Create model for training data
par0 <- c(0, 0.01, 0, 0.1)
mod_train <- depmix(logret ~ 1, 
              data = train, 
              nstates = 2, 
              family = gaussian(), 
              respstart = par0)

# Fit model on training data
fit_train <- fit(mod_train)

# Create model for new/test data
mod_new <- depmix(logret ~ 1, 
                  data = test, 
                  nstates = 2, 
                  family = gaussian())

# Copy estimated parameters from trained model
mod_new <- setpars(mod_new, getpars(mod_train))

# Get most likely state sequence for new data
vit_new <- viterbi(mod_new)

# Plot log-returns in new data, coloured by states
plot(test$logret, col = vit_new$state, 
     type = "h", ylab = "log-return",
     main = "Decoded states for new data")
points(test$logret, col = vit_new$state, 
       pch = 20, cex = 0.5)