#------------------------------------------------------ # # Three state MCMC (cfr Andrieu at al.) # -> solution by eigenvector # # GdA 3/4/2016 #------------------------------------------------------ T <- rbind( c(0, 1, 0), c(0, 0.1, 0.9), c(0.6, 0.4, 0) ) # 1) very compact-------------------------------------------------- # find index of eigenvector 1 ind <- which (round(abs(eigen(t(T))$values), 10) == 1) # (*) # rescale eigenvector such the the sum is unitary p <- abs(eigen(t(T))$vector[,ind] / sum(eigen(t(T))$vector[,ind]) ) # 2) expanded form vals <- eigen(t(T))$values vecs <- abs( eigen(t(T))$vectors ) ind <- which (round(abs(vals), 10) == 1) # (*) (p <- vecs[,ind] / sum (vecs[,ind]) ) # (*) the rounding is needed because the eigenvector found # is not an integer!