Indagamos por todo el mundo online para así traerte la solución para tu problema, si tienes alguna pregunta puedes dejarnos tu inquietud y contestamos sin falta.
Solución:
Bueno, para usar la descomposición Eigen, necesitamos trabajar con t(P).
La definición de una matriz de probabilidad de transición difiere entre probabilidad/estadística y álgebra lineal. En estadísticas todas las filas de P suma a 1, mientras que en álgebra lineal, todas las columnas de P suma a 1. Así que en lugar de eigen(P)nosotros necesitamos eigen(t(P)):
e <- Re(eigen(t(P))$vectors[, 1])
e / sum(e)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708
Como podemos ver, solo hemos usado el primer vector Eigen, es decir, el vector Eigen del valor Eigen más grande. Por lo tanto, no hay necesidad de calcular todos los valores propios/vectores usando eigen. El método de potencia se puede utilizar para encontrar un vector propio del valor propio más grande. Implementemos esto en R:
stydis1 <- function (A)
n <- dim(A)[1L]
## checking
if (any(.rowSums(A, n, n) != 1))
stop (" 'A' is not a Markov matrix")
## implement power method
e <- runif (n)
oldnorm <- sqrt(c(crossprod(e)))
repeat
e <- crossprod(A, e)
newnorm <- sqrt(c(crossprod(e)))
if (abs(newnorm / oldnorm - 1) < 1e-8) break
e <- e / newnorm
oldnorm <- newnorm
## rescale `e` so that it sums up to 1
c(e / sum(e))
stydis1 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708
Y el resultado es correcto.
De hecho, no tenemos que explotar la descomposición de Eigen. Podemos ajustar el método utilizado en su segunda pregunta vinculada. Allí, tomamos energía de matriz que es costosa como comentaste; pero ¿por qué no volver a convertirlo en una multiplicación matriz-vector?
stydis2 <- function (A)
n <- dim(A)[1L]
## checking
if (any(.rowSums(A, n, n) != 1))
stop (" 'A' is not a Markov matrix")
## direct computation
b <- A[1, ]
oldnorm <- sqrt(c(crossprod(b)))
repeat
b <- crossprod(A, b)
newnorm <- sqrt(c(crossprod(b)))
if (abs(newnorm / oldnorm - 1) < 1e-8) break
oldnorm <- newnorm
## return stationary distribution
c(b)
stydis2 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708
Partimos de una distribución inicial arbitraria, digamos A[1, ], y aplicar iterativamente la matriz de transición hasta que la distribución converja. De nuevo, el resultado es correcto.
tu vector y = Re(eigen(P)$vectors[, 1]) no es una distribución (ya que no suma uno) y resuelve P'y = yno x'P = x. La solución de sus preguntas y respuestas vinculadas resuelve aproximadamente lo último:
x = c(0.00259067357512953, 0.0259067357512953, 0.116580310880829,
0.310880829015544, 0.272020725388601, 0.272020725388601)
all(abs(x %*% P - x) < 1e-10) # TRUE
Al transponer P, puede usar su enfoque de valor propio:
x2 = Re(eigen(t(P))$vectors[, 1])
x2 <- x2/sum(x2)
(function(x) all(abs(x %*% P - x) < 1e-10))(
x2
) # TRUE
Sin embargo, está encontrando un vector estacionario diferente en este caso.
Nos encantaría que puedieras comunicar esta división si te valió la pena.