1a Discrete-time Markov chain

A Markov chain on a finite state space $\mathcal{S}=\{1,\dots,n\}$ is governed by a row-stochastic transition matrix $P$, where $P_{ij}=\Pr(X_{t+1}=j \mid X_t=i)$. The Markov property says the future depends on the present only — not the path taken to get there:

If the row vector $\mu_t$ holds the distribution over states at time $t$, then one step forward is just a vector–matrix product, $\mu_{t+1} = \mu_t P$, and $k$ steps is $\mu_{t+k} = \mu_t P^{k}$. For an irreducible, aperiodic chain this converges to a unique stationary distribution $\pi$ — the row vector fixed by $P$:

1b Exponential inter-arrivals & the Poisson process

The queue runs in continuous time. Customers arrive as a Poisson process of rate $\lambda$, which is equivalent to saying the gaps between arrivals are i.i.d. $\mathrm{Exponential}(\lambda)$. Service times are also $\mathrm{Exponential}(\mu)$. The exponential is the only continuous distribution that is memoryless:

Memorylessness is what makes the queue a (continuous-time) Markov chain on the number of customers in the system: how long the current arrival or service has already taken tells you nothing about how much longer it will take.

1c The M/M/1 queue & Little's law

M/M/1 = Markovian arrivals, Markovian service, 1 server, infinite FIFO buffer. Define the utilization (traffic intensity) $\rho=\lambda/\mu$. The queue is stable only when $\rho<1$; otherwise work arrives faster than it can be cleared and the line grows without bound. For a stable M/M/1 queue the steady-state metrics have clean closed forms:

$L$ is the mean number in the system, $L_q$ the mean number waiting, $W$ the mean time in system, and $W_q$ the mean waiting time before service. These four are not independent — they are stitched together by Little's law, which holds for essentially any stable queueing system, M/M/1 or not:

3.1 Build the chain & solve $\pi$ by eigenvector

Define the transition matrix and assert each row sums to 1. The stationary distribution is the eigenvector of $P^\top$ whose eigenvalue is closest to 1, renormalized to sum to one.

# mm1_markov.py — part 1 of 3

import numpy as np

# 4-page "web navigation" chain: P[i,j] = P(go to page j | on page i)
P = np.array([
    [0.10, 0.40, 0.30, 0.20],
    [0.20, 0.10, 0.40, 0.30],
    [0.25, 0.25, 0.10, 0.40],
    [0.30, 0.30, 0.30, 0.10],
])
assert np.allclose(P.sum(axis=1), 1.0)   # rows are probability distributions

def stationary_eig(P):
    """pi P = pi  <=>  pi is the left eigenvector of P for eigenvalue 1."""
    vals, vecs = np.linalg.eig(P.T)            # left eigvecs of P = eigvecs of P^T
    i = np.argmin(np.abs(vals - 1.0))         # pick eigenvalue closest to 1
    pi = np.real(vecs[:, i])
    return pi / pi.sum()                        # normalise to a distribution

pi_eig = stationary_eig(P)

3.2 Solve $\pi$ again by power iteration

Start from any distribution and repeatedly apply $\mu \leftarrow \mu P$. For an irreducible aperiodic chain this converges geometrically to $\pi$. This is the algorithm behind PageRank.

# mm1_markov.py — part 1 (cont.)

def stationary_power(P, tol=1e-14, max_iter=100_000):
    """Repeatedly push a distribution through P until it stops moving."""
    n = P.shape[0]
    mu = np.full(n, 1.0 / n)                  # start uniform
    for _ in range(max_iter):
        nxt = mu @ P
        if np.max(np.abs(nxt - mu)) < tol:
            return nxt
        mu = nxt
    return mu

pi_pow = stationary_power(P)

print("pi (eigenvector) :", np.round(pi_eig, 6))
print("pi (power iter)  :", np.round(pi_pow, 6))
print("agree?           :", np.allclose(pi_eig, pi_pow))
print("pi P == pi ?     :", np.allclose(pi_pow @ P, pi_pow))

3.3 Analytic M/M/1 metrics

The closed forms from section 1c, written once:

# mm1_markov.py — part 2 of 3

def mm1_theory(lam, mu):
    """Closed-form steady-state metrics for a stable M/M/1 queue (rho < 1)."""
    assert lam < mu, "unstable: arrivals outpace service (rho >= 1)"
    rho = lam / mu
    return {
        "rho": rho,                       # utilization  = P(server busy)
        "L":   rho / (1 - rho),            # mean number in system
        "Lq":  rho**2 / (1 - rho),         # mean number waiting
        "W":   1 / (mu - lam),             # mean time in system
        "Wq":  rho / (mu - lam),           # mean wait before service
    }

lam, mu = 1.2, 1.5                          # rho = 0.8: busy but stable
theory = mm1_theory(lam, mu)

3.4 Simulate the queue (event-driven, Lindley recursion)

Draw all inter-arrival and service times, accumulate arrival epochs, then walk the Lindley recursion to get each customer's service-start and departure. From those we read $W$ and $W_q$ directly; for $L$ we integrate the number-in-system over a merged, sorted timeline of $+1$ (arrival) and $-1$ (departure) events.

# mm1_markov.py — part 3 of 3

def simulate_mm1(lam, mu, n_cust=200_000, seed=42):
    """Event-driven FIFO single-server queue. Returns measured metrics."""
    rng = np.random.default_rng(seed)
    inter = rng.exponential(1.0 / lam, n_cust)   # gaps ~ Exp(lam)  (Poisson arrivals)
    serve = rng.exponential(1.0 / mu,  n_cust)   # service ~ Exp(mu)
    arrive = np.cumsum(inter)                       # arrival epochs a_k

    start  = np.empty(n_cust)
    depart = np.empty(n_cust)
    for k in range(n_cust):                        # Lindley: d_k = max(a_k, d_{k-1}) + s_k
        s = arrive[k] if k == 0 else max(arrive[k], depart[k - 1])
        start[k]  = s
        depart[k] = s + serve[k]

    W  = (depart - arrive).mean()                   # time in system
    Wq = (start  - arrive).mean()                   # time waiting

    # time-average L: integrate #-in-system over a merged event timeline
    t = np.concatenate([arrive, depart])
    d = np.concatenate([np.ones(n_cust), -np.ones(n_cust)])
    order = np.argsort(t, kind="mergesort")     # stable sort of all events
    t, d = t[order], d[order]
    n_in = np.cumsum(d)                             # number in system after each event
    dt   = np.diff(t)
    L    = np.sum(n_in[:-1] * dt) / (t[-1] - t[0])
    util = serve.sum() / depart[-1]              # fraction of time server is busy
    return {"rho": util, "L": L, "Lq": L - util, "W": W, "Wq": Wq}

sim = simulate_mm1(lam, mu)

# --- compare sim vs theory, and verify Little's law from measured numbers ---
print(f"{'metric':<6}{'theory':>12}{'sim':>12}")
for k in ["rho", "L", "Lq", "W", "Wq"]:
    print(f"{k:<6}{theory[k]:>12.4f}{sim[k]:>12.4f}")
print(f"Little's law  L = lam*W : {lam * sim['W']:.4f} vs L_sim {sim['L']:.4f}")

4a Stationary distribution of the chain

pi (eigenvector) : [0.217292 0.256799 0.2714   0.254509]
pi (power iter)  : [0.217292 0.256799 0.2714   0.254509]
agree?           : True
pi P == pi ?     : True

Both methods give the same $\pi$. Page 3 is the most-visited (≈27.1%), page 1 the least (≈21.7%) — the long-run "PageRank" of this little site. Power iteration reached this from a uniform start in a handful of steps:

4b M/M/1 queue — simulation vs. theory

metric      theory         sim
rho         0.8000      0.7990
L           4.0000      3.9757
Lq          3.2000      3.1767
W           3.3333      3.3110
Wq          2.6667      2.6450
Little's law  L = lam*W : 3.9731 vs L_sim 3.9757