如何在R语言中拟合隐马尔可夫模型并估计p1-p6参数

模型结构明确 #

先把你的问题转化为清晰的HMM数学定义，这是参数识别的核心：

• 隐状态：2个（1=红帽，2=蓝帽）
• 观测值：2个（1=ball_1，2=ball_2）
• 状态转移规则：
  • $P(S_t=1|S_{t-1}=1)=p1$，$P(S_t=2|S_{t-1}=1)=1-p1$
  • $P(S_t=2|S_{t-1}=2)=p2$，$P(S_t=1|S_{t-1}=2)=1-p2$
• 观测生成规则（观测依赖当前隐状态+前一次观测）：
  • 红帽（$S_t=1$）下：
    • $P(O_t=1|O_{t-1}=1, S_t=1)=p3$，$P(O_t=2|O_{t-1}=1, S_t=1)=1-p3$
    • $P(O_t=2|O_{t-1}=2, S_t=1)=p4$，$P(O_t=1|O_{t-1}=2, S_t=1)=1-p4$
  • 蓝帽（$S_t=2$）下：
    • $P(O_t=1|O_{t-1}=1, S_t=2)=p5$，$P(O_t=2|O_{t-1}=1, S_t=2)=1-p5$
    • $P(O_t=2|O_{t-1}=2, S_t=2)=p6$，$P(O_t=1|O_{t-1}=2, S_t=2)=1-p6$

这是带自回归观测的HMM（AR-HMM），标准HMM包默认不支持，需要用depmixS4自定义响应模型实现。

用depmixS4实现步骤 #

1. 数据预处理 #

把观测序列整理成包含前序观测的数据框（因为观测依赖前一次结果）：

# 假设你的观测序列是obs_seq（比如模拟生成的100个值）
obs_df <- data.frame(
  obs = obs_seq,
  lag_obs = c(NA, obs_seq[-length(obs_seq)])  # 新增前一次观测列
)
obs_df <- obs_df[-1, ]  # 去掉无前置观测的第一行

2. 定义自定义HMM模型 #

用depmix()指定2个隐状态，以lag_obs为协变量构建二分类响应模型：

library(depmixS4)

# 构建模型：响应变量obs，协变量lag_obs，2个隐状态，二项logit链接
model <- depmix(obs ~ lag_obs, 
                data = obs_df,
                nstates = 2,
                family = binomial("logit"))

3. 拟合模型并映射参数到p1-p6 #

拟合后提取参数，通过logit逆变换（plogis()）转换为概率：

# EM算法拟合模型
fit_model <- fit(model)

# 提取状态转移参数（对应p1、p2）
trans_params <- getpars(fit_model)[1:4]
p1 <- plogis(trans_params[1])  # 红帽→红帽的概率
p2 <- plogis(trans_params[4])  # 蓝帽→蓝帽的概率

# 提取观测模型参数（对应p3-p6）
obs_params <- getpars(fit_model)[5:8]
# 红帽下的观测参数
intercept1 <- obs_params[1]
coef_lag1 <- obs_params[2]
p3 <- plogis(intercept1 + coef_lag1*1)  # 红帽中前一次取ball_1，当前取ball_1的概率
p4 <- 1 - plogis(intercept1 + coef_lag1*2)  # 红帽中前一次取ball_2，当前取ball_2的概率

# 蓝帽下的观测参数
intercept2 <- obs_params[3]
coef_lag2 <- obs_params[4]
p5 <- plogis(intercept2 + coef_lag2*1)  # 蓝帽中前一次取ball_1，当前取ball_1的概率
p6 <- 1 - plogis(intercept2 + coef_lag2*2)  # 蓝帽中前一次取ball_2，当前取ball_2的概率

# 输出估计结果
cat(sprintf("估计参数：\np1=%.4f, p2=%.4f\np3=%.4f, p4=%.4f\np5=%.4f, p6=%.4f\n",
            p1, p2, p3, p4, p5, p6))

关键注意事项 #

1. 参数识别问题：如果拟合结果不稳定，可能是隐状态可互换（红/蓝帽标签反转），可以给模型设置合理初始参数：

   # 示例初始参数：p1=0.7, p2=0.8, p3=0.6, p4=0.7, p5=0.5, p6=0.6
   init_pars <- c(
     qlogis(0.7), qlogis(0.3),  # 红帽转移的logit值
     qlogis(0.2), qlogis(0.8),  # 蓝帽转移的logit值
     qlogis(0.6) - 1*qlogis(0.6/0.4), qlogis(0.6/0.4),  # 红帽观测参数
     qlogis(0.5) - 1*qlogis(0.5/0.5), qlogis(0.5/0.5)   # 蓝帽观测参数
   )
   model_init <- setpars(model, init_pars)
   fit_model <- fit(model_init)
   

2. 初始观测处理：如果不想丢弃第一行，可以给初始观测单独设定均匀分布（符合你模拟时的初始0.5概率），通过depmix()的initdata参数配置，不过对100样本影响极小。
3. HMM包的局限：HMM包仅支持标准HMM（观测仅依赖当前隐状态），无法处理你的自回归观测逻辑，因此不推荐使用。

内容的提问来源于stack exchange，提问作者stats_noob