Chapter 10 隐马尔可夫模型

chap10_1

chap10_1

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
chap10_2

chap10_2

# 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] "红" "白"
chap10_3

chap10_3

chap10_4

chap10_4

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
  }
}
chap10_5

chap10_5

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
chap10_6

chap10_6

# 维特比算法封装进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)
}
chap10_7

chap10_7

# 例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"
chap10_8

chap10_8

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