# 马科夫模型和马科夫链

### 战立侃 · 2019-08-13

suppressMessages(library(markovchain))
suppressMessages(library(igraph)) # getSlots("markovchain")
suppressMessages(library(expm))   #  %^%
suppressMessages(require(matlab)) # canonicForm
library(ggplot2)

X0 <- c(1, 0)
TI2 <- matrix(c(0.8, 0.2, 1, 0), nrow = 2, byrow = TRUE)
(X1 <- X0 %*% TI2)
##      [,1] [,2]
## [1,]  0.8  0.2
(X2 <- X1 %*% TI2)
##      [,1] [,2]
## [1,] 0.84 0.16

round(TI2 %^% 4, 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
round(TI2 %^% 100, 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17

# 非周期不可化简的(irreducible and aperiodic)
(XI2 <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
round(XI2 %*% (TI2 %^% 20), 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
round(XI2 %*% (TI2 %^% 21), 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
# 周期性的(periodic)
(TP2 <- matrix(c(0, 1, 1, 0), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    0    1
## [2,]    1    0
XI2 %*% (TP2 %^% 20)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TP2 %^% 21)
##      [,1] [,2]
## [1,]    0    1
## [2,]    1    0
# 可化简的(reducible)
(TR2 <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TR2 %^% 20)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TR2 %^% 21)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

# 三种状态、非周期不可化简的(irreducible and aperiodic)
(XI3 <- matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1),
nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
(TI3 <- matrix(c(0.8, 0.2, 0, 0.8, 0, 0.2, 0.1, 0, 0.9),
nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]  0.8  0.2  0.0
## [2,]  0.8  0.0  0.2
## [3,]  0.1  0.0  0.9
round(XI3 %*% (TI3 %^% 50), 3)
##       [,1]  [,2] [,3]
## [1,] 0.625 0.125 0.25
## [2,] 0.625 0.125 0.25
## [3,] 0.625 0.125 0.25

# 三种状态、可化简的
(TR3 <- matrix(c(0.8, 0.2, 0, 0.8, 0, 0.2, 0, 0, 1),
nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]  0.8  0.2  0.0
## [2,]  0.8  0.0  0.2
## [3,]  0.0  0.0  1.0
round(XI3 %*% (TR3 %^% 63), 6)
##          [,1]     [,2]     [,3]
## [1,] 0.094597 0.019592 0.885811
## [2,] 0.078367 0.016230 0.905403
## [3,] 0.000000 0.000000 1.000000
round(XI3 %*% (TR3 %^% 153), 2)
##      [,1] [,2] [,3]
## [1,]    0    0    1
## [2,]    0    0    1
## [3,]    0    0    1

MTR3 <- as(TR3, "markovchain")
States <- c("G", "M", "A")
MTR3@states <- States
attr(MTR3@transitionMatrix, "dimnames") <- list(States, States)
ExtMatr <- extractMatrices(MTR3)
# names(ExtMatr)
FN <- ExtMatr[["N"]]
FN
##    G M
## G 25 5
## M 20 5
FN %*% c(1, 1)
##   [,1]
## G   30
## M   25

FM <- function(T0){
LL <- sapply(1:2500,
function(i) grep("A", rmarkovchain(60, MTR3, t0 = T0, include.t0 = TRUE))) # markovchainSequence
LL <- LL[lengths(LL) > 0L]
#LL[lengths(LL) == 0L] <- 60
LLS <- sapply(LL, min)
mean(LLS)
}
FM("G")
## [1] 22.57484
FM("M")
## [1] 18.30045
FM("G") - FM("M")
## [1] 3.792273
Simulate <- function(start) {
LL <- lapply(1:2500, function(i) rmarkovchain(60, MTR3, t0 = start))
LL <- do.call(rbind, LL)
LL <- gsub("A", 1, LL)
LL <- gsub("G|M", 0, LL)
mode(LL) <- "numeric"
Probs <- colMeans(LL) # colSums(LL) / 2500
SS <- data.frame(
Prob  = Probs, Steps = 1:60,
Start = start, stringsAsFactors = FALSE)
return(SS)
}

DS1 <- Simulate("G")
mean(DS1 $Prob * DS1$ Steps)
## [1] 21.73269
DS2 <- Simulate("M")
mean(DS2 $Prob * DS2$ Steps)
## [1] 23.40962
DS2 $Steps <- DS2$ Steps + 5
mean(DS2 $Prob * DS2$ Steps)
## [1] 26.66659
DS  <- rbind(DS1, DS2)

MTR3@transitionMatrix[2, ] <- c(0.6, 0, 0.4)
ExtMatr <- extractMatrices(MTR3)
# names(ExtMatr)
FN <- ExtMatr[["N"]]
FN
##      G   M
## G 12.5 2.5
## M  7.5 2.5
FN %*% c(1, 1)
##   [,1]
## G   15
## M   10`

Grewal, J. K., Krzywinski, M., & Altman, N. (2019). Points of significance: Markov models—Markov chains. Nature Methods, 16(8), 661. doi:10.1038/s41592-019-0463-2