# First Passage Time Moments,
# Howard, Dynamic Probabilistic Systems, Vol. I, Markov Models, 1971, p.305
# (5.1.22), (5.2.63), (5.2.73)
# R-code by PCM 2015/05/07
# limiting multistep transition probability matrix phi (5.2.71), p.305
my_phi <- function(P)
# P = one-step transition matrix
{
my_pi <- function(P)
# R-code adapted from Matloff, The art of R programming, 2011, p.200f
# P = one-step transition matrix
{n <- nrow(P)
I_minus_P_transp <- diag(n) - t(P)
I_minus_P_transp[n,] <- rep(1,n)
rhs <- c(rep(0,n-1),1)
pi_vector <- solve(I_minus_P_transp,rhs)
return(pi_vector)
}
pi_vec <- my_pi(P)
n <- length(pi_vec)
phi <- matrix(pi_vec,nrow=n,ncol=n,byrow=TRUE)
return(phi)
}
# transient sum matrix Tg1 (5.1.122), (5.2.73)
my_Tg1 <- function(P, PHI)
# P = one-step transition matrix
# PHI = limiting multistep transition probability matrix
{
n <- nrow(P)
I_min_P <- diag(1,nrow=n,ncol=n) - P
ifelse(n <= 5, {print("[I - P]"); print(I_min_P)}, FALSE)
I_min_P_plus_PHI <- I_min_P + PHI
ifelse(n <= 5, {print("[I - P + PHI]"); print(I_min_P_plus_PHI)}, FALSE)
I_min_P_plus_PHI_inv <- solve(I_min_P_plus_PHI)
ifelse(n <= 5, {print("[I - P + PHI]^-1"); print(I_min_P_plus_PHI_inv)}, FALSE)
ifelse(n <= 5, {print("Check1:[I-P+PHI][I-P+PHI]^-1"); print(I_min_P_plus_PHI %*% I_min_P_plus_PHI_inv)}, FALSE)
ifelse(n <= 5, {print("Check2:[I-P+PHI]^-1[I-P+PHI]"); print(I_min_P_plus_PHI_inv %*% I_min_P_plus_PHI)}, FALSE)
ifelse(n <= 5, {print("Transient Sum Matrix Tg1 (5.2.73)"); print(I_min_P_plus_PHI_inv - PHI)}, FALSE)
return(I_min_P_plus_PHI_inv - PHI)
}
# Mean first passage time (MFPT) matrix Theta_bar (5.2.63) p.302
my_theta_bar <- function(PHI, Tg1)
{
n <- nrow(Tg1)
Tg1_box_I <- diag(diag(Tg1), nrow=n, ncol=n)
ifelse(n <= 5, {print("Tg1_box_I"); print(Tg1_box_I)}, FALSE)
U <- matrix(rep(1,n*n), nrow=n, ncol=n)
ifelse(n <= 5, {print("U"); print(U)} , FALSE)
U_Tg1_box_I <- U %*% Tg1_box_I
ifelse(n <= 5, {print("U_Tg1_box_I"); print(U_Tg1_box_I)}, FALSE)
first_theta_bar <- diag(1, nrow=n, ncol=n) + U_Tg1_box_I -Tg1
ifelse(n <= 5, {print("first_theta_bar"); print(first_theta_bar)}, FALSE)
PHI_box_I_inv <- diag(1.0/diag(PHI), nrow=n, ncol=n)
ifelse(n <= 5, {print("PHI_box_I_inv"); print(PHI_box_I_inv)}, FALSE)
theta_bar <- first_theta_bar %*% PHI_box_I_inv
cat("theta_bar = Mean First Passage Time (MFPT) matrix (5.2.63)\n")
return(theta_bar)
}
P <- matrix(c(0.3,0.2,0.5,0.1,0.8,0.1,0.4,0.4,0.2),nrow=3,byrow=TRUE) # (5.2.70), p.305
P
## [,1] [,2] [,3]
## [1,] 0.3 0.2 0.5
## [2,] 0.1 0.8 0.1
## [3,] 0.4 0.4 0.2
PHI <- my_phi(P)
ifelse(nrow(P) <= 5, my_Tg1(P, PHI), FALSE)
## [1] "[I - P]"
## [,1] [,2] [,3]
## [1,] 0.7 -0.2 -0.5
## [2,] -0.1 0.2 -0.1
## [3,] -0.4 -0.4 0.8
## [1] "[I - P + PHI]"
## [,1] [,2] [,3]
## [1,] 0.9 0.4 -0.3
## [2,] 0.1 0.8 0.1
## [3,] -0.2 0.2 1.0
## [1] "[I - P + PHI]^-1"
## [,1] [,2] [,3]
## [1,] 1.3 -0.7666667 0.4666667
## [2,] -0.2 1.4000000 -0.2000000
## [3,] 0.3 -0.4333333 1.1333333
## [1] "Check1:[I-P+PHI][I-P+PHI]^-1"
## [,1] [,2] [,3]
## [1,] 1.000000e+00 5.551115e-17 0.000000e+00
## [2,] 0.000000e+00 1.000000e+00 2.775558e-17
## [3,] -1.110223e-16 -5.551115e-17 1.000000e+00
## [1] "Check2:[I-P+PHI]^-1[I-P+PHI]"
## [,1] [,2] [,3]
## [1,] 1.000000e+00 2.081668e-16 0
## [2,] -1.387779e-17 1.000000e+00 0
## [3,] -5.551115e-17 0.000000e+00 1
## [1] "Transient Sum Matrix Tg1 (5.2.73)"
## [,1] [,2] [,3]
## [1,] 1.1 -1.366667 0.2666667
## [2,] -0.4 0.800000 -0.4000000
## [3,] 0.1 -1.033333 0.9333333
## [1] 1.1
my_theta_bar(PHI, my_Tg1(P,my_phi(P)))
## [1] "[I - P]"
## [,1] [,2] [,3]
## [1,] 0.7 -0.2 -0.5
## [2,] -0.1 0.2 -0.1
## [3,] -0.4 -0.4 0.8
## [1] "[I - P + PHI]"
## [,1] [,2] [,3]
## [1,] 0.9 0.4 -0.3
## [2,] 0.1 0.8 0.1
## [3,] -0.2 0.2 1.0
## [1] "[I - P + PHI]^-1"
## [,1] [,2] [,3]
## [1,] 1.3 -0.7666667 0.4666667
## [2,] -0.2 1.4000000 -0.2000000
## [3,] 0.3 -0.4333333 1.1333333
## [1] "Check1:[I-P+PHI][I-P+PHI]^-1"
## [,1] [,2] [,3]
## [1,] 1.000000e+00 5.551115e-17 0.000000e+00
## [2,] 0.000000e+00 1.000000e+00 2.775558e-17
## [3,] -1.110223e-16 -5.551115e-17 1.000000e+00
## [1] "Check2:[I-P+PHI]^-1[I-P+PHI]"
## [,1] [,2] [,3]
## [1,] 1.000000e+00 2.081668e-16 0
## [2,] -1.387779e-17 1.000000e+00 0
## [3,] -5.551115e-17 0.000000e+00 1
## [1] "Transient Sum Matrix Tg1 (5.2.73)"
## [,1] [,2] [,3]
## [1,] 1.1 -1.366667 0.2666667
## [2,] -0.4 0.800000 -0.4000000
## [3,] 0.1 -1.033333 0.9333333
## [1] "Tg1_box_I"
## [,1] [,2] [,3]
## [1,] 1.1 0.0 0.0000000
## [2,] 0.0 0.8 0.0000000
## [3,] 0.0 0.0 0.9333333
## [1] "U"
## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 1 1
## [3,] 1 1 1
## [1] "U_Tg1_box_I"
## [,1] [,2] [,3]
## [1,] 1.1 0.8 0.9333333
## [2,] 1.1 0.8 0.9333333
## [3,] 1.1 0.8 0.9333333
## [1] "first_theta_bar"
## [,1] [,2] [,3]
## [1,] 1.0 2.166667 0.6666667
## [2,] 1.5 1.000000 1.3333333
## [3,] 1.0 1.833333 1.0000000
## [1] "PHI_box_I_inv"
## [,1] [,2] [,3]
## [1,] 5 0.000000 0
## [2,] 0 1.666667 0
## [3,] 0 0.000000 5
## theta_bar = Mean First Passage Time (MFPT) matrix (5.2.63)
## [,1] [,2] [,3]
## [1,] 5.0 3.611111 3.333333
## [2,] 7.5 1.666667 6.666667
## [3,] 5.0 3.055556 5.000000
# MacKay's Example and Rule of Thumb
# P = transition matrix P(State=j at t+1 | State=i at t) = P[i, j]
P <- matrix(c(
.50,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 01
.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 02
.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 03
.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 04
.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 05
.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 06
.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 07
.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 08
.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 09
.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 10
.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00,.00, # row 11
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00,.00, # row 12
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00,.00, # row 13
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00,.00, # row 14
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00,.00, # row 15
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00,.00, # row 16
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00,.00, # row 17
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00,.00, # row 18
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50,.00, # row 19
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.00,.50, # row 20
.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.00,.50,.50) # row 21
, nrow=21, ncol=21, byrow=TRUE)
#
PHI <- my_phi(P)
ifelse(nrow(P) <= 5, my_Tg1(P, PHI), FALSE)
## [1] FALSE
theta_bar <- my_theta_bar(PHI, my_Tg1(P, PHI))
## theta_bar = Mean First Passage Time (MFPT) matrix (5.2.63)
theta_bar
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 21 2 6 12 20 30 42 56 72 90 110 132 156
## [2,] 40 21 4 10 18 28 40 54 70 88 108 130 154
## [3,] 78 38 21 6 14 24 36 50 66 84 104 126 150
## [4,] 114 74 36 21 8 18 30 44 60 78 98 120 144
## [5,] 148 108 70 34 21 10 22 36 52 70 90 112 136
## [6,] 180 140 102 66 32 21 12 26 42 60 80 102 126
## [7,] 210 170 132 96 62 30 21 14 30 48 68 90 114
## [8,] 238 198 160 124 90 58 28 21 16 34 54 76 100
## [9,] 264 224 186 150 116 84 54 26 21 18 38 60 84
## [10,] 288 248 210 174 140 108 78 50 24 21 20 42 66
## [11,] 310 270 232 196 162 130 100 72 46 22 21 22 46
## [12,] 330 290 252 216 182 150 120 92 66 42 20 21 24
## [13,] 348 308 270 234 200 168 138 110 84 60 38 18 21
## [14,] 364 324 286 250 216 184 154 126 100 76 54 34 16
## [15,] 378 338 300 264 230 198 168 140 114 90 68 48 30
## [16,] 390 350 312 276 242 210 180 152 126 102 80 60 42
## [17,] 400 360 322 286 252 220 190 162 136 112 90 70 52
## [18,] 408 368 330 294 260 228 198 170 144 120 98 78 60
## [19,] 414 374 336 300 266 234 204 176 150 126 104 84 66
## [20,] 418 378 340 304 270 238 208 180 154 130 108 88 70
## [21,] 420 380 342 306 272 240 210 182 156 132 110 90 72
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 182 210 240 272 306 342 380 420
## [2,] 180 208 238 270 304 340 378 418
## [3,] 176 204 234 266 300 336 374 414
## [4,] 170 198 228 260 294 330 368 408
## [5,] 162 190 220 252 286 322 360 400
## [6,] 152 180 210 242 276 312 350 390
## [7,] 140 168 198 230 264 300 338 378
## [8,] 126 154 184 216 250 286 324 364
## [9,] 110 138 168 200 234 270 308 348
## [10,] 92 120 150 182 216 252 290 330
## [11,] 72 100 130 162 196 232 270 310
## [12,] 50 78 108 140 174 210 248 288
## [13,] 26 54 84 116 150 186 224 264
## [14,] 21 28 58 90 124 160 198 238
## [15,] 14 21 30 62 96 132 170 210
## [16,] 26 12 21 32 66 102 140 180
## [17,] 36 22 10 21 34 70 108 148
## [18,] 44 30 18 8 21 36 74 114
## [19,] 50 36 24 14 6 21 38 78
## [20,] 54 40 28 18 10 4 21 40
## [21,] 56 42 30 20 12 6 2 21
cat("MFPT(x(0)=10 -> x(MFPT)=20) = \n")
## MFPT(x(0)=10 -> x(MFPT)=20) =
print(theta_bar[11,21])
## [1] 310
cat("MFPT(x(0)=0 -> x(MFPT)=20) = \n")
## MFPT(x(0)=0 -> x(MFPT)=20) =
print(theta_bar[1,21])
## [1] 420