############################################################################
#         R FUNCTION TO COMPUTE MARKOV CHAIN STATIONARY DISTRIBUTIONS
#
#         (by Jeffrey S. Rosenthal, probability.ca, 2008)
# 
# This simple R function, statdist, computes the stationary distribution of
# a finite Markov chain, given its transition probabilities, input as a
# single n x n array of transition rows.  (Released under GPL.)
############################################################################

# statdist: compute stationary distribution of Markov transition probabilities
#           (the probabilities should be input as a single n x n array)
statdist = function(transrows) {
nn = sqrt(length(transrows))
if (nn != floor(nn))
  stop("Error, length is not a perfect square.")
transmatt = matrix(transrows, nrow=nn)
transmat <<- t( transmatt )   # (exported)
for (i in 1:nn)
  if (sum(transmat[i,]) != 1)
    stop("Error, row ", i, " does not sum to 1.")
v = eigen(transmatt)$values
for (i in nn:1)
  if (abs(v[i]-1) < 1e-10)
    unitindex = i;
u = Re( eigen(transmatt)$vectors[,unitindex] )
return( u/sum(u) )
}

# EXAMPLE:
mytransrows = c(1/2, 1/2, 1/3, 2/3)
pi = statdist(mytransrows)
P = transmat
print(pi)
print(pi %*% P)