MarkovChain_Tutorial.md

Tutorial: Introduction to Markov Chains

1. What is a Markov Chain?

A Markov chain is a mathematical model for a system that evolves step by step, where the future depends only on the present state, not on the full history. This property is called the Markov property:

$$ Pr(X_{t+1} \mid X_t, X_{t-1}, \ldots, X_0) = \Pr(X_{t+1} \mid X_t) $$

• States: possible conditions the system can be in.
• Transitions: probabilities of moving from one state to another.

2. Transition Matrix

We represent a Markov chain with a transition matrix P: $$ P = \begin{bmatrix} P_{11} & P_{12} & \dots & P_{1n} \ P_{21} & P_{22} & \dots & P_{2n} \ \vdots & \vdots & \ddots & \vdots \ P_{n1} & P_{n2} & \dots & P_{nn} \end{bmatrix}, $$ where $P_{ij} = \Pr(X_{t+1}=j \mid X_t=i)$ • Each row sums to 1 (probabilities). • If we start with a distribution vector x_t, the next step is:

$ x_{t+1} = x_t P$

3. Example: Weather Model

Suppose weather can be Sunny (S) or Rainy (R). • If it’s sunny today, tomorrow is sunny with probability 0.9, rainy with 0.1. • If it’s rainy today, tomorrow is sunny with 0.5, rainy with 0.5.

Transition matrix: $$ P = \begin{bmatrix} 0.9 & 0.1 \ 0.5 & 0.5 \end{bmatrix}. $$ If we start with $x_0 = [1,0]$ (100% sunny), then: $$ x_1 = x_0 P = [1,0] \begin{bmatrix} 0.9 & 0.1 \ 0.5 & 0.5 \end{bmatrix} = [0.9, 0.1]. $$ So tomorrow: 90% chance sunny, 10% rainy.

4. Steady State (Stationary Distribution)

In many Markov chains, as time goes on, the distribution converges to a fixed vector $\pi$ such that: $$ \pi = \pi P, \quad \sum_i \pi_i = 1. $$ This is the long-run probability of being in each state.

For the weather example, solving gives: $\pi = [0.83,, 0.17].$ (So in the long run: ~83% sunny, ~17% rainy, regardless of the initial condition.)

5. What Problems Do Markov Chains Solve?

They are used when you want to model random processes that evolve over time: • Queues and systems: waiting lines, network traffic, customer arrivals. • Random walks: modeling diffusion, stock prices (simplified). • Google PageRank: web surfing is a Markov chain over websites. • Biology: DNA sequence modeling. • Speech recognition / NLP: word sequences modeled as Markov processes. • Operations research: logistics, inventory, and reliability models.

Core idea: Markov chains help us predict the future distribution of states, and understand long-term behavior.

6. How to Compute in Practice
   1. Build transition matrix P from data or assumptions.
   2. Forecast future states: compute $x_t P^k$.
   3. Find steady state: solve $\pi = \pi P$. • Either as eigenvector of $P^\top$ with eigenvalue 1. • Or iterate numerically until convergence.

# markov_demo.py
import numpy as np
import matplotlib.pyplot as plt

def simulate_chain(P, x0, steps):
    x = np.asarray(x0, dtype=float)
    s = x.sum()
    x = x if s == 0 else x / s
    hist = [x.copy()]
    for _ in range(steps):
        x = x @ P
        total = x.sum()
        if total != 0:
            x = x / total
        hist.append(x.copy())
    return np.array(hist)

def power_iteration_stationary(P, steps=1000, x0=None):
    n = P.shape[0]
    if x0 is None:
        x = np.ones(n) / n
    else:
        x = np.asarray(x0, dtype=float)
        s = x.sum()
        x = x if s == 0 else x / s
    for _ in range(steps):
        x = x @ P
        total = x.sum()
        if total != 0:
            x = x / total
    return x

def left_eigenvector_stationary(P):
    vals, vecs = np.linalg.eig(P.T)
    idx = np.argmin(np.abs(vals - 1))
    v = np.real(vecs[:, idx])
    v = np.abs(v)
    return v / v.sum()

# -------- Weather example --------
P_weather = np.array([[0.9, 0.1],
                      [0.5, 0.5]], dtype=float)
x0_weather = np.array([1.0, 0.0])

hist_weather = simulate_chain(P_weather, x0_weather, steps=30)

plt.figure(figsize=(6,4))
plt.plot(hist_weather[:, 0], label="Sunny")
plt.plot(hist_weather[:, 1], label="Rainy")
plt.title("Weather: Distribution Over Time")
plt.xlabel("Step")
plt.ylabel("Probability")
plt.legend()
plt.show()

pi_power_w = power_iteration_stationary(P_weather, steps=1000)
pi_eig_w = left_eigenvector_stationary(P_weather)
print("Weather stationary (power):", pi_power_w)
print("Weather stationary (eig):  ", pi_eig_w)

plt.figure(figsize=(4,3))
plt.bar(["Sunny","Rainy"], pi_power_w)
plt.title("Weather: Stationary Distribution")
plt.ylabel("Probability")
plt.show()

# -------- Parcel example (A–E) --------
nodes = ["A","B","C","D","E"]
P_parcel = np.array([
    [0.0, 0.3, 0.7, 0.0, 0.0],
    [0.0, 0.0, 0.0, 1.0, 0.0],
    [0.0, 0.0, 0.0, 0.6, 0.4],
    [0.0, 0.0, 0.0, 0.0, 1.0],
    [1.0, 0.0, 0.0, 0.0, 0.0],
], dtype=float)

x0_counts = np.array([100, 10, 20, 0, 0], dtype=float)
x0_prob = x0_counts / x0_counts.sum()

hist_parcel = simulate_chain(P_parcel, x0_prob, steps=25)

plt.figure(figsize=(7,4))
for i, name in enumerate(nodes):
    plt.plot(hist_parcel[:, i], label=name)
plt.title("Parcel: Distribution Over Time (Normalized Probabilities)")
plt.xlabel("Step")
plt.ylabel("Probability")
plt.legend()
plt.show()

x1_counts = x0_counts @ P_parcel
print("Parcel x0 (counts):", x0_counts)
print("Parcel x1 (counts):", x1_counts)

plt.figure(figsize=(5,3))
plt.bar(nodes, x0_counts)
plt.title("Parcel: x0 (Counts)")
plt.ylabel("Count")
plt.show()

plt.figure(figsize=(5,3))
plt.bar(nodes, x1_counts)
plt.title("Parcel: x1 after one step (Counts)")
plt.ylabel("Count")
plt.show()

pi_parcel = power_iteration_stationary(P_parcel, steps=2000)
print("Parcel stationary distribution (probabilities):")
for name, p in zip(nodes, pi_parcel):
    print(f"{name}: {p:.4f}")

plt.figure(figsize=(5,3))
plt.bar(nodes, pi_parcel)
plt.title("Parcel: Stationary Distribution")
plt.ylabel("Probability")
plt.show()

Weather stationary (power): [0.83333333 0.16666667]
Weather stationary (eig):   [0.83333333 0.16666667]

Parcel x0 (counts): [100.  10.  20.   0.   0.]
Parcel x1 (counts): [ 0. 30. 70. 22.  8.]

Parcel stationary distribution (probabilities):
A: 0.2688
B: 0.0806
C: 0.1882
D: 0.1935
E: 0.2688