I have difficulties to determine the stationary distribution of a markov model. I start to understand the theory and connections: Given a stochastic matrix, to dermine the stationary distribution we need to find the eigen vector for the largest eigenvalue (that is 1)
I started with generating a stochastic matrix
set.seed(6534)
stoma <- matrix(abs(rnorm(25)), nrow=5, ncol=5)
stoma <- (stoma)/rowSums(stoma) # that should make it a stochastic matrix rowSums(stoma) == 1
Afterward I use the R eigen function
ew <- eigen(stoma)
But I dont understand the results
> ew
$values
[1] 1.000000e+00+0.000000e+00i -6.038961e-02+0.000000e+00i -3.991160e-17+0.000000e+00i
[4] -1.900754e-17+1.345763e-17i -1.900754e-17-1.345763e-17i
$vectors
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4472136+0i 0.81018968+0i 0.3647755+0i -0.0112889+0.1658253i -0.0112889-0.1658253i
[2,] -0.4472136+0i 0.45927081+0i -0.7687393+0i 0.5314923-0.1790588i 0.5314923+0.1790588i
[3,] -0.4472136+0i 0.16233945+0i 0.2128250+0i -0.7093859+0.0000000i -0.7093859+0.0000000i
[4,] -0.4472136+0i -0.09217315+0i 0.4214660+0i -0.1305497-0.1261247i -0.1305497+0.1261247i
[5,] -0.4472136+0i -0.31275073+0i -0.2303272+0i 0.3197321+0.1393583i 0.3197321-0.1393583i
The vector for the largest value (1) has all the same component value "-0.4472136". Even if I change the seed, to draw different number, I get the same values again. What do i miss ? Why the components of the eigenvector are all eqaul? Why they do not sum up to 1 - since this should be a stationary distribtuion?
Thank you for your help!
It's a confusion between the left and right eigenvectors. Try eigen(t(stoma)).
I found it useful to try out the example from the wikipedia article on stochastic matrices:
p <- matrix(c(0,0,1/4,0,0,0,0,1/4,0,0,
1/2,1,0,1/2,0,0,0,1/4,0,0,1/2,0,1/4,1/2,1),
ncol=5)
p
[,1] [,2] [,3] [,4] [,5]
## [1,] 0.00 0.00 0.5 0.00 0.50
## [2,] 0.00 0.00 1.0 0.00 0.00
## [3,] 0.25 0.25 0.0 0.25 0.25
## [4,] 0.00 0.00 0.5 0.00 0.50
## [5,] 0.00 0.00 0.0 0.00 1.00
Check it's a stochastic matrix:
rowSums(p)
## [1] 1 1 1 1 1
Eigenvalues:
zapsmall(eigen(p)$values)
## [1] 1.0000000 0.7071068 -0.7071068 0.0000000 0.0000000
Eigenvectors:
print(zapsmall(eigen(p)$vectors),digits=3)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.447 0.354 0.354 -0.802 -0.609
## [2,] 0.447 0.707 0.707 0.535 -0.167
## [3,] 0.447 0.500 -0.500 0.000 0.000
## [4,] 0.447 0.354 0.354 0.267 0.776
## [5,] 0.447 0.000 0.000 0.000 0.000
(Results as yours but sign arbitrary flipped. R scales eigenvectors so that the sum of squares of each column is 1: sqrt(1/5) is approx. 0.447 ...)
You're looking for the other (left?) eigenvectors:
print(zapsmall(eigen(t(p))$vectors),digits=3)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 -0.149 0.3011 -0.5 0.707
## [2,] 0 -0.149 0.3011 0.5 0.000
## [3,] 0 -0.422 -0.8517 0.0 0.000
## [4,] 0 -0.149 0.3011 -0.5 -0.707
## [5,] 1 0.869 -0.0517 0.5 0.000