rmarkov-chainsmarkov

Source code for calculation of stationary distribution in R


Take a look at this link.

I am trying to understand the following source code meant for finding stationary distribution of a matrix:

# Stationary distribution of discrete-time Markov chain
# (uses eigenvectors)
stationary <- function(mat)
{
    x = eigen(t(mat))$vectors[,1]
    as.double(x/sum(x))
}

I tested the following source code myself:

> rm(list=ls())
>
> P <- matrix(c(0.66, 0.34,
+               0.66, 0.34), nrow=2, ncol=2, byrow = TRUE)
>
> x <- eigen(t(P))
> x$values
[1] 1 0

$vectors
          [,1]       [,2]
[1,] 0.8889746 -0.7071068
[2,] 0.4579566  0.7071068

> y <- x$vectors[,1]
> y
[1] 0.8889746 0.4579566
>  

looks like the command

y <- x$vectors[,1]

is selecting the 1st column of the matrix.

Why wasn't that simply written like the following?

# Stationary distribution of discrete-time Markov chain
# (uses eigenvectors)
stationary <- function(mat)
{
    x = eigen(t(mat))
    y = x[,1]
    as.double(y/sum(y))
}

What was the reason for introduction of a dollar sign and vector keyword?


Solution

  • Let's test out your proposal:

    > P <- matrix(c(0.66, 0.34, 0.66, 0.34), nrow=2, ncol=2, byrow = TRUE)
    > x <- eigen(t(P))
    > print(x)
    eigen() decomposition
    $values
    [1] 1 0
    
    $vectors
              [,1]       [,2]
    [1,] 0.8889746 -0.7071068
    [2,] 0.4579566  0.7071068
    
    > y = x[,1]
    

    This would produce the following error message:

    Error in x[, 1] : incorrect number of dimensions
    

    eigen returns a named list, with eigenvalues named values and eigenvectors named vectors. To access this component of the list. we use the dollar sign. Hence, that is why the code x$vectors which extract the matrix.