> P <- matrix(data = 0, nrow = 2, ncol = 2)
> a <- 1/3
> P[1, ] <- c(a, 1 - a)
> P[2, ] <- c(1, 0)
> temp <- P
> for (i in 1:100) {
+     temp <- temp %*% P
+ }
> round(temp[1, 1] - 1/(2 - a), 2)
[1] 0
> P <- matrix(data = 0, nrow = 2, ncol = 2)
> a <- 1/4
> P[1, ] <- c(a, 1 - a)
> P[2, ] <- c(1, 0)
> temp <- P
> for (i in 1:100) {
+     temp <- temp %*% P
+ }
> round(temp[1, 1] - 1/(2 - a), 2)
[1] 0
> P <- matrix(data = 0, nrow = 2, ncol = 2)
> a <- 1/2
> P[1, ] <- c(a, 1 - a)
> P[2, ] <- c(1, 0)
> temp <- P
> for (i in 1:100) {
+     temp <- temp %*% P
+ }
> round(temp[1, 1] - 1/(2 - a), 2)
[1] 0
> P <- matrix(data = 0, nrow = 3, ncol = 3)
> P[1, ] <- c(0, 2/3, 1/3)
> P[2, ] <- c(1/4, 0, 3/4)
> P[3, ] <- c(4/5, 1/5, 0)
> coefs <- t(P) - diag(3)
> wts <- solve(rbind(c(1, 1, 1), coefs)[1:3, ], c(1, 0, 0))
> W <- matrix(data = rep(wts, 3), nrow = 3, ncol = 3, byrow = T)
> Z <- solve(diag(3) - P + W)
> M <- matrix(data = 0, nrow = 3, ncol = 3)
> for (i in 1:3) {
+     for (j in 1:3) {
+         M[i, j] <- (Z[j, j] - Z[i, j])/wts[j]
+     }
+ }
> M
         [,1]     [,2] [,3]
[1,] 0.000000 1.818182  2.0
[2,] 2.058824 0.000000  1.5
[3,] 1.411765 2.454545  0.0
> C <- matrix(data = rep(1, 9), nrow = 3, ncol = 3)
> I <- diag(3)
> R <- round(C - (I - P) %*% M)
> print(R)
     [,1] [,2] [,3]
[1,]    3    0    0
[2,]    0    3    0
[3,]    0    0    3