在 R 中学习隐马尔可夫模型
Learning hidden markov model in R
隐马尔可夫模型 (HMM) 是您观察一系列观察结果但不知道模型生成观察结果所经历的状态序列的模型。隐马尔可夫模型的分析试图从观察到的数据中恢复隐藏状态的序列。
我有包含观察值和隐藏状态（观察值是连续值）的数据，其中隐藏状态由专家标记。我想训练一个 HMM，它能够 - 基于（以前看不见的）观察序列 - 恢复相应的隐藏状态。
是否有任何 R 包可以做到这一点？研究现有包（depmixS4、HMM、seqHMM - 仅用于分类数据）允许您仅指定一些隐藏状态。
编辑：
示例：
data.tagged.by.expert = data.frame(
hidden.state = c("Wake", "REM", "REM", "NonREM1", "NonREM2", "REM", "REM", "Wake"),
sensor1 = c(1,1.2,1.2,1.3,4,2,1.78,0.65),
sensor2 = c(7.2,5.3,5.1,1.2,2.3,7.5,7.8,2.1),
sensor3 = c(0.01,0.02,0.08,0.8,0.03,0.01,0.15,0.45)
)
data.newly.measured = data.frame(
sensor1 = c(2,3,4,5,2,1,2,4,5,8,4,6,1,2,5,3,2,1,4),
sensor2 = c(2.1,2.3,2.2,4.2,4.2,2.2,2.2,5.3,2.4,1.0,2.5,2.4,1.2,8.4,5.2,5.5,5.2,4.3,7.8),
sensor3 = c(0.23,0.25,0.23,0.54,0.36,0.85,0.01,0.52,0.09,0.12,0.85,0.45,0.26,0.08,0.01,0.55,0.67,0.82,0.35)
)
我想创建一个具有离散时间 t 的 HMM，其中随机变量 x(t) 表示时间 t, x(t) 
{"Wake", "REM", "NonREM1", "NonREM2"}, 和3个连续随机变量sensor1(t), sensor2(t), sensor3(t)代表时间的观测值t。
model.hmm = learn.model(data.tagged.by.user)
然后我想使用创建的模型来估计负责新测量观察的隐藏状态
hidden.states = estimate.hidden.states(model.hmm, data.newly.measured)
数据（training/testing）
为了能够运行学习朴素贝叶斯分类器的方法，我们需要更长的数据集
states = c("NonREM1", "NonREM2", "NonREM3", "REM", "Wake")
artificial.hypnogram = rep(c(5,4,1,2,3,4,5), times = c(40,150,200,300,50,90,30))
data.tagged.by.expert = data.frame(
hidden.state = states[artificial.hypnogram],
sensor1 = log(artificial.hypnogram) + runif(n = length(artificial.hypnogram), min = 0.2, max = 0.5),
sensor2 = 10*artificial.hypnogram + sample(c(-8:8), size = length(artificial.hypnogram), replace = T),
sensor3 = sample(1:100, size = length(artificial.hypnogram), replace = T)
)
hidden.hypnogram = rep(c(5,4,1,2,4,5), times = c(10,10,15,10,10,3))
data.newly.measured = data.frame(
sensor1 = log(hidden.hypnogram) + runif(n = length(hidden.hypnogram), min = 0.2, max = 0.5),
sensor2 = 10*hidden.hypnogram + sample(c(-8:8), size = length(hidden.hypnogram), replace = T),
sensor3 = sample(1:100, size = length(hidden.hypnogram), replace = T)
)
解决方案
在解决方案中，我们使用了维特比算法——结合朴素贝叶斯分类器。
在每个时钟时间t，一个隐马尔可夫模型由
组成
一个未观察到的状态（在这种情况下表示为 hidden.state）采用有限数量的状态
states = c("NonREM1", "NonREM2", "NonREM3", "REM", "Wake")
一组观测变量（本例中为传感器 1、传感器 2、传感器 3）
转移矩阵
根据转移概率分布进入新状态
 "P(hidden.state_{t+1}|hidden.state_t)")
（转换矩阵）。这可以很容易地从 data.tagged.by.expert 计算出来，例如使用
library(markovchain)
emit_p <- markovchainFit(data.tagged.by.expert$hidden.state)$estimate
发射矩阵
每次转换后，根据条件概率分布 "P(sensor_{i,t}|hidden.state_t = H)")
（发射矩阵）生成一个观察值（sensor_i），这取决于仅 hidden.state 的当前状态 H。我们将用朴素贝叶斯分类器替换发射矩阵。
library(caret)
library(klaR)
library(e1071)
model = train(hidden.state ~ .,
data = data.tagged.by.expert,
method = 'nb',
trControl=trainControl(method='cv',number=10)
)
维特比算法
为了解决这个问题，我们使用 Viterbi algorithm，"Wake" 状态的初始概率为 1，否则为 0。 （我们希望患者在实验开始时是清醒的）
# we expect the patient to be awake in the beginning
start_p = c(NonREM1 = 0,NonREM2 = 0,NonREM3 = 0, REM = 0, Wake = 1)
# Naive Bayes model
model_nb = model$finalModel
# the observations
observations = data.newly.measured
nObs <- nrow(observations) # number of observations
nStates <- length(states) # number of states
# T1, T2 initialization
T1 <- matrix(0, nrow = nStates, ncol = nObs) #define two 2-dimensional tables
row.names(T1) <- states
T2 <- T1
Byj <- predict(model_nb, newdata = observations[1,])$posterior
# init first column of T1
for(s in states)
T1[s,1] = start_p[s] * Byj[1,s]
# fill T1 and T2 tables
for(j in 2:nObs) {
Byj <- predict(model_nb, newdata = observations[j,])$posterior
for(s in states) {
res <- (T1[,j-1] * emit_p[,s]) * Byj[1,s]
T2[s,j] <- states[which.max(res)]
T1[s,j] <- max(res)
}
}
# backtract best path
result <- rep("", times = nObs)
result[nObs] <- names(which.max(T1[,nObs]))
for (j in nObs:2) {
result[j-1] <- T2[result[j], j]
}
# show the result
result
# show the original artificial data
states[hidden.hypnogram]
参考资料
要阅读有关该问题的更多信息，请参阅 Vomlel Jiří，Kratochvíl Václav：用于睡眠阶段分类的动态贝叶斯网络，第 11 届不确定性处理研讨会 (WUPES'18) 的论文集，第 12 页。 205-215，Eds：Kratochvíl Václav，Vejnarová Jiřina，不确定性处理研讨会 (WUPES'18)，（Třeboň，CZ，2018/06/06）[2018] Download
隐马尔可夫模型 (HMM) 是您观察一系列观察结果但不知道模型生成观察结果所经历的状态序列的模型。隐马尔可夫模型的分析试图从观察到的数据中恢复隐藏状态的序列。
我有包含观察值和隐藏状态（观察值是连续值）的数据，其中隐藏状态由专家标记。我想训练一个 HMM，它能够 - 基于（以前看不见的）观察序列 - 恢复相应的隐藏状态。
是否有任何 R 包可以做到这一点？研究现有包（depmixS4、HMM、seqHMM - 仅用于分类数据）允许您仅指定一些隐藏状态。
编辑：
示例：
data.tagged.by.expert = data.frame(
hidden.state = c("Wake", "REM", "REM", "NonREM1", "NonREM2", "REM", "REM", "Wake"),
sensor1 = c(1,1.2,1.2,1.3,4,2,1.78,0.65),
sensor2 = c(7.2,5.3,5.1,1.2,2.3,7.5,7.8,2.1),
sensor3 = c(0.01,0.02,0.08,0.8,0.03,0.01,0.15,0.45)
)
data.newly.measured = data.frame(
sensor1 = c(2,3,4,5,2,1,2,4,5,8,4,6,1,2,5,3,2,1,4),
sensor2 = c(2.1,2.3,2.2,4.2,4.2,2.2,2.2,5.3,2.4,1.0,2.5,2.4,1.2,8.4,5.2,5.5,5.2,4.3,7.8),
sensor3 = c(0.23,0.25,0.23,0.54,0.36,0.85,0.01,0.52,0.09,0.12,0.85,0.45,0.26,0.08,0.01,0.55,0.67,0.82,0.35)
)
我想创建一个具有离散时间 t 的 HMM，其中随机变量 x(t) 表示时间 t, x(t) {"Wake", "REM", "NonREM1", "NonREM2"}, 和3个连续随机变量sensor1(t), sensor2(t), sensor3(t)代表时间的观测值t。
model.hmm = learn.model(data.tagged.by.user)
然后我想使用创建的模型来估计负责新测量观察的隐藏状态
hidden.states = estimate.hidden.states(model.hmm, data.newly.measured)
数据（training/testing）
为了能够运行学习朴素贝叶斯分类器的方法，我们需要更长的数据集
states = c("NonREM1", "NonREM2", "NonREM3", "REM", "Wake")
artificial.hypnogram = rep(c(5,4,1,2,3,4,5), times = c(40,150,200,300,50,90,30))
data.tagged.by.expert = data.frame(
hidden.state = states[artificial.hypnogram],
sensor1 = log(artificial.hypnogram) + runif(n = length(artificial.hypnogram), min = 0.2, max = 0.5),
sensor2 = 10*artificial.hypnogram + sample(c(-8:8), size = length(artificial.hypnogram), replace = T),
sensor3 = sample(1:100, size = length(artificial.hypnogram), replace = T)
)
hidden.hypnogram = rep(c(5,4,1,2,4,5), times = c(10,10,15,10,10,3))
data.newly.measured = data.frame(
sensor1 = log(hidden.hypnogram) + runif(n = length(hidden.hypnogram), min = 0.2, max = 0.5),
sensor2 = 10*hidden.hypnogram + sample(c(-8:8), size = length(hidden.hypnogram), replace = T),
sensor3 = sample(1:100, size = length(hidden.hypnogram), replace = T)
)
解决方案
在解决方案中，我们使用了维特比算法——结合朴素贝叶斯分类器。
在每个时钟时间t，一个隐马尔可夫模型由
组成一个未观察到的状态（在这种情况下表示为 hidden.state）采用有限数量的状态
states = c("NonREM1", "NonREM2", "NonREM3", "REM", "Wake")一组观测变量（本例中为传感器 1、传感器 2、传感器 3）
转移矩阵
根据转移概率分布进入新状态
（转换矩阵）。这可以很容易地从 data.tagged.by.expert 计算出来，例如使用
library(markovchain)
emit_p <- markovchainFit(data.tagged.by.expert$hidden.state)$estimate
发射矩阵
每次转换后，根据条件概率分布（发射矩阵）生成一个观察值（sensor_i），这取决于仅 hidden.state 的当前状态 H。我们将用朴素贝叶斯分类器替换发射矩阵。
library(caret)
library(klaR)
library(e1071)
model = train(hidden.state ~ .,
data = data.tagged.by.expert,
method = 'nb',
trControl=trainControl(method='cv',number=10)
)
维特比算法
为了解决这个问题，我们使用 Viterbi algorithm，"Wake" 状态的初始概率为 1，否则为 0。 （我们希望患者在实验开始时是清醒的）
# we expect the patient to be awake in the beginning
start_p = c(NonREM1 = 0,NonREM2 = 0,NonREM3 = 0, REM = 0, Wake = 1)
# Naive Bayes model
model_nb = model$finalModel
# the observations
observations = data.newly.measured
nObs <- nrow(observations) # number of observations
nStates <- length(states) # number of states
# T1, T2 initialization
T1 <- matrix(0, nrow = nStates, ncol = nObs) #define two 2-dimensional tables
row.names(T1) <- states
T2 <- T1
Byj <- predict(model_nb, newdata = observations[1,])$posterior
# init first column of T1
for(s in states)
T1[s,1] = start_p[s] * Byj[1,s]
# fill T1 and T2 tables
for(j in 2:nObs) {
Byj <- predict(model_nb, newdata = observations[j,])$posterior
for(s in states) {
res <- (T1[,j-1] * emit_p[,s]) * Byj[1,s]
T2[s,j] <- states[which.max(res)]
T1[s,j] <- max(res)
}
}
# backtract best path
result <- rep("", times = nObs)
result[nObs] <- names(which.max(T1[,nObs]))
for (j in nObs:2) {
result[j-1] <- T2[result[j], j]
}
# show the result
result
# show the original artificial data
states[hidden.hypnogram]
参考资料
要阅读有关该问题的更多信息，请参阅 Vomlel Jiří，Kratochvíl Václav：用于睡眠阶段分类的动态贝叶斯网络，第 11 届不确定性处理研讨会 (WUPES'18) 的论文集，第 12 页。 205-215，Eds：Kratochvíl Václav，Vejnarová Jiřina，不确定性处理研讨会 (WUPES'18)，（Třeboň，CZ，2018/06/06）[2018] Download