HMM R 包 if (d < delta) { : missing value where TRUE/FALSE needed

Question

我正在尝试在 R 中使用

HMM

包。

我想要 4 个隐藏状态，我的观察值范围是 2 到 15。

我可以毫无问题地初始化隐藏模型：

if (!require(HMM, quietly = TRUE)) {
  install.packages("HMM")
  library(HMM)
} else {
  library(HMM)
}

# Load data
url <- "https://raw.githubusercontent.com/luancvieira/HMM/main/ottawa_2010-2012.csv"
df <- read.csv(url)

observed_data <- df$AvgTemperature

# Define number and names of states
n_states <- 4
state_names <- paste0("State", 1:n_states)

# Sorting symbols
symbol_names <- as.character(sort(unique(observed_data)))
observed_data <- as.character(observed_data)

# Initialize the HMM model with random probabilities
start_probs <- runif(n_states)
trans_probs <- matrix(runif(n_states * n_states),
               nrow = n_states, ncol = n_states)
emission_probs <- matrix(runif(n_states * length(symbol_names)), 
                  nrow = n_states, ncol = length(symbol_names))

# Normalize rows to ensure probabilities sum to 1
start_probs <- start_probs / sum(start_probs)
trans_probs <- trans_probs / rowSums(trans_probs)
emission_probs <- emission_probs / rowSums(emission_probs)

# Initialize the HMM model
hmm_model <- initHMM(States = state_names,
                     Symbols = symbol_names,
                     startProbs = start_probs,
                     transProbs = trans_probs,
                     emissionProbs = emission_probs)

# Print the initialized model
print(hmm_model)

但是，当我尝试运行它时：

bw = baumWelch(hmm = hmm_model,observation = observed_data,
               maxIterations = 100, delta = 0.001)

我明白了

Error in if (d < delta) { : missing value where TRUE/FALSE needed.

如果我将迭代次数减少到 10 次，它运行时不会出现问题，但 10 次迭代不足以收敛。 delta 是算法的停止标准，如下所示：https://cran.r-project.org/web/packages/HMM/HMM.pdf。 delta 的默认值为 1e-9，因此即使没有指定 delta，在运行更多迭代时仍然会返回错误。

Answer 1

很难确切知道为什么会出现错误，但我的猜测是拟合算法没有收敛。这可能是因为：

1. 需要用太少的数据来估计太多的参数。您尝试从大约 1000 个数据点估计 12 个转移概率和大约 50 个发射概率。对于包含潜在组件的模型来说，这是一个难题。

2. 初始参数值不合适。考虑到估计所有参数所需的优化的复杂性，提供仔细选择的起始参数值而不是随机值至关重要。您需要考虑每个符号在每个状态中的预期比例（例如，状态 1 的“2”、“3”和“4”概率较高，但其他符号的概率较低等），并且相应地定义排放概率。同样，可以根据状态持久性的期望来初始化转移概率（即对角线元素可能接近 1）。

话虽这么说，对于温度等连续数据，最好使用连续分布（例如正态分布、伽玛分布等）。这意味着您不必将温度值视为定性符号，而是可以对连续温度测量的分布进行建模。我认为 HMM 包不允许连续分发，所以我建议使用另一个包，例如 depmixS4 或 hmmTMB。我使用下面的 hmmTMB 包提供代码。

模型拟合

library(hmmTMB)

# Load data
url <- "https://raw.githubusercontent.com/luancvieira/HMM/main/ottawa_2010-2012.csv"
df <- read.csv(url)

# Create hidden state model (4 states)
hid <- MarkovChain$new(data = df, 
                       n_states = 4)

# Choose observation distribution (normal)
dists <- list(AvgTemperature = "norm")
# Choose initial parameter values for optimisation
par0 <- list(AvgTemperature = list(mean = c(3, 6, 10, 13),
                                   sd = c(1, 1, 1, 1)))
# Create observation model
obs <- Observation$new(data = df, 
                       n_states = 4, 
                       dists = dists, 
                       par = par0)

# Create and fit HMM
hmm <- HMM$new(obs = obs, hid = hid)
hmm$fit()

预估参数

拟合模型后，可以找到隐藏状态过程的估计参数（转移概率）和观测模型（每个状态的温度平均值和SD）：

> round(hid$tpm()[,,1], 2)
        state 1 state 2 state 3 state 4
state 1    0.96    0.04    0.00    0.00
state 2    0.04    0.91    0.04    0.00
state 3    0.00    0.04    0.93    0.03
state 4    0.00    0.00    0.03    0.97

> round(obs$par()[,,1], 2)
                    state 1 state 2 state 3 state 4
AvgTemperature.mean    3.47    6.75    9.90   13.29
AvgTemperature.sd      1.34    1.19    1.27    1.07

绘图结果

您还可以绘制估计分布，以及按最可能的状态序列着色的数据时间序列图（使用维特比算法预测）。

# Plot estimated state-dependent distributions
hmm$plot_dist("AvgTemperature")

# Plot most likely state sequence
hmm$plot_ts("AvgTemperature")