# Chapter 10 隐马尔可夫模型

``````rm(list = ls())
setClass("HMM", representation(A = "matrix",
B = "matrix",
PI = "numeric",
state_set = "character",
observation_set = "character",
T_length = "integer",
state_series = "character",
observation_series = "character"))

getClass("HMM")``````
``````## Class "HMM" [in ".GlobalEnv"]
##
## Slots:
##
## Name:                   A                  B                 PI
## Class:             matrix             matrix            numeric
##
## Name:           state_set    observation_set           T_length
## Class:          character          character            integer
##
## Name:        state_series observation_series
## Class:          character          character``````
``````HMM_createobser <- function(A, B, PI, T_length,
state_set, observation_set,
#默认是最大概率的方式生成状态序列和观测序列
fun = which.max) {

stopifnot(nrow(A) == length(state_set) &
ncol(A) == length(state_set))
stopifnot(nrow(B) == length(state_set) &
ncol(B) == length(observation_set))
stopifnot(nrow(A) == length(PI) &
ncol(A) == length(PI))
stopifnot(nrow(B) == length(PI))
stopifnot(length(state_set) == length(PI))

stopifnot(all(purrrlyr::by_row(as.data.frame(A), sum,
.collate = "cols")[".out"] == 1))
stopifnot(all(purrrlyr::by_row(as.data.frame(B), sum,
.collate = "cols")[".out"] == 1))
stopifnot(sum(PI) == 1)

A <- as.matrix(A); B <- as.matrix(B)
rownames(A) <- state_set; colnames(A) <- state_set
rownames(B) <- state_set; colnames(B) <- observation_set

len_PI <- length(PI)
PI_ <- vector("numeric", length = len_PI)
PI_ <- PI
names(PI_) <- state_set
PI <- PI_

HMMobj <- new("HMM",
A = A,
B = B,
PI = PI,
T_length = as.integer(T_length),
state_set = as.character(state_set),
observation_set = as.character(observation_set))

first_state <- names(PI[sample(len_PI, 1)])
state_series <- vector("character", length = T_length)
state_series[1] <- first_state

observation_series <- vector("character", length = T_length)
observation_1_idx <- fun(B[which(rownames(B) == first_state), ])
observation_series[1] <- names(observation_1_idx)

state <- first_state
for (t in 2:T_length) {
state_idx <- fun(A[state, ])
state_series[t] <- state <- names(state_idx) # 更新state
observation_idx <- fun(B[which(rownames(B) == state), ])
observation_series[t] <- names(observation_idx)
}

HMMobj@state_series <- state_series
HMMobj@observation_series <- observation_series
HMMobj
}

which.random <- function(input) {
idx <- sample(length(input), 1)
idx <- c(idx)
names(idx) <- names(input)[idx]
return(idx)
}

which.random(mtcars[1, ])``````
``````## mpg
##   1``````
``````# page 173 例1
A1 <- matrix(c(0,1,0,0,
0.4,0,0.6,0,
0,0.4,0,0.6,
0,0,0.5,0.5),
nrow = 4,
ncol = 4,
byrow = TRUE)

B1 <- matrix(c(0.5,0.5,
0.3,0.7,
0.6,0.4,
0.8,0.2),
nrow = 4,
ncol = 2,
byrow = TRUE)

PI1 <- rep(0.25, 4)
T_length1 <- 100
state_set <- paste0("盒子", seq(1, 4))
observation_set <- c("红", "白")

HMM_createobser(A = A1, B = B1,
PI = PI1, T_length = T_length1,
state_set = state_set,
observation_set = observation_set)``````
``````## An object of class "HMM"
## Slot "A":
##       盒子1 盒子2 盒子3 盒子4
## 盒子1   0.0   1.0   0.0   0.0
## 盒子2   0.4   0.0   0.6   0.0
## 盒子3   0.0   0.4   0.0   0.6
## 盒子4   0.0   0.0   0.5   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.3 0.7
## 盒子3 0.6 0.4
## 盒子4 0.8 0.2
##
## Slot "PI":
## 盒子1 盒子2 盒子3 盒子4
##  0.25  0.25  0.25  0.25
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3" "盒子4"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 100
##
## Slot "state_series":
##   [1] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##   [9] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [17] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [25] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [33] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [41] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [49] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [57] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [65] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [73] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [81] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [89] "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3" "盒子4" "盒子3"
##  [97] "盒子4" "盒子3" "盒子4" "盒子3"
##
## Slot "observation_series":
##   [1] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [15] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [29] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [43] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [57] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [71] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [85] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [99] "红" "红"``````
``````HMM_createobser(A = A1, B = B1,
PI = PI1, T_length = T_length1,
state_set = state_set,
observation_set = observation_set,
fun = which.min)``````
``````## An object of class "HMM"
## Slot "A":
##       盒子1 盒子2 盒子3 盒子4
## 盒子1   0.0   1.0   0.0   0.0
## 盒子2   0.4   0.0   0.6   0.0
## 盒子3   0.0   0.4   0.0   0.6
## 盒子4   0.0   0.0   0.5   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.3 0.7
## 盒子3 0.6 0.4
## 盒子4 0.8 0.2
##
## Slot "PI":
## 盒子1 盒子2 盒子3 盒子4
##  0.25  0.25  0.25  0.25
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3" "盒子4"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 100
##
## Slot "state_series":
##   [1] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##   [9] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [17] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [25] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [33] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [41] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [49] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [57] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [65] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [73] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [81] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [89] "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2" "盒子2"
##  [97] "盒子2" "盒子2" "盒子2" "盒子2"
##
## Slot "observation_series":
##   [1] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [15] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [29] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [43] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [57] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [71] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [85] "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红" "红"
##  [99] "红" "红"``````
``````HMM_createobser(A = A1, B = B1,
PI = PI1, T_length = T_length1,
state_set = state_set,
observation_set = observation_set,
fun = which.random)``````
``````## An object of class "HMM"
## Slot "A":
##       盒子1 盒子2 盒子3 盒子4
## 盒子1   0.0   1.0   0.0   0.0
## 盒子2   0.4   0.0   0.6   0.0
## 盒子3   0.0   0.4   0.0   0.6
## 盒子4   0.0   0.0   0.5   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.3 0.7
## 盒子3 0.6 0.4
## 盒子4 0.8 0.2
##
## Slot "PI":
## 盒子1 盒子2 盒子3 盒子4
##  0.25  0.25  0.25  0.25
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3" "盒子4"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 100
##
## Slot "state_series":
##   [1] "盒子2" "盒子3" "盒子3" "盒子3" "盒子2" "盒子2" "盒子1" "盒子4"
##   [9] "盒子3" "盒子3" "盒子3" "盒子4" "盒子4" "盒子3" "盒子4" "盒子1"
##  [17] "盒子2" "盒子4" "盒子1" "盒子4" "盒子2" "盒子1" "盒子3" "盒子4"
##  [25] "盒子3" "盒子3" "盒子4" "盒子3" "盒子3" "盒子1" "盒子4" "盒子3"
##  [33] "盒子1" "盒子2" "盒子1" "盒子1" "盒子2" "盒子3" "盒子1" "盒子1"
##  [41] "盒子4" "盒子1" "盒子1" "盒子4" "盒子2" "盒子4" "盒子1" "盒子2"
##  [49] "盒子2" "盒子4" "盒子1" "盒子4" "盒子4" "盒子1" "盒子4" "盒子1"
##  [57] "盒子4" "盒子2" "盒子1" "盒子1" "盒子3" "盒子4" "盒子2" "盒子1"
##  [65] "盒子4" "盒子1" "盒子3" "盒子3" "盒子3" "盒子4" "盒子2" "盒子4"
##  [73] "盒子3" "盒子1" "盒子2" "盒子2" "盒子4" "盒子4" "盒子1" "盒子2"
##  [81] "盒子4" "盒子4" "盒子4" "盒子2" "盒子2" "盒子4" "盒子3" "盒子4"
##  [89] "盒子2" "盒子3" "盒子1" "盒子2" "盒子4" "盒子3" "盒子3" "盒子2"
##  [97] "盒子4" "盒子3" "盒子2" "盒子2"
##
## Slot "observation_series":
##   [1] "白" "红" "白" "白" "白" "白" "红" "红" "白" "白" "白" "红" "白" "白"
##  [15] "白" "白" "红" "白" "白" "红" "红" "红" "红" "白" "红" "白" "红" "白"
##  [29] "红" "白" "白" "白" "红" "红" "白" "红" "红" "红" "白" "红" "红" "红"
##  [43] "红" "白" "红" "白" "白" "白" "白" "白" "红" "白" "红" "白" "白" "白"
##  [57] "红" "红" "白" "红" "红" "白" "红" "白" "白" "红" "白" "红" "白" "红"
##  [71] "红" "红" "红" "红" "白" "白" "红" "白" "白" "红" "白" "红" "白" "白"
##  [85] "白" "白" "红" "红" "白" "白" "红" "白" "红" "红" "红" "白" "红" "白"
##  [99] "红" "白"``````
``````setClass("HMM_Prob_obser", contains = "HMM",
representation(Prob_obser = "numeric",
alpha_matx = "matrix",
beta_matx = "matrix"))

setGeneric(name = "prob_obs",
def = function(objects1, method, ...) standardGeneric("prob_obs"))``````
``## [1] "prob_obs"``
``````#前向/后向/前向后向
setMethod(f = "prob_obs",
signature = c("HMM_Prob_obser", "character"),
definition =
# method = "forward"
# for 循环中用t遍历, 注意这里不能使用t作为参数名, 改用t_
function(objects1, method, t_) {
PI <- objects1@PI
B <- objects1@B
A <- objects1@A
len_PI <- length(PI)
T_length <- objects1@T_length
state_set <- objects1@state_set
observation_series <- objects1@observation_series

forward_prob_a <- matrix(NA, len_PI, T_length)
rownames(forward_prob_a) <- names(PI)
colnames(forward_prob_a) <- paste0("t_", seq(T_length)) # 可有可无

first_obser <- observation_series[1]
for (state in rownames(forward_prob_a)) {
forward_prob_a[state, 1] <- PI[state]*B[state, first_obser]
}

#PI[state]*B[state, first_obser]
#purrr::map(rownames(forward_prob_a),
#.f = `[`(PI, .)*`[`(B, ., first_obser))

for (t in 2:T_length) {
for (state in rownames(forward_prob_a)) {
forward_prob_a[state, t] <- sum(forward_prob_a[, (t - 1)] *
A[, state]) *
B[state, observation_series[t]]
}
}

backward_prob_b <- matrix(NA, len_PI, T_length)
rownames(backward_prob_b) <- names(PI)
colnames(backward_prob_b) <- paste0("t_", seq(T_length)) # 可有可无

backward_prob_b[, T_length] <- 1
for (t in (T_length - 1):1) {
for (state in rownames(backward_prob_b)) {
backward_prob_b[state, t] <- sum(A[state, ] *
B[, observation_series[t + 1]] *
backward_prob_b[, (t + 1)])
}
}
if (method == "forward" && is.character(t_)) {
P_obser_condi <- sum(forward_prob_a[, T_length])
lis <- list(P_obser_condi, forward_prob_a)
return(lis)
} else if (method == "backward" && is.character(t_)) {
P_obser_condi <- sum(PI*B[, observation_series[1]]*backward_prob_b[, 1])
lis <- list(P_obser_condi, backward_prob_b)
return(lis)
} else if (method == "forwardandbackward" && is.numeric(t_)) {
prob_both <- matrix(NA, len_PI, len_PI)
rownames(prob_both) <- colnames(prob_both) <- names(PI)
# 取t = 1,2,3,4,...,(T_length-1)均可
#t <- 1
for (state_a in rownames(forward_prob_a)) {
for (state_b in rownames(backward_prob_b)) {
prob_both[state_a, state_b] <- forward_prob_a[state_a, t] *
A[state_a, state_b] * B[state_b, observation_series[t + 1]] *
backward_prob_b[state_b, (t + 1)]
}
}
P_obser_condi <- sum(prob_both)
lis <- list(P_obser_condi, forward_prob_a, backward_prob_b)
return(lis)
} else {
return("the method you input is wrong!")
}
})``````
``## [1] "prob_obs"``
``````HMM_obs_Prob <- function(A, B, PI, T_length,
state_set, observation_set,
observation_series,
calcumethod, t) {

stopifnot(nrow(A) == length(state_set) &
ncol(A) == length(state_set))
stopifnot(nrow(B) == length(state_set) &
ncol(B) == length(observation_set))
stopifnot(nrow(A) == length(PI) &
ncol(A) == length(PI))
stopifnot(nrow(B) == length(PI))
stopifnot(length(state_set) == length(PI))

stopifnot(all(purrrlyr::by_row(as.data.frame(A), sum,
.collate = "cols")[".out"] == 1))
stopifnot(all(purrrlyr::by_row(as.data.frame(B), sum,
.collate = "cols")[".out"] == 1))
stopifnot(sum(PI) == 1)
stopifnot(all(observation_series %in% observation_set))

A <- as.matrix(A); B <- as.matrix(B)
rownames(A) <- state_set; colnames(A) <- state_set
rownames(B) <- state_set; colnames(B) <- observation_set

len_PI <- length(PI)
PI_ <- vector("numeric", length = len_PI)
PI_ <- PI
names(PI_) <- state_set
PI <- PI_

HMMProb <- new("HMM_Prob_obser",
A = A,
B = B,
PI = PI,
T_length = as.integer(T_length),
state_set = as.character(state_set),
observation_set = as.character(observation_set),
observation_series = as.character(observation_series))

if (calcumethod == "forward") {
lis <- prob_obs(HMMProb, calcumethod, t)
HMMProb@Prob_obser <- lis[[1]]
HMMProb@alpha_matx <- lis[[2]]
HMMProb
}  else if (calcumethod == "backward") {
lis <- prob_obs(HMMProb, calcumethod, t)
HMMProb@Prob_obser <- lis[[1]]
HMMProb@beta_matx <- lis[[2]]
HMMProb
} else {
lis <- prob_obs(HMMProb, calcumethod, t)
HMMProb@Prob_obser <- lis[[1]]
HMMProb@alpha_matx <- lis[[2]]
HMMProb@beta_matx <- lis[[3]]
HMMProb
}
}``````
``````A2 <- matrix(c(0.5,0.2,0.3,
0.3,0.5,0.2,
0.2,0.3,0.5),
nrow = 3, ncol = 3,
byrow = TRUE)

B2 <- matrix(c(0.5,0.5,
0.4,0.6,
0.7,0.3),
nrow = 3, ncol = 2,
byrow = TRUE)

PI2 <- c(0.2, 0.4, 0.4)
T_length2 <- 3
state_set2 <- paste0("盒子", seq(1, 3))
observation_set2 <- c("红", "白")
obser_series <- c("红", "白", "红")

HMMProb1 <- HMM_obs_Prob(A = A2, B = B2, PI = PI2,
T_length = T_length2,
state_set = state_set2,
observation_set = observation_set2,
observation_series = obser_series,
calcumethod = "forward",
t = "none")
HMMProb1``````
``````## An object of class "HMM_Prob_obser"
## Slot "Prob_obser":
## [1] 0.130218
##
## Slot "alpha_matx":
##        t_1    t_2      t_3
## 盒子1 0.10 0.0770 0.041870
## 盒子2 0.16 0.1104 0.035512
## 盒子3 0.28 0.0606 0.052836
##
## Slot "beta_matx":
## <0 x 0 matrix>
##
## Slot "A":
##       盒子1 盒子2 盒子3
## 盒子1   0.5   0.2   0.3
## 盒子2   0.3   0.5   0.2
## 盒子3   0.2   0.3   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.4 0.6
## 盒子3 0.7 0.3
##
## Slot "PI":
## 盒子1 盒子2 盒子3
##   0.2   0.4   0.4
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 3
##
## Slot "state_series":
## character(0)
##
## Slot "observation_series":
## [1] "红" "白" "红"``````
``HMMProb1@Prob_obser``
``## [1] 0.130218``
``HMMProb1@alpha_matx``
``````##        t_1    t_2      t_3
## 盒子1 0.10 0.0770 0.041870
## 盒子2 0.16 0.1104 0.035512
## 盒子3 0.28 0.0606 0.052836``````
``````HMMProb2 <- HMM_obs_Prob(A = A2, B = B2, PI = PI2,
T_length = T_length2,
state_set = state_set2,
observation_set = observation_set2,
observation_series = obser_series,
calcumethod = "backward",
t = "none")
HMMProb2``````
``````## An object of class "HMM_Prob_obser"
## Slot "Prob_obser":
## [1] 0.130218
##
## Slot "alpha_matx":
## <0 x 0 matrix>
##
## Slot "beta_matx":
##          t_1  t_2 t_3
## 盒子1 0.2451 0.54   1
## 盒子2 0.2622 0.49   1
## 盒子3 0.2277 0.57   1
##
## Slot "A":
##       盒子1 盒子2 盒子3
## 盒子1   0.5   0.2   0.3
## 盒子2   0.3   0.5   0.2
## 盒子3   0.2   0.3   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.4 0.6
## 盒子3 0.7 0.3
##
## Slot "PI":
## 盒子1 盒子2 盒子3
##   0.2   0.4   0.4
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 3
##
## Slot "state_series":
## character(0)
##
## Slot "observation_series":
## [1] "红" "白" "红"``````
``HMMProb2@Prob_obser``
``## [1] 0.130218``
``HMMProb2@beta_matx``
``````##          t_1  t_2 t_3
## 盒子1 0.2451 0.54   1
## 盒子2 0.2622 0.49   1
## 盒子3 0.2277 0.57   1``````
``dim(HMMProb2@alpha_matx)``
``## [1] 0 0``
``````#stopifnot(dim(HMMProb2@alpha_matx) != c(0,0))
stopifnot(dim(HMMProb2@beta_matx)  != c(0,0))

HMMProb3 <- HMM_obs_Prob(A = A2, B = B2, PI = PI2,
T_length = T_length2,
state_set = state_set2,
observation_set = observation_set2,
observation_series = obser_series,
calcumethod = "forwardandbackward",
t = 2)
HMMProb3``````
``````## An object of class "HMM_Prob_obser"
## Slot "Prob_obser":
## [1] 0.130218
##
## Slot "alpha_matx":
##        t_1    t_2      t_3
## 盒子1 0.10 0.0770 0.041870
## 盒子2 0.16 0.1104 0.035512
## 盒子3 0.28 0.0606 0.052836
##
## Slot "beta_matx":
##          t_1  t_2 t_3
## 盒子1 0.2451 0.54   1
## 盒子2 0.2622 0.49   1
## 盒子3 0.2277 0.57   1
##
## Slot "A":
##       盒子1 盒子2 盒子3
## 盒子1   0.5   0.2   0.3
## 盒子2   0.3   0.5   0.2
## 盒子3   0.2   0.3   0.5
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.4 0.6
## 盒子3 0.7 0.3
##
## Slot "PI":
## 盒子1 盒子2 盒子3
##   0.2   0.4   0.4
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 3
##
## Slot "state_series":
## character(0)
##
## Slot "observation_series":
## [1] "红" "白" "红"``````
``HMMProb3@Prob_obser``
``## [1] 0.130218``
``HMMProb3@alpha_matx``
``````##        t_1    t_2      t_3
## 盒子1 0.10 0.0770 0.041870
## 盒子2 0.16 0.1104 0.035512
## 盒子3 0.28 0.0606 0.052836``````
``HMMProb3@beta_matx``
``````##          t_1  t_2 t_3
## 盒子1 0.2451 0.54   1
## 盒子2 0.2622 0.49   1
## 盒子3 0.2277 0.57   1``````
``````setGeneric(name = "gammafunc",
def = function(objects1, t, state, ...) standardGeneric("gammafunc"))``````
``## [1] "gammafunc"``
``````#计算gamma
setMethod(f = "gammafunc",
signature = c("HMM_Prob_obser", "numeric", "character"),
definition =
function(objects1, t, state) {
stopifnot(dim(objects1@alpha_matx) != c(0,0))
stopifnot(dim(objects1@beta_matx)  != c(0,0))

gammaval <-
objects1@alpha_matx[state, t]*objects1@beta_matx[state, t] /
sum(objects1@alpha_matx[, t]*objects1@beta_matx[, t])
return(gammaval)
})``````
``## [1] "gammafunc"``
``gammafunc(HMMProb3, 2, "盒子1")``
``## [1] 0.3193107``
``````#公式10.40分子，但还未对t累加求和
setGeneric("gammafunc_obs_k",
function(objects1, t, state, obs, ...) standardGeneric("gammafunc_obs_k"))``````
``## [1] "gammafunc_obs_k"``
``````setMethod(f = "gammafunc_obs_k",
signature = c("HMM_Prob_obser", "numeric",
"character", "character"),
definition =
function(objects1, t, state, obs) {
stopifnot(dim(objects1@alpha_matx) != c(0,0))
stopifnot(dim(objects1@beta_matx)  != c(0,0))

state_set <- objects1@state_set
observation_set <- objects1@observation_set
state_obs_all <- matrix(NA, length(state_set), length(observation_set))
rownames(state_obs_all) <- state_set
colnames(state_obs_all) <- observation_set

for (j in state_set) {
for (k in observation_set) {
state_obs_all[j, k] <- objects1@alpha_matx[j, t] *
objects1@beta_matx[j, t] * objects1@B[j, k]
}
}
denominator <- sum(state_obs_all)
gammaval_obs <- objects1@alpha_matx[state, t] *
objects1@beta_matx[state, t] *
objects1@B[state, obs] / denominator
return(gammaval_obs)
})``````
``## [1] "gammafunc_obs_k"``
``gammafunc_obs_k(HMMProb3, 2, "盒子1", "红")``
``## [1] 0.1596553``
``````setGeneric(name = "ksifunc",
def = function(objects1, t, state_t, state_tplus1, ...) standardGeneric("ksifunc"))``````
``## [1] "ksifunc"``
``````#计算ksi
setMethod(f = "ksifunc",
signature = c("HMM_Prob_obser", "numeric", "character", "character"),
definition =
function(objects1, t, state_t, state_tplus1) {
stopifnot(dim(objects1@alpha_matx) != c(0, 0))
stopifnot(dim(objects1@beta_matx)  != c(0, 0))

obs <- objects1@observation_series[t + 1]
lis <- prob_obs(objects1, "forwardandbackward", t)

ksival <-
objects1@alpha_matx[state_t, t] * objects1@A[state_t, state_tplus1] *
objects1@B[state_tplus1, obs] *
objects1@beta_matx[state_tplus1, t + 1] / lis[[1]]
return(ksival)
})``````
``## [1] "ksifunc"``
``ksifunc(HMMProb3, 1, "盒子3", "盒子3")``
``## [1] 0.1838456``
``ksifunc(HMMProb3, 2, "盒子3", "盒子3")``
``## [1] 0.1628807``
``````setGeneric("expectstate",
function(objects1, state, ...) standardGeneric("expectstate"))``````
``## [1] "expectstate"``
``````#3(1) page180
setMethod("expectstate",
c("HMM_Prob_obser", "character"),
function(objects1, state) {
T_len <- objects1@T_length

val <- 0
for (t in 1:T_len) {
val_ <- gammafunc(objects1, t, state)
val <- val + val_
}
return(val)
})``````
``## [1] "expectstate"``
``expectstate(HMMProb3, "盒子3")``
``## [1] 1.160623``
``````#公式10.40分子，对t累加求和
setGeneric("expectstate_obs_k",
function(objects1, state, obs, ...) standardGeneric("expectstate_obs_k"))``````
``## [1] "expectstate_obs_k"``
``````setMethod("expectstate_obs_k",
c("HMM_Prob_obser", "character", "character"),
function(objects1, state, obs) {
T_len <- objects1@T_length

val <- 0
for (t in 1:T_len) {
val_ <- gammafunc_obs_k(objects1, t, state, obs)
val <- val + val_
}
return(val)
})``````
``## [1] "expectstate_obs_k"``
``expectstate_obs_k(HMMProb3, "盒子3", "红")``
``## [1] 0.8124361``
``````setGeneric("expect_fromstate",
function(objects1, state, ...) standardGeneric("expect_fromstate"))``````
``## [1] "expect_fromstate"``
``````#3(2) page180
setMethod("expect_fromstate",
c("HMM_Prob_obser", "character"),
function(objects1, state) {
T_len <- objects1@T_length

val <- 0
for (t in 1:(T_len-1)) {
val_ <- gammafunc(objects1, t, state)
val <- val + val_
}
return(val)
})``````
``## [1] "expect_fromstate"``
``expect_fromstate(HMMProb3, "盒子3")``
``## [1] 0.7548726``
``````setGeneric("expect_state_itoj",
function(objects1, statei, statej, ...) standardGeneric("expect_state_itoj"))``````
``## [1] "expect_state_itoj"``
``````#3(3) page180
setMethod("expect_state_itoj",
c("HMM_Prob_obser", "character", "character"),
function(objects1, statei, statej) {
T_len <- objects1@T_length

val <- 0
for (t in 1:(T_len - 1)) {
val_ <- ksifunc(objects1, t, statei, statej)
val <- val + val_
}
return(val)
})``````
``## [1] "expect_state_itoj"``
``expect_state_itoj(HMMProb3, "盒子3", "盒子1")``
``## [1] 0.1626503``
``````# 维特比算法封装进viterbi_bestpath()函数
setClass("newHMM", representation(A = "matrix",
B = "matrix",
PI = "numeric",
state_set = "character",
observation_set = "character",
T_length = "numeric",
state_series = "character",
observation_series = "character"))

viterbi_bestpath <- function(input_A, input_B, input_PI,
input_state_set, input_observation_set,
input_T_length, input_observation_series) {

input_A <- as.matrix(input_A); input_B <- as.matrix(input_B)
rownames(input_A) <- input_state_set
colnames(input_A) <- input_state_set
rownames(input_B) <- input_state_set
colnames(input_B) <- input_observation_set

len_PI <- length(input_PI)
PI_ <- vector("numeric", length = len_PI)
PI_ <- input_PI
names(PI_) <- input_state_set
input_PI <- PI_

newHMM_obj <- new("newHMM",
A = input_A,
B = input_B,
PI = input_PI,
state_set = input_state_set,
observation_set = input_observation_set,
T_length = input_T_length,
observation_series = input_observation_series)

delta_matx <- matrix(NA, length(input_state_set), input_T_length)
rownames(delta_matx) <- input_state_set
colnames(delta_matx) <- paste0("t_", seq(input_T_length))
#delta_matx

psai_matx <- matrix(NA, length(input_state_set), input_T_length)
rownames(psai_matx) <- input_state_set
colnames(psai_matx) <- paste0("t_", seq(input_T_length))
#psai_matx

for (i in input_state_set) {
delta_matx[i, 1] <- newHMM_obj@PI[i] *
newHMM_obj@B[i, input_observation_series[1]]
}
for (i in input_state_set) {
psai_matx[i, 1] <- 0
}
for (t in 2:input_T_length) {
for (i in input_state_set) {
delta_matx[i, t] <- max(delta_matx[, t - 1] *
newHMM_obj@A[, i]) *
newHMM_obj@B[i, input_observation_series[t]]
psai_matx[i, t] <- input_state_set[which.max(delta_matx[, t - 1] *
newHMM_obj@A[, i])]
}
}

Pstar <- max(delta_matx[, input_T_length])

best_path <- vector("character", length = input_T_length)
best_path[input_T_length] <- input_state_set[which.max(delta_matx[, input_T_length])]

for (t in (input_T_length - 1):1) {
best_path[t] <- psai_matx[best_path[t + 1], t + 1]
}
newHMM_obj@state_series <- best_path
return(newHMM_obj)
}``````
``````# 例10.3
A0 <- matrix(c(0.5,0.2,0.3,
0.3,0.5,0.2,
0.2,0.3,0.5),
nrow = 3, ncol = 3,
byrow = TRUE)

B0 <- matrix(c(0.5,0.5,
0.4,0.6,
0.7,0.3),
nrow = 3, ncol = 2,
byrow = TRUE)

PI0 <- c(0.2, 0.4, 0.4)
obs <- c("红", "白", "红")
T_leng <- length(obs)

states_set <- paste0("盒子", seq(1, 3))
obs_set <- c("红", "白")

viterbi_bestpath(A0, B0, PI0,
states_set, obs_set,
T_leng, obs) -> HMM_obj
HMM_obj@state_series``````
``## [1] "盒子3" "盒子3" "盒子3"``
``````# page189 10.3练习题
A0 <- matrix(c(0.5,0.2,0.3,
0.3,0.5,0.2,
0.2,0.3,0.5),
nrow = 3, ncol = 3,
byrow = TRUE)

B0 <- matrix(c(0.5,0.5,
0.4,0.6,
0.7,0.3),
nrow = 3, ncol = 2,
byrow = TRUE)

PI0 <- c(0.2, 0.4, 0.4)
obs <- c("红", "白", "红", "白")
T_leng <- length(obs)

states_set <- paste0("盒子", seq(1, nrow(A0)))
obs_set <- c("红", "白")

viterbi_bestpath(A0, B0, PI0,
states_set, obs_set,
T_leng, obs) -> HMM_obj1
HMM_obj1@state_series``````
``## [1] "盒子3" "盒子2" "盒子2" "盒子2"``
``````# page188 10.1练习题
HMM_obs_Prob(A0, B0, PI0, T_leng,
states_set, obs_set,
obs, "backward", "none") -> HMM_obs_Prob1
HMM_obs_Prob1@Prob_obser``````
``## [1] 0.0600908``
``````# page189 10.2练习题
A0 <- matrix(c(0.5,0.1,0.4,
0.3,0.5,0.2,
0.2,0.2,0.6),
nrow = 3, ncol = 3,
byrow = TRUE)

B0 <- matrix(c(0.5,0.5,
0.4,0.6,
0.7,0.3),
nrow = 3, ncol = 2,
byrow = TRUE)

PI0 <- c(0.2, 0.3, 0.5)
obs <- c("红", "白", "红", "红", "白", "红", "白", "白")
T_leng <- length(obs)

states_set <- paste0("盒子", seq(1, nrow(A0)))
obs_set <- c("红", "白")

HMM_obs_Prob(A0, B0, PI0, T_leng,
states_set, obs_set,
obs, "forwardandbackward", 4) -> HMM_obs_Prob2
HMM_obs_Prob2``````
``````## An object of class "HMM_Prob_obser"
## Slot "Prob_obser":
## [1] 0.003476709
##
## Slot "alpha_matx":
##        t_1    t_2      t_3        t_4        t_5         t_6         t_7
## 盒子1 0.10 0.0780 0.040320 0.02086680 0.01143210 0.005496686 0.002810565
## 盒子2 0.12 0.0840 0.026496 0.01236192 0.01019392 0.003383726 0.002459520
## 盒子3 0.35 0.0822 0.068124 0.04361112 0.01109573 0.009288344 0.002534528
##               t_8
## 盒子1 0.001325022
## 盒子2 0.001210633
## 盒子3 0.000941054
##
## Slot "beta_matx":
##               t_1        t_2        t_3        t_4      t_5    t_6  t_7
## 盒子1 0.006325692 0.01482964 0.02556442 0.04586531 0.105521 0.1861 0.43
## 盒子2 0.006847065 0.01227214 0.02343448 0.05280909 0.100883 0.2415 0.51
## 盒子3 0.005778550 0.01568294 0.02678985 0.04280618 0.111934 0.1762 0.40
##       t_8
## 盒子1   1
## 盒子2   1
## 盒子3   1
##
## Slot "A":
##       盒子1 盒子2 盒子3
## 盒子1   0.5   0.1   0.4
## 盒子2   0.3   0.5   0.2
## 盒子3   0.2   0.2   0.6
##
## Slot "B":
##        红  白
## 盒子1 0.5 0.5
## 盒子2 0.4 0.6
## 盒子3 0.7 0.3
##
## Slot "PI":
## 盒子1 盒子2 盒子3
##   0.2   0.3   0.5
##
## Slot "state_set":
## [1] "盒子1" "盒子2" "盒子3"
##
## Slot "observation_set":
## [1] "红" "白"
##
## Slot "T_length":
## [1] 8
##
## Slot "state_series":
## character(0)
##
## Slot "observation_series":
## [1] "红" "白" "红" "红" "白" "红" "白" "白"``````
``gammafunc(HMM_obs_Prob2, 4, states_set[3])``
``## [1] 0.5369518``
``````# 有待改进
# 10.4 Baum-Welch算法
A0 <- matrix(c(0.5,0.2,0.3,
0.3,0.5,0.2,
0.2,0.3,0.5),
nrow = 3, ncol = 3,
byrow = TRUE)

B0 <- matrix(c(0.5,0.5,
0.4,0.6,
0.7,0.3),
nrow = 3, ncol = 2,
byrow = TRUE)

PI0 <- c(0.2, 0.4, 0.4)
obs <- c("红", "白", "红")
T_leng <- length(obs)

states_set <- paste0("盒子", seq(1, 3))
obs_set <- c("红", "白")

HMMProb4 <- HMM_obs_Prob(A = A0, B = B0, PI = PI0,
T_length = T_leng,
state_set = states_set,
observation_set = obs_set,
observation_series = obs,
calcumethod = "forwardandbackward",
t = 1) # t 取 1 至 T_leng均可

AA <- matrix(NA, length(states_set), length(states_set))
rownames(AA) <- colnames(AA) <- states_set; AA
BB <- matrix(NA, length(states_set), length(obs_set))
rownames(BB) <- states_set; colnames(BB) <- obs_set; BB
PIPI <- vector("numeric", length = length(states_set))
names(PIPI) <- states_set; PIPI

diff_statemtx <- diff_emissmtx <- Inf
eps <- 0.0001
step <- 1
while (diff_statemtx > eps || diff_emissmtx > eps) {

if (step <= 100) {
for (i in states_set) {
for (j in states_set) {
AA[i, j] <- expect_state_itoj(HMMProb4, i, j) /
expect_fromstate(HMMProb4, i)
}
}
diff_statemtx <- norm(AA - HMMProb4@A, type = "2")

for (j in states_set) {
for (k in obs_set) {
BB[j, k] <- expectstate_obs_k(HMMProb4, j, k) /
expectstate(HMMProb4, j)
}
}
diff_emissmtx <- norm(BB - HMMProb4@B, type = "2")

for (i in states_set) {
PIPI[i] <- gammafunc(HMMProb4, 1, i)
}

HMMProb4@A <- AA
HMMProb4@B <- BB
HMMProb4@PI <- PIPI

print(AA)
cat("\n")
print(BB)
cat("\n")
print(PIPI)
cat("\n")
} else {
print(AA)
cat("\n")
print(BB)
cat("\n")
print(PIPI)
cat("\n")
}
step <- step + 1
}``````