## illustration of how to find stationary distributions via eigendecomposition

# transition matrix for a two-state chain
P_twostate <- function(p, q) {
  matrix(c(1 - p, p,
           q, 1 - q),
         nrow = 2, byrow = TRUE)
}


# non-oscillating two-state chain
P1 <- P_twostate(0.8, 0.4)

# stationary distribution (via eigendecomposition)
eigen(t(P1))
v1 <- eigen(t(P1))$vectors[, 1]
pi1 <- v1/sum(v1)
pi1

# verify
pi1 %*% P1

# note this matches limiting distribution (analytic, from notes)
c(0.8/(0.8 + 0.4), 0.4/(0.8 + 0.4))

# oscillating two-state chain (no limit)
P2 <- P_twostate(1, 1)

# stationary distribution is easy to find: (0.5, 0.5)
c(0.5, 0.5) %*% P2

# via eigendecomposition (matches previous result)
eigen(t(P2))
v2 <- eigen(t(P2))$vectors[,1]
pi2 <- v2/sum(v2)
pi2