Markov Chains

import numpy as np

P = np.array([[0.9, 0.1],
              [0.2, 0.8]])

# Method 1: Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eig(P.T)
idx = np.argmin(np.abs(eigenvalues - 1))
stationary = np.real(eigenvectors[:, idx])
stationary = stationary / stationary.sum()
print(f"Stationary distribution: {stationary}")
# Output: [0.6667 0.3333]

# Method 2: Power iteration
pi = np.array([0.5, 0.5])
for _ in range(100):
    pi = pi @ P
print(f"Stationary distribution: {pi}")
# Output: [0.6667 0.3333]

import numpy as np
from typing import List, Tuple, Optional

class MarkovChain:
    """A complete Markov chain implementation."""

    def __init__(self, transition_matrix: np.ndarray, state_names: Optional[List[str]] = None):
        self.P = np.array(transition_matrix, dtype=float)
        n = self.P.shape[0]
        assert self.P.shape == (n, n), "Matrix must be square"
        assert np.all(self.P >= 0), "All entries must be non-negative"
        assert np.allclose(self.P.sum(axis=1), 1), "Rows must sum to 1"
        self.n = n
        self.state_names = state_names or [str(i) for i in range(n)]

    def step(self, state: int) -> int:
        """Take one step from the given state."""
        return np.random.choice(self.n, p=self.P[state])

    def simulate(self, start: int, steps: int) -> List[int]:
        """Simulate a path of given length starting from state."""
        path = [start]
        for _ in range(steps):
            path.append(self.step(path[-1]))
        return path

    def n_step_matrix(self, n: int) -> np.ndarray:
        """Compute the n-step transition matrix P^n."""
        return np.linalg.matrix_power(self.P, n)

    def stationary_distribution(self) -> np.ndarray:
        """Find the stationary distribution via eigendecomposition."""
        eigenvalues, eigenvectors = np.linalg.eig(self.P.T)
        idx = np.argmin(np.abs(eigenvalues - 1))
        stationary = np.real(eigenvectors[:, idx])
        return stationary / stationary.sum()

    def is_irreducible(self) -> bool:
        """Check if the chain is irreducible."""
        reach = np.linalg.matrix_power(self.P, self.n - 1)
        return np.all(reach > 0)

    def is_aperiodic(self) -> bool:
        """Check if the chain is aperiodic."""
        from math import gcd
        for i in range(self.n):
            periods = []
            P_n = np.eye(self.n)
            for k in range(1, self.n * self.n + 1):
                P_n = P_n @ self.P
                if P_n[i, i] > 0:
                    periods.append(k)
            if periods:
                g = periods[0]
                for p in periods[1:]:
                    g = gcd(g, p)
                if g > 1:
                    return False
        return True

    def is_ergodic(self) -> bool:
        """Check if the chain is ergodic (irreducible + aperiodic)."""
        return self.is_irreducible() and self.is_aperiodic()

    def absorbing_states(self) -> List[int]:
        """Find all absorbing states."""
        return [i for i in range(self.n) if self.P[i, i] == 1.0]

    def fundamental_matrix(self) -> np.ndarray:
        """Compute the fundamental matrix for absorbing chains."""
        absorbing = self.absorbing_states()
        transient = [i for i in range(self.n) if i not in absorbing]
        if not transient:
            return np.eye(0)
        Q = self.P[np.ix_(transient, transient)]
        return np.linalg.inv(np.eye(len(transient)) - Q)

    def steady_state_distribution(self, start: int, steps: int = 10000) -> np.ndarray:
        """Estimate stationary distribution via simulation."""
        counts = np.zeros(self.n)
        state = start
        for _ in range(steps):
            state = self.step(state)
            counts[state] += 1
        return counts / counts.sum()


# === Example Usage ===

# Weather model
P_weather = [[0.6, 0.3, 0.1],
             [0.3, 0.4, 0.3],
             [0.2, 0.3, 0.5]]

mc = MarkovChain(P_weather, ["Sunny", "Cloudy", "Rainy"])

print("Stationary distribution:", mc.stationary_distribution())
print("Is ergodic:", mc.is_ergodic())
print("10-step transitions:\n", mc.n_step_matrix(10))

# Simulate 30 days
path = mc.simulate(0, 30)  # Start from Sunny
weather_counts = [path.count(i) for i in range(3)]
print("Weather distribution:", [c/len(path) for c in weather_counts])

Stationary distribution: Solve $\pi P = \pi$ with $\pi_1 + \pi_2 = 1$:

• $\pi_1 = 0.7\pi_1 + 0.4\pi_2 \implies 0.3\pi_1 = 0.4\pi_2 \implies \pi_1 = \frac{4}{3}\pi_2$
• $\frac{4}{3}\pi_2 + \pi_2 = 1 \implies \frac{7}{3}\pi_2 = 1 \implies \pi_2 = \frac{3}{7}$
• $\boldsymbol{\pi} = (4/7, 3/7) \approx (0.571, 0.429)$

3-step transition: $P^3 = P \cdot P \cdot P = \begin{pmatrix} 0.583 & 0.417 \\ 0.556 & 0.444 \end{pmatrix}$

Transition matrix (states 1, 2, 3, 4):

$Q = \begin{pmatrix} 0 & 0.5 & 0.5 \\ 0.5 & 0 & 0 \\ 0.5 & 0 & 0 \end{pmatrix}$, $R = \begin{pmatrix} 0 \\ 0.5 \\ 0.5 \end{pmatrix}$

Probability of absorption into cell 4 from cell 1: $0.5$

Transition matrix:

Solving $\pi P = \pi$: By symmetry, $\pi_0 = \pi_3$ and $\pi_1 = \pi_2$. From balance equations: $\pi_0 = 0.5\pi_1$ and $\pi_1 = \pi_0 + 0.5\pi_2 = \pi_0 + 0.5\pi_1$, so $\pi_1 = 2\pi_0$.

Normalization: $\pi_0 + 2\pi_0 + 2\pi_0 + \pi_0 = 1 \implies 6\pi_0 = 1$

$\boldsymbol{\pi} = (1/6, 1/3, 1/3, 1/6)$

Adjacency matrix: $L = \begin{pmatrix} 0 & 1 & 1 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}$

Column-stochastic: $M = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}$

PageRank equation: $\vec{r} = 0.05 \cdot \frac{1}{3}\vec{1} + 0.85 \cdot M\vec{r}$

Solving: $\vec{r} = (1/3, 1/3, 1/3)$ (uniform due to the cycle A?B?C?A)