使用R语言HMM包时baumWelch函数报错：缺失TRUE/FALSE所需值

问题重现 #

使用R语言的HMM包构建含4个隐状态的模型，观测值范围为2-15，初始化过程无异常，但运行baumWelch函数时，当迭代次数设置为100（或较大值）时出现如下报错：

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

减少迭代次数到10时可正常运行，但迭代次数不足无法让模型收敛，即使使用默认的delta=1e-9仍会报错。

初始化代码如下：

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)

报错的baumWelch调用代码：

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

报错原因 #

该错误本质是迭代过程中计算的对数似然差值d出现了NA值，通常由以下原因导致：

• 数值下溢：随机初始化的概率矩阵可能出现极小值，后续计算中多次相乘后数值趋近于0，取对数时产生NaN。
• 极端初始参数：完全随机的初始概率可能导致模型在迭代中陷入不稳定状态，似然值计算出现异常。
• 零概率问题：发射或转移矩阵中存在0概率值，计算log概率时直接产生-Inf，最终导致差值计算出现NA。

解决方案 #

1. 带平滑的参数初始化 #

在初始化概率矩阵时，给每个元素添加极小的平滑项（如1e-8），避免出现0概率，同时保证归一化后概率和为1：

# 带平滑的初始概率初始化
epsilon <- 1e-8

start_probs <- runif(n_states) + epsilon
start_probs <- start_probs / sum(start_probs)

trans_probs <- matrix(runif(n_states * n_states) + epsilon,
                      nrow = n_states, ncol = n_states)
trans_probs <- trans_probs / rowSums(trans_probs)

emission_probs <- matrix(runif(n_states * length(symbol_names)) + epsilon, 
                         nrow = n_states, ncol = length(symbol_names))
emission_probs <- emission_probs / rowSums(emission_probs)

2. 基于数据分布初始化发射概率 #

完全随机的发射概率可能偏离数据真实分布，可基于观测数据的分位数分组来初始化，让初始参数更合理：

# 基于观测数据分位数初始化发射概率
# 将观测值按4个分位数分组（对应4个状态）
observed_numeric <- as.numeric(observed_data)
quantiles <- quantile(observed_numeric, probs = seq(0, 1, length.out = n_states + 1))

emission_probs <- matrix(0, nrow = n_states, ncol = length(symbol_names))
for (i in 1:n_states) {
  # 找到当前状态对应的观测值范围
  lower <- quantiles[i]
  upper <- quantiles[i+1]
  # 计算该范围内每个观测值的频率
  mask <- observed_numeric >= lower & observed_numeric <= upper
  freq <- table(observed_numeric[mask])
  # 将频率映射到symbol_names的位置
  for (sym in names(freq)) {
    col_idx <- which(symbol_names == sym)
    emission_probs[i, col_idx] <- freq[sym]
  }
  # 添加平滑项并归一化
  emission_probs[i, ] <- emission_probs[i, ] + epsilon
  emission_probs[i, ] <- emission_probs[i, ] / sum(emission_probs[i, ])
}

3. 分阶段迭代 #

先运行少量迭代得到相对稳定的模型，再以此为初始值继续迭代，避免一次性迭代过多步出现数值异常：

# 第一阶段：少量迭代得到稳定模型
bw_stage1 <- baumWelch(hmm = hmm_model, observation = observed_data,
                       maxIterations = 20, delta = 0.001)

# 第二阶段：用第一阶段的结果继续迭代
bw_final <- baumWelch(hmm = bw_stage1$hmm, observation = observed_data,
                      maxIterations = 80, delta = 0.001)

4. 跟踪迭代过程 #

开启trace=TRUE参数，查看每一步的似然值变化，定位异常出现的迭代步骤：

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

修改后的完整示例代码 #

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
observed_numeric <- as.numeric(observed_data)

# 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)

# 初始化参数：带平滑+基于数据的发射概率
epsilon <- 1e-8

# 初始状态概率
start_probs <- runif(n_states) + epsilon
start_probs <- start_probs / sum(start_probs)

# 转移概率矩阵
trans_probs <- matrix(runif(n_states * n_states) + epsilon,
                      nrow = n_states, ncol = n_states)
trans_probs <- trans_probs / rowSums(trans_probs)

# 基于分位数的发射概率初始化
quantiles <- quantile(observed_numeric, probs = seq(0, 1, length.out = n_states + 1))
emission_probs <- matrix(epsilon, nrow = n_states, ncol = length(symbol_names))

for (i in 1:n_states) {
  lower <- quantiles[i]
  upper <- quantiles[i+1]
  mask <- observed_numeric >= lower & observed_numeric <= upper
  freq <- table(observed_numeric[mask])
  
  for (sym in names(freq)) {
    col_idx <- which(symbol_names == sym)
    emission_probs[i, col_idx] <- freq[sym] + epsilon
  }
  
  emission_probs[i, ] <- emission_probs[i, ] / sum(emission_probs[i, ])
}

# 初始化HMM模型
hmm_model <- initHMM(States = state_names,
                     Symbols = symbol_names,
                     startProbs = start_probs,
                     transProbs = trans_probs,
                     emissionProbs = emission_probs)

# 分阶段迭代
bw_stage1 <- baumWelch(hmm = hmm_model, observation = observed_data,
                       maxIterations = 20, delta = 0.001)
bw_final <- baumWelch(hmm = bw_stage1$hmm, observation = observed_data,
                      maxIterations = 80, delta = 0.001, trace = TRUE)

# 查看最终模型
print(bw_final$hmm)

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