How to find when a matrix converges with a cycle

Question

I have been given a matrix:

P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE) 

Using functions, mpow , rows_equal , matrices_equal . I want to find when P^n converges, in other words, what n is when all rows are equal in the matrix and when P^n = P^(n+1) .

Just looking at the functions i , it was possible to deduce that around n=19-21 matrix will converge.

Although, I want to find the correct n using a loop. Here are the mpow , rows_equal and matrices_equal functions. I know that they can be written in different ways, but please keep them as they are.

 mpow <- function(P, n, d=4) { if (n == 0) diag(nrow(P))) else if (n== 1) P else P %*% mpow(P, n - 1)) } rows_equal <- function(P, d = 4) { P_new <- trunc(P * 10^d) for (k in 2:nrow(P_new)) { if (!all(P_new[1, ] == P_new[k, ])) { return(FALSE)} } return(TRUE) } matrices_equal <- function(A, B, d = 4) { A_new <- trunc(A * 10^d) B_new <-trunc(B * 10^d) if (all(A_new == B_new)) TRUE else FALSE } 

Now, to write a loop, we have to do something line by line:

First, create this function:

 when_converged <- function(P) {...} 

and

 for (n in 1:50) 

Try when t.ex n = 50.

Although I do not know how to write code for this correctly, can someone help me with this?

Thanks for reading my question.

Answer 1

Actually, a much better way to do this is:

 ## transition probability matrix P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE) ## a function to find stationary distribution stydis <- function(P, tol = 1e-16) { n <- 1; e <- 1 P0 <- P ## transition matrix P0 while(e > tol) { P <- P %*% P0 ## resulting matrix P e <- max(abs(sweep(P, 2, colMeans(P)))) n <- n + 1 } cat(paste("convergence after",n,"steps\n")) P[1, ] } 

Then when you call the function:

 stydis(P) # convergence after 71 steps # [1] 0.002590674 0.025906736 0.116580311 0.310880829 0.272020725 0.272020725 

The stydis function is essentially continuous:

 P <- P %*% P0 

until convergence P reached. Convergence is numerically determined by the norm L1 of the mismatch matrix:

 sweep(P, 2, colMeans(P)) 

The norm L1 is the maximum, absolute value of all matrix elements. When the norm L1 falls below 1e-16 , convergence occurs.

As you can see, convergence takes 71 steps. Now we can get faster "convergence" by controlling tol (tolerance):

 stydis(P, tol = 1e-4) # convergence after 17 steps # [1] 0.002589361 0.025898057 0.116564506 0.310881819 0.272068444 0.271997814 

But if you check:

 mpow(P, 17) # [,1] [,2] [,3] [,4] [,5] [,6] # [1,] 0.002589361 0.02589806 0.1165645 0.3108818 0.2720684 0.2719978 # [2,] 0.002589415 0.02589722 0.1165599 0.3108747 0.2720749 0.2720039 # [3,] 0.002589738 0.02589714 0.1165539 0.3108615 0.2720788 0.2720189 # [4,] 0.002590797 0.02590083 0.1165520 0.3108412 0.2720638 0.2720515 # [5,] 0.002592925 0.02592074 0.1166035 0.3108739 0.2719451 0.2720638 # [6,] 0.002588814 0.02590459 0.1166029 0.3109419 0.2720166 0.2719451 

Only the first 4 digits match, since you put tol = 1e-4 .

A floating point number has a maximum of 16 digits, so I suggest you use tol = 1e-16 for a reliable convergence test.