Wie zu finden, wenn eine Matrix mit einer Schleife konvergiert

Question

Ich habe eine Matrix gegeben worden:

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) 

die Funktionen verwenden, mpow, rows_equal, matrices_equal. Ich möchte herausfinden, wenn P^n konvergiert, mit anderen Worten, was n ist, wenn alle Zeilen in der Matrix gleich sind und wenn P^n = P^(n+1).

Durch nur Blick auf die Funktionen i haben es geschafft zu schließen, dass um n=19-21 die Matrix konvergieren wird.

Obwohl ich die richtige n mit einer Schleife finden wollen. Hier sind die Funktionen mpow, rows_equal und matrices_equal. Ich weiß, dass sie anders geschrieben werden können, aber behalte sie so, wie sie sind.

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 
    } 

Nun, die Schleife zu schreiben, sollten wir es etwas entlang der Linien von tun:

Zuerst eine Funktion wie so erstellen:

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

und

for (n in 1:50) 

Um zu versuchen, wenn t.ex n = 50.

Obwohl ich nicht weiß, wie man den Code richtig schreibt, kann mir dann jemand helfen?

Vielen Dank für das Lesen meiner Frage.

Answer 1

Eigentlich ein viel besserer Weg, dies zu tun:

## 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, ] 
    } 

Dann, wenn Sie die Funktion aufrufen:

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

Die Funktion stydis, im Wesentlichen kontinuierlich tut:

P <- P %*% P0 

bis Konvergenz von P erreicht ist. Konvergenz wird numerisch durch die Norm L1 der Diskrepanz Matrix bestimmt:

sweep(P, 2, colMeans(P)) 

Die L1-Norm ist der maximale, absolute Wert aller Matrixelemente. Wenn die L1-Norm unter 1e-16 fällt, tritt eine Konvergenz auf.

Wie Sie sehen können, dauert die Konvergenz 71 Schritte. Jetzt können wir schneller „Konvergenz“ erhalten durch Steuern tol (Toleranz):

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

Aber wenn Sie überprüfen:

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 

Nur die ersten vier Ziffern sind die gleichen, wie Sie tol = 1e-4 setzen.

Eine Gleitkommazahl hat maximal 16 Stellen. Daher würde ich Ihnen empfehlen, tol = 1e-16 für einen zuverlässigen Konvergenztest zu verwenden.