1. In Example 11.8, make states 0 and 4 into absorbing states. Find the fundamental matrix N, and also Nc and NR, for the resulting absorbing chain. Interpret the results.

Example 11.8 (Ehrenfest Model) The following is a special case of a model, called the Ehrenfest model,3 that has been used to explain diā†µusion of gases. The general model will be discussed in detail in Section 11.5. We have two urns that, between them, contain four balls. At each step, one of the four balls is chosen at random and moved from the urn that it is in into the other urn. We choose, as states, the number of balls in the first urn. The transition matrix is then

probability_left = 1/3
probability_right = 2/3
s = c(probability_left, probability_right)
ehrenfest = matrix(c(0, .25, 0, 0, 0, 1, 0, .5, 0, 0, 0, .75, 0, .75, 0, 0, 0, .5, 0, 1, 0, 0, 0, .25, 0), nrow=5, ncol=5)
ehrenfest
##      [,1] [,2] [,3] [,4] [,5]
## [1,] 0.00  1.0 0.00  0.0 0.00
## [2,] 0.25  0.0 0.75  0.0 0.00
## [3,] 0.00  0.5 0.00  0.5 0.00
## [4,] 0.00  0.0 0.75  0.0 0.25
## [5,] 0.00  0.0 0.00  1.0 0.00

A state si of a Markov chain is called absorbing if it is impossible to leave it (i.e., pii = 1). States 0 and 4 are already absorbing states and thus do not need to be transformed. We will alter this matrix so it will be in canonical form by bringing columns 2 and 4 (labeling 1-5) to the right and brining rows 1 and 5 (labeling 1-5, but are states 0 and 4) to the bottom. We get the resultant matrix.

canonical = matrix(c(.25, 0, 0, 0, 0, .75, 0, .75, 0, 0, 0, 0, .25, 0, 0, 0, .5, 0, 1, 0, 0, .5, 0, 0, 1), nrow=5, ncol=5)
canonical
##      [,1] [,2] [,3] [,4] [,5]
## [1,] 0.25 0.75 0.00  0.0  0.0
## [2,] 0.00 0.00 0.00  0.5  0.5
## [3,] 0.00 0.75 0.25  0.0  0.0
## [4,] 0.00 0.00 0.00  1.0  0.0
## [5,] 0.00 0.00 0.00  0.0  1.0

From this, we can obtain the matrix Q, shown below

Q = matrix(c(.25, 0, 0, .75, 0, .75, 0, 0, .25), nrow=3, ncol=3)
Q
##      [,1] [,2] [,3]
## [1,] 0.25 0.75 0.00
## [2,] 0.00 0.00 0.00
## [3,] 0.00 0.75 0.25

As well as the matrix 1-Q, which can be used to identify N:

matrix_1_Q = matrix(c(.75, 1, 1, -.25, 1, -.25, 1, 1, .75), nrow=3, ncol=3)
matrix_1_Q
##      [,1]  [,2] [,3]
## [1,] 0.75 -0.25 1.00
## [2,] 1.00  1.00 1.00
## [3,] 1.00 -0.25 0.75

To identify N, we take take the -1 power of the 1-Q matrix

library(matrixcalc) 
N = matrix.power(matrix_1_Q, -1)
round(N, 2)
##       [,1] [,2]  [,3]
## [1,] -1.78 0.11  2.22
## [2,] -0.44 0.78 -0.44
## [3,]  2.22 0.11 -1.78

t = Nc

c = matrix(c(1, 1, 1))
t = N%*%c
round(t, 2)
##       [,1]
## [1,]  0.56
## [2,] -0.11
## [3,]  0.56

From the canonical form, R =

R = matrix(c(.5, 0, 0, 0, 0, .5), nrow=3, ncol=2)
R
##      [,1] [,2]
## [1,]  0.5  0.0
## [2,]  0.0  0.0
## [3,]  0.0  0.5

B = NR, therefore B =

B = N %*% R
round(B, 2)
##       [,1]  [,2]
## [1,] -0.89  1.11
## [2,] -0.22 -0.22
## [3,]  1.11 -0.89