如何在R中绘制非齐次隐马尔可夫模型(NHMM)各状态的分布

你拟合的模型包含4个高斯分布的响应变量、4个隐状态，可按以下步骤提取参数并绘图：

步骤1：提取各状态的高斯分布参数 #

depmixS4拟合结果中高斯响应的标准差默认以对数形式存储，需转换为原始值：

library(ggplot2)
library(depmixS4)

n_states <- NHMM.fit@nstates
n_responses <- length(NHMM.fit@response)

# 构建各状态、各响应的参数表
state_params <- data.frame()
for (state in 1:n_states) {
  for (resp in 1:n_responses) {
    mean_val <- NHMM.fit@response[[state]][[resp]]@parameters$coefficients
    sd_val <- exp(NHMM.fit@response[[state]][[resp]]@parameters$sd)
    resp_name <- names(NHMM.fit@response[[state]])[resp]
    state_params <- rbind(state_params, data.frame(
      state = factor(state),
      response = resp_name,
      mean = mean_val,
      sd = sd_val
    ))
  }
}

步骤2：多响应分面绘制各状态分布 #

该方法可以同时展示4个响应变量下不同状态的概率密度分布：

# 生成密度计算用的网格数据
plot_data <- data.frame()
for (i in 1:nrow(state_params)) {
  param_row <- state_params[i, ]
  x_vals <- seq(param_row$mean - 3*param_row$sd, param_row$mean + 3*param_row$sd, length.out = 100)
  density_vals <- dnorm(x_vals, mean = param_row$mean, sd = param_row$sd)
  plot_data <- rbind(plot_data, data.frame(
    state = param_row$state,
    response = param_row$response,
    x = x_vals,
    density = density_vals
  ))
}

# 绘图
ggplot(plot_data, aes(x = x, y = density, color = state, fill = state)) +
  geom_line(linewidth = 1) +
  geom_area(alpha = 0.2, position = "identity") +
  facet_wrap(~response, scales = "free") +
  labs(x = "变量取值", y = "概率密度", color = "隐状态", fill = "隐状态") +
  theme_bw()

可选：单变量快速绘图（base R版本） #

如果仅需快速查看单个响应的状态分布，可使用基础绘图函数实现：

# 以第一个响应变量Clust1.North为例
resp_id <- 1
state_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3")

# 确定x轴范围
x_lim <- range(sapply(1:n_states, function(s) {
  m <- NHMM.fit@response[[s]][[resp_id]]@parameters$coefficients
  sd <- exp(NHMM.fit@response[[s]][[resp_id]]@parameters$sd)
  c(m - 3*sd, m + 3*sd)
}))

# 绘制画布
plot(NA, xlim = x_lim, ylim = c(0, max(sapply(1:n_states, function(s) {
  m <- NHMM.fit@response[[s]][[resp_id]]@parameters$coefficients
  sd <- exp(NHMM.fit@response[[s]][[resp_id]]@parameters$sd)
  dnorm(m, mean = m, sd = sd)
}))), xlab = "Clust1.North取值", ylab = "概率密度")

# 逐个添加状态分布曲线
for (s in 1:n_states) {
  m <- NHMM.fit@response[[s]][[resp_id]]@parameters$coefficients
  sd <- exp(NHMM.fit@response[[s]][[resp_id]]@parameters$sd)
  curve(dnorm(x, mean = m, sd = sd), add = TRUE, col = state_colors[s], lwd = 2)
}

# 添加图例
legend("topright", legend = paste0("状态", 1:n_states), col = state_colors, lwd = 2)

可选：双变量联合状态分布绘制 #

如果需要查看两个响应变量在各状态下的联合分布，可使用以下代码：

# 以Clust1.North和Clust2.North两个响应为例
target_resps <- c("Clust1.North", "Clust2.North")
target_params <- subset(state_params, response %in% target_resps)
# 转换为宽表
wide_params <- reshape(target_params, idvar = "state", timevar = "response", direction = "wide")

joint_data <- data.frame()
for (s in 1:n_states) {
  sp <- subset(wide_params, state == s)
  # 生成二维网格
  x_vals <- seq(sp$mean.Clust1.North - 3*sp$sd.Clust1.North, sp$mean.Clust1.North + 3*sp$sd.Clust1.North, length.out = 50)
  y_vals <- seq(sp$mean.Clust2.North - 3*sp$sd.Clust2.North, sp$mean.Clust2.North + 3*sp$sd.Clust2.North, length.out = 50)
  grid <- expand.grid(x = x_vals, y = y_vals)
  # 计算联合密度（depmixS4默认多响应变量独立）
  grid$density <- dnorm(grid$x, sp$mean.Clust1.North, sp$sd.Clust1.North) * dnorm(grid$y, sp$mean.Clust2.North, sp$sd.Clust2.North)
  grid$state <- factor(s)
  joint_data <- rbind(joint_data, grid)
}

# 分面绘制热力图
ggplot(joint_data, aes(x = x, y = y, fill = density)) +
  geom_tile() +
  scale_fill_viridis_c() +
  facet_wrap(~state) +
  labs(x = "Clust1.North", y = "Clust2.North", fill = "联合密度") +
  theme_bw()

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