xx <- matrix(c(16,3,2,13,5,10,11,8,9,6,7,12,4,15,14,1),nrow = 4,ncol = 4,byrow = T)
xxx <- xx/rowSums(xx)
xeig <- eigen(xxx)
xvecs <- xeig$vectors
lvecs <- MASS::ginv(xvecs)
pi_eq <- lvecs[1,]/sum(lvecs[1,])
pi_eq