CONTINUOUS-TIME MARKOV CHAINS

Chapter 1: From Probabilities to Rates

We just decided to throw away the tick. But if there is no tick, what number describes “how likely is it to fail”? You cannot say “the chance it fails this instant,” because an instant has zero width, and the probability of a jump in zero time is exactly zero. The probability question has no good answer at an instant.

The number that does survive at an instant is a rate. A rate q is the instantaneous tendency to leave a state, measured per unit time, with units of 1/time. “The machine fails at a rate of 0.2 per hour” is a sentence that makes sense without ever mentioning a tick.

The analogy to lean on: rate is to probability what speed is to distance. You never ask “how far did the car travel in this instant?” — the answer is zero. You ask “how fast is it going?” Speed is the slope of distance; rate is the slope of accumulated jump-probability. Distance accumulates over a window; probability accumulates over a window; speed and rate are the slopes at an instant.

Here is the bridge between the two languages. Over a tiny window of width dt, the probability that the jump has happened is approximately q·dt. This is the linear approximation, and it is the single most-used formula in the whole subject. But it is only the leading term — the honest, all-dt answer is a curve.

That curve is P(jumped by time t) = 1 − e−qt. It starts at 0 (nothing has happened yet), rises, and approaches 1 as t→∞ (eventually it surely happens). Its slope at the origin is exactly q — which is why q·dt works for small dt: you are riding the tangent line.

Where does that curve come from? From the simplest differential equation in the world. Let x(t) be the probability the machine is still up at time t. It leaks toward DOWN at rate q, so it loses a q-fraction of itself per unit time: dx/dt = −q·x. The solution is x(t) = e−qt, and “already jumped” is the complement, 1 − e−qt. Hold onto that scalar ODE — in Chapter 6 it becomes dP/dt = PQ and in Chapter 7 its solution becomes the matrix exponential exp(Qt). The whole lesson is this one-line idea, grown up.

Deriving the curve: separate the variables

Let us actually solve dx/dt = −q·x from scratch so the exponential is earned, not asserted. The equation says “the rate of change of x is proportional to x itself,” and the trick for any such equation is separation of variables: gather all the x’s on one side and all the t’s on the other.

Now integrate both sides. The left integrates to ln x, the right to −q·t + C for an integration constant C: ln x = −q·t + C. Exponentiate both sides to undo the log: x(t) = eC·e−qt.

Fix the constant with the initial condition. At t = 0 the machine is certainly still up, so x(0) = 1. Plugging in: 1 = eC·e0 = eC, hence eC = 1. The survival curve is therefore x(t) = e−qt exactly, with no free parameters left. The “already jumped” probability is the complement F(t) = 1 − e−qt. That is the entire derivation — one separation, one integral, one initial condition.

It is worth naming why this matters beyond Chapter 1. The same three moves — separate, integrate, fix the constant — will reappear in Chapter 6 with matrices, where dP/dt = PQ “separates” into P(t) = exp(Qt). A scalar ODE you can solve by hand becomes a matrix ODE you solve by eigenvalues. Same skeleton, more bones.

Why the hazard is exactly q at every instant

The hazard is the instantaneous jump-rate given you have not yet jumped. It is defined as the density of jump times divided by the survival probability: h(t) = f(t) / S(t), where f(t) is the derivative of F(t) = 1 − e−qt and S(t) = e−qt.

Differentiate F: f(t) = d/dt (1 − e−qt) = q·e−qt. Divide by survival: h(t) = q·e−qt / e−qt = q. The exponentials cancel completely — the hazard is the constant q for all t. This is the “constant hazard” phrase made literal: no matter how long you have already waited, your per-instant chance of jumping is the same q. That cancellation is unique to the exponential, and it is exactly what makes the holding time in Chapter 2 memoryless.

Worked example: q = 0.2 per hour, two windows

Take the canonical fail rate q = 0.2 per hour. We test the linear rule against the truth at two windows.

Small window, dt = 0.1 hr. Linear estimate: q·dt = 0.2 × 0.1 = 0.0200. True probability: 1 − e−0.2×0.1 = 1 − e−0.02. Now e−0.02 = 0.980199, so the true value is 1 − 0.980199 = 0.0198. The estimate 0.0200 versus the truth 0.0198 — an absolute error of 0.0002, about 1%. Excellent.

Large window, dt = 2.0 hr. Linear estimate: q·dt = 0.2 × 2.0 = 0.40. True probability: 1 − e−0.2×2.0 = 1 − e−0.4. Here e−0.4 = 0.6703, so the truth is 1 − 0.6703 = 0.3297. The linear rule said 0.40 but the truth is 0.3297 — it overshoots by ~0.07, a 21% error. Worse, a big enough q·dt can exceed 1, which is nonsense for a probability.

Conclusion, stated numerically: q·dt is the jump probability only for small dt. For finite windows you must use the curve 1 − e−qt. The error of the linear rule shrinks like dt² as dt→0 — halve dt and the error quarters.

Why the error is exactly second-order

Where does that “quarters when you halve” rule come from? Taylor-expand the true curve at the origin. e−x = 1 − x + x²/2 − …, so with x = q·dt:

The first term q·dt is exactly the linear rule. So the error — true minus linear — is everything after it: −(q·dt)²/2 + …, whose leading term is (q·dt)²/2 in magnitude. The error scales as dt², so halving dt multiplies the error by (1/2)² = 1/4 — the quartering, derived.

Check it against the worked numbers. At q = 0.2, dt = 0.1: predicted error ≈ (0.2×0.1)²/2 = (0.02)²/2 = 0.0004/2 = 0.0002, matching the observed 0.0002 to the digit. The sign is negative (the linear rule overshoots), which is why q·dt > 1 − e−q·dt always.

Units check, worked example #1: a slow rate, a long window

Before the quiz, a second units example to cement the dimensional reasoning. A capacitor discharges (a one-way “jump” from charged to empty) at rate q = 0.5 per second; we ask about a window of dt = 2 s. Dimensionally q·dt = (1/s)·(s) = dimensionless, as any probability must be. Numerically the linear rule gives 0.5×2 = 1.0 — which is already the maximum a probability can be, a screaming red flag that the window is far too wide for the tangent. The truth is 1 − e−0.5×2 = 1 − e−1 = 1 − 0.3679 = 0.6321. The linear rule claimed certainty; the curve says only 63%. Wide windows break the tangent.

Units check, worked example #2 (the quiz, worked)

A reaction has rate q = 4 per second and we look at dt = 0.01 s. The product q·dt has units (1/s)·(s) = dimensionless — good, a probability must be unitless. Numerically 4 × 0.01 = 0.0400 (true value 1 − e−0.04 = 0.0392), a perfectly valid small probability. But the rate 4 itself is not a probability: it has units and it exceeds 1. The window dt = 0.01 s is so short that q·dt = 0.04 << 1, so the tangent and the curve agree to three decimals — the opposite regime from example #1.

python
# jump_prob two ways: exact 1-exp(-q*dt) vs the linear q*dt, plus the error.
import math

def jump_prob(q, dt):
    exact = 1.0 - math.exp(-q * dt)   # the true curve at this window
    linear = q * dt                  # tangent at the origin, slope q
    return exact, linear, abs(exact - linear)

q = 0.2
for dt in [0.1, 2.0]:
    e, lin, err = jump_prob(q, dt)
    print(f"dt={dt}: linear={lin:.4f}  true={e:.4f}  err={err:.4f}")
# dt=0.1: linear=0.0200  true=0.0198  err=0.0002   (great, 1% off)
# dt=2.0: linear=0.4000  true=0.3297  err=0.0703   (linear overshoots badly)

# Watch the error die as O(dt^2) when dt shrinks, and compare to the
# predicted leading term (q*dt)^2 / 2 from the Taylor expansion:
for dt in [0.4, 0.2, 0.1, 0.05]:
    _, _, err = jump_prob(q, dt)
    predicted = (q * dt) ** 2 / 2          # leading error term
    print(f"dt={dt}: err={err:.6f}  predicted={predicted:.6f}")
# dt=0.4: err=0.003228 predicted=0.003200   (halving dt quarters BOTH)
# dt=0.2: err=0.000797 predicted=0.000800
# dt=0.1: err=0.000200 predicted=0.000200   <- matches the worked 0.0002

# Dimensional sanity: q*dt is unitless; q alone is not a probability.
for q_, dt_ in [(0.5, 2.0), (4.0, 0.01)]:
    e, lin, _ = jump_prob(q_, dt_)
    flag = "OVER 1!" if lin > 1 else "ok"
    print(f"q={q_} dt={dt_}: linear={lin:.4f} ({flag})  true={e:.4f}")
# q=0.5 dt=2.0: linear=1.0000 (OVER 1!)  true=0.6321  <- wide window breaks tangent
# q=4.0 dt=0.01: linear=0.0400 (ok)       true=0.0392  <- short window, agree

The library one-liner you would actually use is just the exact form — p = 1 - math.exp(-q*dt) (or 1 - np.exp(-q*dt) for arrays). There is no fancier function to call; the whole content is “use the curve, not the tangent.”

Reconciling with Chapter 0’s tick

Chapter 0 ran a discrete chain on a grid of width dt and flipped a coin with probability pstep each tick. We can now say exactly what that coin probability must be to match a rate-q process: pstep = 1 − e−q·dt, the true per-window jump probability. Many textbooks instead write the approximation pstep ≈ q·dt, and Chapter 0’s “wasted steps” story is exactly the cost of taking dt small enough that the approximation is safe.

Run the numbers for q = 0.2/hr at the coarse tick dt = 1 hr. The honest coin is 1 − e−0.2 = 1 − 0.8187 = 0.1813, while the lazy linear coin is 0.2. They differ by ~0.019 — a 10% relative error baked into the simulation at every single tick. Now take dt = 0.001 hr: honest coin 1 − e−0.0002 = 0.00019998 vs linear 0.0002, agreeing to seven decimals. This is precisely the accuracy-versus-waste trade Chapter 0 made tangible — and the whole reason we will replace the tick entirely with exp(Qt).

So the rate is not a new, exotic object — it is the dt→0 limit of “per-tick jump probability divided by dt.” Take any tick probability, divide by the tick width, shrink the tick: the ratio converges to q. Rate is the per-tick probability with the tick taken out.

From-scratch: sampling the jump time and recovering q

The cleanest proof that q really governs the process: simulate jump times, then recover q from them in two independent ways. We sample the time-to-jump by the inverse-CDF trick (solve U = 1 − e−qt for t, giving t = −ln(1−U)/q), then check that the sample mean is 1/q and that the empirical hazard is flat at q.

python
# Sample jump times for a rate-q process and recover q two ways.
import math, random

q = 0.2
def sample_jump_time(q):
    U = random.random()              # uniform in (0,1)
    return -math.log(1.0 - U) / q   # invert 1 - e^(-q t) = U

times = [sample_jump_time(q) for _ in range(200000)]
mean_t = sum(times) / len(times)
print(f"mean jump time = {mean_t:.3f}  (theory 1/q = {1/q})")   # ~5.0

# Empirical hazard h(t) = (jumps in [t, t+w]) / (still-waiting at t) / w, flat at q.
w = 0.5
for t0 in [0.0, 2.0, 5.0, 10.0]:
    alive = [x for x in times if x >= t0]
    jumped = [x for x in alive if x < t0 + w]
    h = (len(jumped) / len(alive)) / w if alive else 0
    print(f"hazard near t={t0}: {h:.3f}  (theory q = {q})")   # all ~0.2, FLAT

The flat hazard table is the punchline: at t = 0, 2, 5, 10 the empirical per-instant jump-rate is the same ~0.2. The clock does not age. That flatness is the constant hazard, and it is what we cash in for memorylessness in the next chapter.

Side-by-side: rate vs probability, at a glance

Read every row as the same warning: a rate is a slope, a probability is an accumulated area, and you convert between them with a time window. Skip the window and you are comparing a speed to a distance.

From-scratch: the honest discrete-tick coin (matching Chapter 0)

To close the loop with Chapter 0, here is the discrete simulator done correctly — the per-tick coin uses the exact 1 − e−q·dt, not the linear q·dt. We then verify that the simulated mean time-to-jump matches the continuous theory 1/q, and that shrinking dt does not change the answer (only the wasted-step count).

python
# Discrete tick simulator with the EXACT per-tick coin; recover 1/q.
import math, random

def time_to_jump_discrete(q, dt):
    p_step = 1.0 - math.exp(-q * dt)   # honest per-tick jump prob
    t, steps = 0.0, 0
    while random.random() > p_step:    # keep ticking until the coin fires
        t += dt; steps += 1
    return t + dt, steps          # jump lands at end of the firing tick

q = 0.2
for dt in [1.0, 0.1, 0.01]:
    N = 50000
    res = [time_to_jump_discrete(q, dt) for _ in range(N)]
    mean_t = sum(t for t, _ in res) / N
    mean_steps = sum(s for _, s in res) / N
    print(f"dt={dt}: mean_time={mean_t:.2f} (theory {1/q})  wasted_steps/run~{mean_steps:.0f}")
# dt=1.0:  mean_time~5.0  wasted_steps/run~4    (coarse: cheap, but blocky times)
# dt=0.1:  mean_time~5.0  wasted_steps/run~49   (finer: same answer, 10x more steps)
# dt=0.01: mean_time~5.0  wasted_steps/run~499  (fine: same answer, 100x the waste)

The mean time is ~5.0 in every row because the honest coin encodes the true rate; only the wasted step count grows as dt shrinks. That is Chapter 0’s dilemma in one experiment, and it is exactly why the continuous-time machinery (no tick at all) wins. The rate q is the one number all three rows share.

The cumulative hazard, and why it adds up

One last lens that pays off later. Define the cumulative hazard Λ(t) = ∫0t h(u) du — the running total of instantaneous risk. For our constant hazard h(u) = q, this integral is trivially Λ(t) = q·t, a straight line. The survival probability is then S(t) = e−Λ(t) = e−qt — the exponential, recovered from the accumulated hazard.

This view generalizes beautifully. If the hazard were not constant (some other process), you would still get S(t) = e−Λ(t), just with a bent Λ. The CTMC is the special, blessed case where Λ is a straight line of slope q. And when several independent risks act at once (Chapter 3’s race), their cumulative hazards simply add: Λtotal(t) = (q1 + q2 + …)·t, so the combined survival is e−(Σq)t — an exponential with the summed rate. That single fact is the “min of exponentials is exponential” result we will prove geometrically in Chapter 3, here falling out of one addition.

So rates add, and adding rates means multiplying survival probabilities — the logarithm turns the product of e−qit factors into the sum of exponents. Keep that “rates add” reflex; it is the engine behind the generator matrix Q in Chapter 5, whose off-diagonal entries are exactly these additive rates.

From-scratch: the cumulative hazard, numerically

python
# Cumulative hazard Lambda(t) = integral of h(u); for constant h=q it's q*t.
# Recover S(t) = exp(-Lambda(t)) and confirm it matches 1 - F(t).
import math

q = 0.2
def cumulative_hazard(q, t, n=10000):
    du = t / n
    return sum(q for _ in range(n)) * du   # integral of the FLAT hazard q

for t in [1.0, 3.0, 5.0]:
    Lam = cumulative_hazard(q, t)
    S = math.exp(-Lam)
    print(f"t={t}: Lambda={Lam:.3f} (=q*t={q*t})  S=exp(-Lambda)={S:.4f}")
# t=3.0: Lambda=0.600 (=q*t=0.6)  S=exp(-Lambda)=0.5488  -> the Ch2 survival number

# Two independent risks: cumulative hazards ADD, survivals MULTIPLY.
q1, q2 = 0.2, 0.05
t = 4.0
combined = math.exp(-(q1 + q2) * t)
product = math.exp(-q1 * t) * math.exp(-q2 * t)
print(combined, product)   # identical -> sum-of-rates = product-of-survivals (Ch3)

The two printed numbers are identical because adding the rates in the exponent is the same as multiplying the separate survival factors. This tiny experiment is the seed of Chapter 3’s race and Chapter 5’s additive generator — the same “rates add” arithmetic, scaled up.

Recap: the four things a rate is

Before moving on, pin down the four equivalent readings of a single rate q, because the rest of the lesson leans on all of them. (1) Slope: q is the slope of the jump-probability curve at the origin — the tangent we drew. (2) Reciprocal mean: 1/q is the average wait, the mean holding time (Chapter 2). (3) Hazard: q is the constant per-instant jump-chance, the flat hazard we derived. (4) Generator entry: q is an off-diagonal of the matrix Q (Chapter 5).

These are not four facts to memorize separately — they are one number seen from four angles, and being able to flip between them on demand is most of the fluency this lesson is after. When Chapter 5 stacks rates into Q, you should hear all four meanings at once in every entry.

A final sanity anchor for the whole chapter: at q = 0.2/hr, the slope is 0.2, the mean wait is 5 hr, the hazard is a flat 0.2, and the generator entry will be 0.2. One process, one number, four faces — and the curve 1 − e−qt is the bridge that ties the slope at the origin to the full finite-window probability. Everything downstream is this idea wearing more dimensions.

If you internalize only one sentence from this chapter, make it this: multiply a rate by a small time to get a probability, but never confuse the two. The rate lives outside any window and can be huge; the probability lives on a window and is capped at 1. Chapter 2 takes this single rate and asks how long you actually wait before the jump it predicts finally happens.

Those three relations — the tangent, the exact curve, and the reciprocal mean — are the complete toolkit for a single rate. The rest of the lesson is what happens when many rates act at once.

A closing caution worth repeating: the linear rule q·dt is a tool for small windows only, and its failure for large windows is not a rounding issue but a structural one — it can exceed 1 and become a non-probability. Whenever a window is not small compared to 1/q, reach for the exact 1 − e−qt. The discipline of knowing which regime you are in is half of using rates correctly.

python
# Rule of thumb: the linear rule is safe while q*dt < ~0.1 (error < ~0.5%).
import math
for qdt in [0.01, 0.1, 0.5, 1.0]:
    true = 1 - math.exp(-qdt)
    rel_err = abs(qdt - true) / true
    print(f"q*dt={qdt}: linear={qdt}  true={true:.4f}  rel_err={rel_err*100:.1f}%")
# q*dt=0.01: rel_err=0.5%   (linear fine)
# q*dt=0.10: rel_err=5.2%   (borderline)
# q*dt=1.00: rel_err=58%    (linear useless -- use the curve)

The relative error of the linear rule grows from 0.5% at q·dt = 0.01 to 58% at q·dt = 1 — a clean rule of thumb: trust the tangent only while q·dt stays well under ~0.1, and switch to 1 − e−q·dt otherwise.

Chapter 2: Exponential Holding Times & Memorylessness

Chapter 1 gave us a rate q — the constant per-instant tendency to jump. Now ask the natural follow-up: how long do we actually sit in a state before jumping? That waiting time is a random variable called the holding time, and for a CTMC it has one specific, beautiful distribution.

Start with intuition that feels wrong but is true. You have waited 10 minutes for a bus that comes “on average every 15 minutes.” How much longer? Intuition screams “almost here!” The honest answer for this kind of process is: still 15 minutes on average. The clock has no memory of the 10 minutes you already spent.

That bizarre, useful fact is the memoryless property, and the only continuous distribution that has it is the exponential distribution. Its density is λ·e−λt, where λ is the state’s total outflow rate (for a state with a single exit, λ is just that exit’s rate q).

Why exponential and not something else? Because Chapter 1’s constant hazard forces it. If the chance of jumping in the next dt is always λ·dt regardless of how long you have waited, then the survival probability S(t) = P(still here at t) obeys dS/dt = −λ·S — the same scalar ODE from Chapter 1 — whose solution is S(t) = e−λt. The holding time being exponential is not an assumption bolted on; it is the unavoidable consequence of a constant rate.

The mean holding time is 1/λ. This is the most useful single number in the chapter: a state you leave at rate 0.2/hr you sit in for about 5 hours on average; a state you leave at rate 4/s you sit in for a quarter of a second. High rate, short stay; low rate, long stay.

And here is the memoryless property written out. Given that you have already survived a time s in the state, the distribution of the remaining time is identical to the original distribution from zero — same shape, same mean 1/λ. The time already spent is genuinely forgotten. This is exactly the property that makes the future depend only on the present state (the Markov property), which is why CTMCs are built on it.

Deriving the mean holding time = 1/λ

We keep asserting “the mean is 1/λ.” Let us derive it, because the derivation reveals why a faster rate means a shorter stay. The mean of a positive random variable can be computed as the integral of its survival function: E[T] = ∫0∞ S(t) dt = ∫0∞ e−λt dt.

That integral is elementary. The antiderivative of e−λt is −(1/λ)e−λt. Evaluate from 0 to ∞: at the top e−∞ = 0, at the bottom e0 = 1, giving E[T] = 0 − (−(1/λ)·1) = 1/λ. The mean is the reciprocal of the rate, full stop.

The intuition the integral makes precise: S(t) = e−λt is the “fraction still waiting at time t,” and the area under it is the average wait. A bigger λ makes that curve plunge faster, shrinking the area, shrinking the mean. High rate, small area, short stay — geometry, not jargon.

The spread tells the same story. The variance of an exponential is 1/λ², so the standard deviation is 1/λ — equal to the mean. Exponential waits are highly variable: the typical wait is about as big as the average wait itself, which is why sample holding times scatter so widely in the histogram below.

Memorylessness, proven in one line

The memoryless property is the equation P(T > s + r | T > s) = P(T > r): given you have already waited s, the chance of waiting r more equals the chance a fresh clock waits r. Prove it directly from the definition of conditional probability.

The exponents subtract: e−λ(s+r)+λs = e−λr = P(T > r). The s cancels completely — the time already spent leaves no trace. That single cancellation, e−λ(s+r) / e−λs = e−λr, is the whole memoryless property, and it works only because the exponent is linear in t (i.e. constant hazard). No other continuous distribution has this.

Worked example: UP state, λ = 0.2/hr

The machine is UP and leaves at total rate λ = 0.2 per hour. Three questions, fully worked.

Mean stay. 1/λ = 1/0.2 = 5.0 hours. On average the machine stays up 5 hours before something happens.

Survival to t = 3. P(still up at 3) = e−λ·3 = e−0.2×3 = e−0.6. Numerically e−0.6 = 0.5488. So there is a 54.88% chance it is still up after 3 hours.

Memoryless check. Given it has already been up 3 hours, what is the chance it survives 2 more? Naively you might discount for “age,” but the formula is P(survive 2 more | survived 3) = e−λ·2 = e−0.4 = 0.6703. Now compare: a brand-new machine surviving its first 2 hours has probability e−0.4 = 0.6703 too. Identical. The 3 hours already spent are completely forgotten.

That equality — 0.6703 = 0.6703 — is the memoryless property in two numbers. The conditional remaining-time distribution does not depend on s, so its mean is always 1/λ = 5 no matter how long you have waited.

Worked example: spread, not just center

Same UP state, λ = 0.2/hr. The variance is 1/λ² = 1/0.04 = 25 hr², so the standard deviation is √25 = 5 hr — identical to the mean of 5 hr. Concretely, this means a state with a 5-hour average stay routinely produces stays of 1 hour and stays of 12 hours; the wait is not “about 5 hours,” it is “5 hours give or take 5 hours.”

The median tells the same story of skew. Solve 1 − e−λt = 0.5 for t: e−λt = 0.5, so t = ln 2 / λ = 0.6931/0.2 = 3.47 hr. The median (3.47 hr) sits well below the mean (5 hr) — a long right tail drags the average up. Half the holding times finish before 3.47 hr, but the rare long ones pull the mean to 5.

Why “average every 15 minutes” felt wrong (the bus paradox)

Return to the bus. If buses arrive as an exponential process with mean gap 15 min, and you show up at a random instant, the time until the next bus has mean 15 min — not 7.5. Memorylessness says the process does not care that some time has elapsed since the last bus; from your arrival the clock restarts fresh. The 10 minutes you waited bought you nothing because the wait never “ages.”

This is genuinely counterintuitive and worth sitting with: a deterministic 15-minute schedule would give a 7.5-minute expected wait, but pure exponential randomness doubles it back to 15. Constant hazard is exactly the regime where patience is never rewarded.

python
# Sample exponential holding times by inverse-CDF: t = -ln(U)/lambda, U ~ Uniform(0,1).
import math, random

def exp_sample(lam):
    U = random.random()          # uniform in (0,1)
    return -math.log(U) / lam     # invert 1-e^(-lam t) = U

lam = 0.2
samples = [exp_sample(lam) for _ in range(100000)]
mean = sum(samples) / len(samples)
print(f"sample mean = {mean:.3f}  (theory 1/lam = {1/lam})")   # ~5.0

# Memoryless check: keep only holds that survived past s=3, subtract 3, re-average.
s = 3.0
residual = [x - s for x in samples if x > s]
print(f"residual mean = {sum(residual)/len(residual):.3f}")   # still ~5.0, unchanged!

# Variance and median, both from the sample, both matching theory.
mu = mean
var = sum((x - mu) ** 2 for x in samples) / len(samples)
samples.sort(); median = samples[len(samples) // 2]
print(f"var={var:.2f} (theory 1/lam^2={1/lam**2})  std={var**0.5:.2f} (=mean)")
print(f"median={median:.2f} (theory ln2/lam={math.log(2)/lam:.2f})")   # ~3.47 < mean

# Survival-function integral: E[T] = integral_0^inf S(t) dt, done numerically.
def survival_integral(lam, T=100.0, n=100000):
    dt = T / n
    return sum(math.exp(-lam * (i + 0.5) * dt) for i in range(n)) * dt
print(f"area under S(t) = {survival_integral(lam):.3f}  (= 1/lam = {1/lam})")  # ~5.0

The library one-liner: numpy.random.exponential(scale=1/lam, size=N) — note NumPy parameterizes by the scale 1/λ (the mean), not the rate, a classic gotcha.

Exponential vs the alternatives, at a glance

The exponential is the only row with a flat hazard, and the flat hazard is exactly what lets the future depend on the present state alone — the Markov property. Every other model needs to remember how long it has waited, which means the “state” must secretly include a clock. CTMCs buy their simplicity by paying with the exponential assumption.

From-scratch: the bus paradox, simulated

Let us make the “still 15 minutes” claim undeniable by simulation. We generate a long stream of exponential inter-arrival gaps (mean 15), drop a random observer into the timeline, and measure the wait until the next arrival. If memorylessness holds, that wait should average ~15, not ~7.5.

python
# Bus paradox: random observer's wait-for-next under exponential gaps.
import math, random

mean_gap = 15.0; lam = 1.0 / mean_gap
# build a long timeline of arrival instants
arrivals, t = [], 0.0
while t < 1_000_000:
    t += -math.log(random.random()) / lam
    arrivals.append(t)

# drop N observers at uniformly random times, measure wait to next arrival
import bisect
waits = []
for _ in range(200000):
    obs = random.uniform(0, arrivals[-2])
    idx = bisect.bisect_right(arrivals, obs)   # first arrival after obs
    waits.append(arrivals[idx] - obs)
print(f"observer wait = {sum(waits)/len(waits):.2f} min  (naive guess 7.5, truth ~15)")
# observer wait ~ 15.0  -> the wait did NOT halve; memorylessness in the wild

# Contrast: deterministic 15-min schedule gives exactly 7.5 (linear ramp).
det = [(15 - (x % 15)) for x in (random.uniform(0, 15000) for _ in range(200000))]
print(f"deterministic wait = {sum(det)/len(det):.2f} min  (~7.5)")

The two numbers — ~15 for exponential, ~7.5 for the fixed schedule — are the bus paradox quantified. Randomness with constant hazard doubles your expected wait relative to a clockwork schedule, because the big gaps are both longer and more likely to contain your random arrival. That is the price (and the simplicity) of memorylessness.

The flip side: exponential gaps make a Poisson process

If holding times are exponential and independent, the count of jumps in a fixed window is Poisson-distributed. The two views are the same process seen from different axes: exponential measures “time between events,” Poisson measures “number of events in a window.” A rate-λ exponential gap process produces, over a window of length T, a Poisson count with mean λT.

This is why CTMCs and the Poisson process are inseparable. Each state’s exit clock is a Poisson alarm; the whole chain is a tangle of competing Poisson processes (Chapter 3’s race). When Chapter 5’s uniformization reinterprets the chain as ticking at a fast Poisson rate μ, it is cashing in exactly this exponential⇔Poisson duality. The exponential holding time is not just convenient — it is the atom of the entire continuous-time-Markov world.

From-scratch: histogram convergence to the exponential density

python
# Bin sampled holding times; the histogram should trace lam*exp(-lam t).
import math, random

lam = 0.2; N = 200000; binw = 2.0
hist = {}
for _ in range(N):
    t = -math.log(random.random()) / lam      # inverse-CDF sample
    b = int(t / binw)
    hist[b] = hist.get(b, 0) + 1

print(f"{'t':>6} {'empirical':>10} {'theory':>10}")
for b in range(6):
    t_mid = (b + 0.5) * binw
    emp = hist.get(b, 0) / N / binw            # empirical density
    theory = lam * math.exp(-lam * t_mid)        # lam e^(-lam t)
    print(f"{t_mid:>6.1f} {emp:>10.4f} {theory:>10.4f}")
#      t  empirical    theory
#    1.0     0.1810     0.1810   (peak near 0, density = lam at t=0)
#    5.0     0.0809     0.0809   (decaying exponentially)
#    9.0     0.0362     0.0362   (long right tail -> mean > median)

The empirical and theoretical columns match across the whole range: the histogram is the exponential density λe−λt, peaking at λ when t = 0 and decaying with the long right tail that drags the mean above the median. Sampling confirms the formula, and the formula explains the shape.

Worked example: reading the CDF, survival, and quantiles

One rate, four questions, all from the same curve at λ = 0.2/hr. CDF (jumped by t=4): 1 − e−0.8 = 1 − 0.4493 = 0.5507 — just past even odds. Survival (still here at t=4): the complement, 0.4493. Median: ln2/λ = 3.47 hr (solved earlier). 90th percentile: solve 1 − e−λt = 0.9 → t = −ln(0.1)/0.2 = 2.3026/0.2 = 11.5 hr.

That last number is the eye-opener: 90% of holding times finish within 11.5 hr, but the long exponential tail means the remaining 10% can run far longer, which is what pulls the mean (5 hr) well above the median (3.47 hr). A single rate fixes the entire distribution — every quantile, the mean, the spread — with no extra parameters.

python
# All the exponential summaries from one rate: CDF, survival, quantiles.
import math

lam = 0.2
cdf      = lambda t: 1 - math.exp(-lam * t)        # P(jumped by t)
survival = lambda t: math.exp(-lam * t)            # P(still here at t)
quantile = lambda p: -math.log(1 - p) / lam      # inverse CDF

print(f"CDF(4)      = {cdf(4):.4f}")        # 0.5507
print(f"survival(4) = {survival(4):.4f}")   # 0.4493
print(f"mean        = {1/lam:.3f}")         # 5.000
print(f"median      = {quantile(0.5):.3f}")  # 3.466
print(f"90th pct    = {quantile(0.9):.3f}")  # 11.513
# median < mean < 90th-pct: the long right tail in three numbers.

The three summary statistics — median 3.47, mean 5.0, 90th percentile 11.5 — spread out exactly as a right-skewed distribution should, all generated by the single rate 0.2. This is why a CTMC needs only rates: the rate is the full holding-time distribution.

Worked example: comparing two states’ clocks

Put the UP and DOWN clocks side by side to feel the rate→dwell relationship. UP leaves at λ = 0.2/hr: mean dwell 1/0.2 = 5 hr, and after 5 hr the survival is e−1 = 0.368 (37% still up). DOWN leaves at λ = 1.0/hr: mean dwell 1/1.0 = 1 hr, and after 1 hr survival is also e−1 = 0.368.

The shared 0.368 is no accident: every exponential has the same survival e−1 at exactly one mean-lifetime, regardless of rate. It is the “1/e” landmark of exponential decay — after one mean dwell, a fraction 1/e ≈ 37% of the original mass remains. A fast state and a slow state look identical when you measure time in units of their own mean dwell. Rate sets the timescale; the shape is universal.

python
# Two states, different rates, same shape when scaled by their mean dwell.
import math

for name, lam in [("UP", 0.2), ("DOWN", 1.0)]:
    mean = 1 / lam
    surv_at_mean = math.exp(-lam * mean)        # always e^-1
    surv_at_2mean = math.exp(-lam * 2 * mean)    # always e^-2
    print(f"{name:4s} lam={lam}: mean={mean:.1f}h  S(mean)={surv_at_mean:.3f}  S(2*mean)={surv_at_2mean:.3f}")
# UP   lam=0.2: mean=5.0h  S(mean)=0.368  S(2*mean)=0.135
# DOWN lam=1.0: mean=1.0h  S(mean)=0.368  S(2*mean)=0.135   <- SAME numbers!
# Measured in its own mean-dwell units, every exponential is identical.

Both states show S(mean) = 0.368 and S(2·mean) = 0.135 — the rate only rescales the time axis, never the shape. This universality is why a single rate fully specifies a holding time: once you know the timescale 1/λ, the whole curve is determined.

Why memorylessness is the load-bearing assumption

It is worth being explicit about how much rides on the exponential holding time. The entire Markov property of a CTMC — “the future depends only on the current state” — requires memorylessness. If holding times had memory (if a state “aged”), then knowing the current state would not be enough; you would also need to know how long you had been there, and the clean matrix machinery of Chapters 5–8 would collapse.

This is the price of admission. By assuming exponential holding times we buy the ability to summarize the whole process in one matrix Q, to compute P(t) in closed form, and to find π from a linear solve. Drop the assumption and none of that survives intact. So when you model a real system as a CTMC, you are asserting that its states have constant hazard — a claim to check, not take for granted.

When the claim fails — machines that wear out, customers who grow impatient — the fix is to enrich the state space until the future depends only on the (enlarged) present. Add an “age” or “phase” dimension and the process becomes Markov again, at the cost of more states. This phase-type trick is how practitioners model non-exponential waits while keeping the CTMC toolkit; it is the standard escape hatch, and it all rests on the memorylessness we derived here.

python
# Direct memoryless test: residual wait after surviving s is independent of s.
import math, random

lam = 0.2; N = 500000
samples = [-math.log(random.random()) / lam for _ in range(N)]

for s in [0.0, 3.0, 8.0, 15.0]:
    residuals = [x - s for x in samples if x > s]   # wait remaining, given survived s
    mean_res = sum(residuals) / len(residuals)
    print(f"survived s={s:4.1f}: mean residual wait = {mean_res:.3f}  (always ~1/lam = {1/lam})")
# survived s= 0.0: mean residual wait ~ 5.0
# survived s= 8.0: mean residual wait ~ 5.0   <- the 8 hours waited are FORGOTTEN
# survived s=15.0: mean residual wait ~ 5.0   <- still 5, no matter how long

# Contrast a NON-memoryless wait (uniform on [0,10]): residual DOES shrink.
unif = [random.uniform(0, 10) for _ in range(N)]
for s in [0.0, 8.0]:
    res = [x - s for x in unif if x > s]
    print(f"uniform, survived s={s}: mean residual = {sum(res)/len(res):.2f}  (shrinks! not memoryless)")
# uniform survived 0: ~5.0 ; survived 8: ~1.0  -> waiting DID use up the clock
# Only the exponential keeps a CONSTANT residual -- that's why CTMCs use it.
# A uniform "holding time" would force the state to secretly remember its age,
# breaking the Markov property; the exponential is the unique memoryless choice.
print("exp residual is flat; uniform residual shrinks toward 0 as s grows.")
print("=> only the exponential lets the CURRENT STATE alone predict the future.")
# That predictive sufficiency of the present state IS the Markov property.
# It is the one assumption a CTMC keeps; the discretization grid is the one it drops.
# When a real wait isn't memoryless, enlarge the state (phase-type) to restore it.

The mean residual wait is ~5 hr regardless of whether you have already survived 0, 3, 8, or 15 hours — the elapsed time is genuinely forgotten. This is the memoryless property measured directly from samples, and it is the empirical face of the constant hazard we derived: waiting never makes the next jump any closer.

Chapter 3: Competing Clocks — The Race

Chapter 2 handled a state with a single exit. But a real state can jump to several places. From UP, the machine can FAIL (rate 0.2/hr) or be pulled for scheduled MAINTENANCE (rate 0.05/hr). How does a CTMC decide both when it leaves and where it goes?

The central image of the whole subject: one exponential alarm clock per destination. Set clock FAIL to a random time drawn from exp(0.2) and clock MAINT to a random time drawn from exp(0.05). They tick down simultaneously. Whichever rings first wins — and that single event decides everything.

The winning clock’s ring time is the holding time (the when). The identity of the winning clock is the destination (the where). You do not need a separate “which event happened” coin: the clocks themselves choose, by racing. One mechanism answers both questions at once.

Two facts fall out of this race, and they are the heart of the chapter. Fact 1 (the when): the minimum of independent exponentials is itself exponential, with rate equal to the sum of the individual rates. So the holding time in UP is exp(0.2 + 0.05) = exp(0.25). Fact 2 (the where): the probability that a particular clock wins is its rate divided by the total rate, qi / Σq.

Both facts have clean intuitions. The when is faster than any single clock because more clocks racing means someone rings sooner — two alarms beat one. The where is proportional to rate because a faster clock draws smaller times more often, so it wins more often, in exact proportion to its rate’s share of the total.

Notice what this buys us: a multi-exit state is fully described by its outgoing rates. The total rate sets the clock; the ratios set the destination. We will package exactly these two pieces, separately, in Chapter 4, and stack all the rates into one matrix in Chapter 5.

Proving Fact 1: the min of exponentials is exponential(sum)

Why is the holding time exponential with the summed rate? Work with survival probabilities, because for a minimum they multiply. The minimum of the clocks exceeds t exactly when every clock exceeds t. For independent clocks, that joint probability is the product of the individual survivals.

The exponents add (Chapter 1’s cumulative-hazard fact). And e−(q1+q2)t is precisely the survival function of an exponential with rate q1+q2. So the minimum is exponential with the summed rate — no approximation, exact. For our numbers, min(exp(0.2), exp(0.05)) is exp(0.25), mean 1/0.25 = 4 hr.

Proving Fact 2: P(clock i wins) = q_i/Σq

Why does the winning probability equal the rate share? Condition on the winning instant. Clock i wins at time t if it fires in [t, t+dt] — probability qi·e−qit·dt — while every other clock has survived past t. Integrate over all possible winning times t.

The product of all the survivals collapses to e−(Σq)t, the integral of that is 1/Σq, and the leftover qi out front gives qi/Σq. Plug in: P(fail) = 0.2/0.25 = 0.80. The faster clock wins more often, in exact proportion to its rate — derived, not asserted.

Notice the two facts came from the same collapse ∏ e−qjt = e−(Σq)t. The when reads it as a survival curve; the where integrates it against one firing clock. One identity, two answers — which is exactly why a single race resolves both questions.

Worked example: FAIL vs MAINTENANCE from UP

From UP, the rates are qfail = 0.2 and qmaint = 0.05 per hour. Total leave-rate Σ = 0.2 + 0.05 = 0.25 per hour.

The when. The holding time is exponential with rate Σ = 0.25, so its mean is 1/Σ = 1/0.25 = 4.0 hours. Note this is shorter than the 5 hours you would get from fail alone (1/0.2 = 5) — the second clock pulls the expected exit earlier.

The where. Probability the jump goes to FAIL: qfail/Σ = 0.2/0.25 = 0.80. Probability it goes to MAINT: qmaint/Σ = 0.05/0.25 = 0.20. They sum to 1.00, as a destination distribution must.

So in plain words: from UP we typically sit about 4 hours, then 80% of the time we have failed and 20% of the time we are in maintenance. One race, two answers.

The quiz, worked

A state has outgoing rates 1, 2, 3 to A, B, C. Total Σ = 1 + 2 + 3 = 6. Then P(C) = 3/6 = 0.5, P(B) = 2/6 = 0.3333, P(A) = 1/6 = 0.1667. Mean holding time = 1/6 = 0.1667 time units — faster than any single clock because all three race together.

python
# A race: draw one exponential time per destination; min is WHEN, argmin is WHERE.
import math, random

def race(rates):
    times = [-math.log(random.random()) / q for q in rates]  # one clock each
    t_win = min(times)                 # the holding time (when)
    i_win = times.index(t_win)        # the destination index (where)
    return i_win, t_win

rates = [0.2, 0.05]   # fail, maint
N = 100000
wins = [0, 0]; total_t = 0.0
for _ in range(N):
    i, t = race(rates)
    wins[i] += 1; total_t += t
print(f"P(fail) emp = {wins[0]/N:.3f}  theory = {0.2/0.25:.3f}")  # ~0.80
print(f"mean hold emp = {total_t/N:.3f}  theory = {1/0.25:.3f}")    # ~4.00

The library one-liner draws all the clocks at once and uses argmin/min directly: times = np.random.exponential(1/np.array(rates)); i_win, t_win = times.argmin(), times.min().

From-scratch: verifying both facts at once, three clocks

Let us stress the result with three competing rates and confirm both the winner shares and the mean hold against theory in one experiment. We use the quiz’s rates 1, 2, 3 (total 6), expecting P(A,B,C) = (1,2,3)/6 and mean hold 1/6 = 0.1667.

python
# Three-clock race: confirm P(win) = rate/sum AND mean hold = 1/sum.
import math, random

def race(rates):
    times = [-math.log(random.random()) / q for q in rates]
    t_win = min(times)
    return times.index(t_win), t_win

rates = [1.0, 2.0, 3.0]; S = sum(rates)
N = 300000
wins = [0, 0, 0]; total_t = 0.0
for _ in range(N):
    i, t = race(rates); wins[i] += 1; total_t += t
for i, q in enumerate(rates):
    print(f"clock {i}: emp P(win)={wins[i]/N:.3f}  theory={q/S:.3f}")
print(f"mean hold emp={total_t/N:.4f}  theory 1/sum={1/S:.4f}")
# clock 0: emp P(win)~0.167  theory 0.167   (rate 1 / sum 6)
# clock 1: emp P(win)~0.333  theory 0.333   (rate 2 / sum 6)
# clock 2: emp P(win)~0.500  theory 0.500   (rate 3 / sum 6)
# mean hold emp~0.1667  theory 0.1667       (faster than ANY single clock)

# Independence check: the winning TIME does not depend on WHO won.
import statistics
by_winner = {0: [], 1: [], 2: []}
for _ in range(150000):
    i, t = race(rates); by_winner[i].append(t)
for i in by_winner:
    print(f"mean hold | winner={i}: {statistics.mean(by_winner[i]):.3f}")
# all ~0.167 regardless of winner -> WHEN and WHERE are independent (Chapter 4)

The last block is the quiet bombshell that powers Chapter 4: the winning time is the same ~0.167 no matter which clock won. The when and the where are statistically independent, which is precisely why you may sample them separately and in either order.

Worked example: adding a clock speeds up the exit

Concretely watch the “more clocks ring sooner” effect. With fail alone (q = 0.2), the mean exit time is 1/0.2 = 5 hr. Add the maintenance clock (q = 0.05): the total rate becomes 0.25 and the mean exit drops to 1/0.25 = 4 hr. Add a third clock, say a reboot at q = 0.15: total 0.4, mean exit 1/0.4 = 2.5 hr. Each new exit, however slow, can only shorten the expected stay, because the minimum of more clocks is smaller.

And the destination split rebalances each time. With all three (fail 0.2, maint 0.05, reboot 0.15, total 0.4): P(fail) = 0.2/0.4 = 0.50, P(maint) = 0.05/0.4 = 0.125, P(reboot) = 0.15/0.4 = 0.375. They sum to 1, and adding the reboot clock stole probability mass from the others in proportion to its rate. The race is a zero-sum competition for the shares qi/Σq.

The geometry of the race

Here is a picture that makes both facts obvious. Imagine each clock’s sampled time as a point thrown down on the time axis, with faster clocks throwing points closer to 0 on average. The winner is whichever point is leftmost (the min), and a faster clock’s cloud of points sits further left, so it wins more often — exactly in proportion to its rate. The winning time is the position of that leftmost point, which gets pulled toward 0 as you add more clouds, because more points means a smaller leftmost one.

This is why the two facts are inseparable: both are statements about the minimum of a pile of exponential points. “Who is leftmost” (argmin) gives the destination; “how far left” (min) gives the time. One geometric act, read two ways — the unifying image behind the generator matrix Q we build next.

From-scratch: the min really is exp(sum)

Fact 1 deserves a direct empirical check: draw the minimum of two exponentials thousands of times and confirm its distribution is exponential with the summed rate. We compare the sample mean to 1/Σq and the sample survival curve to e−Σq·t.

python
# Verify: min(exp(q1), exp(q2)) ~ exp(q1+q2), in mean AND in survival.
import math, random

q1, q2 = 0.2, 0.05; S = q1 + q2     # total rate 0.25
def exp_draw(q): return -math.log(random.random()) / q

mins = [min(exp_draw(q1), exp_draw(q2)) for _ in range(300000)]
print(f"mean of min = {sum(mins)/len(mins):.3f}  (theory 1/S = {1/S:.3f})")   # ~4.0

# Survival of the min should match exp(-S t) at several t.
for t in [2.0, 4.0, 8.0]:
    emp = sum(1 for x in mins if x > t) / len(mins)
    print(f"P(min > {t}) emp={emp:.3f}  theory exp(-S t)={math.exp(-S*t):.3f}")
# P(min > 4) emp~0.368  theory 0.368   -> min is exp(0.25), NOT exp(0.2) or exp(0.05)

# And argmin shares: faster clock wins proportionally more.
fail_wins = sum(1 for _ in range(300000) if exp_draw(q1) < exp_draw(q2))
print(f"P(fail wins) = {fail_wins/300000:.3f}  (theory q1/S = {q1/S:.3f})")   # ~0.80

The mean of the minimum is ~4.0 (matching 1/0.25, not 1/0.2 or 1/0.05), the survival curve traces e−0.25t, and the fail clock wins ~80% of races — all three verified facts in one experiment. The minimum of clocks is genuinely a new, faster exponential, and the winner share is genuinely the rate share.

The efficient race: draw the time first, then the winner

The naive race draws one clock per exit, but Facts 1 and 2 give a cheaper recipe that scales to any number of destinations: draw the holding time once from exp(Σq), then pick the winner from the categorical distribution qi/Σq. Two draws total, regardless of how many exits. This is the engine of the Gillespie simulator (Chapter 9), and it works precisely because when and where are independent (Chapter 4).

python
# Efficient race: ONE exp(sum) draw for time, ONE categorical draw for winner.
import math, random

def efficient_race(rates):
    S = sum(rates)
    t_win = -math.log(random.random()) / S      # WHEN: exp(sum of rates)
    u = random.random() * S; cum = 0.0             # WHERE: categorical by rate share
    for i, q in enumerate(rates):
        cum += q
        if u <= cum: return i, t_win
    return len(rates) - 1, t_win

rates = [0.2, 0.05, 0.15]            # fail, maint, reboot (sum 0.4)
N = 300000; wins = [0] * 3; tt = 0.0
for _ in range(N):
    i, t = efficient_race(rates); wins[i] += 1; tt += t
S = sum(rates)
print("win shares:", [round(w / N, 3) for w in wins])      # [0.5 0.125 0.375]
print("theory    :", [round(q / S, 3) for q in rates])    # [0.5 0.125 0.375]
print("mean hold :", round(tt / N, 3), "theory", round(1 / S, 3))  # 2.5

The efficient two-draw race reproduces the same shares [0.5, 0.125, 0.375] and mean hold 2.5 as the naive one-clock-per-exit version — but with constant cost regardless of the number of exits. This is exactly why real simulators never draw a clock per edge; they draw the time and the destination separately, the recipe Chapter 4 formalizes.

The race is the whole continuous-time story

It is worth pausing on how much this one mechanism explains. Every transition in any CTMC, no matter how many states, is a race of exponential clocks — one clock per outgoing edge of the current state. The generator matrix Q (Chapter 5) is just a ledger of all these clock rates; the simulator (Chapter 9) is just “run the race, jump, repeat.” If you understand this single race, you understand the dynamics of the entire subject.

The two facts also resolve a question that trips people up: “does adding a rare exit change anything if it almost never fires?” Yes — even a slow clock shortens the expected holding time (Fact 1, rates add) and steals a sliver of destination probability (Fact 2, shares). A rarely-taken edge still speeds up your average departure, because the minimum of more clocks is always smaller. There is no “free” exit; every rate you add to a state makes it leave sooner.

python
# Adding a rare exit: how much does it shorten the stay and shift the shares?
def race_stats(rates):
    S = sum(rates)
    return 1 / S, [q / S for q in rates]   # mean hold, win shares

# before: just fail (0.2) and maint (0.05)
hold0, shares0 = race_stats([0.2, 0.05])
print(f"2 exits : mean hold={hold0:.3f}  shares={[round(x,3) for x in shares0]}")

# after: add a RARE reboot exit at rate 0.02
hold1, shares1 = race_stats([0.2, 0.05, 0.02])
print(f"3 exits : mean hold={hold1:.3f}  shares={[round(x,3) for x in shares1]}")
print(f"stay shortened by {hold0 - hold1:.3f}h even though reboot fires only {shares1[2]*100:.0f}% of the time")
# 2 exits : mean hold=4.000  shares=[0.8, 0.2]
# 3 exits : mean hold=3.704  shares=[0.741, 0.185, 0.074]
# the rare 0.02 reboot still cut ~0.3h off the average stay -- no free exit
print(f"fail share dropped {shares0[0]:.3f} -> {shares1[0]:.3f} (reboot stole some)")
# even fail's share fell from 0.80 to 0.74 -- the new clock competes for ALL of it
print(f"all 3 shares sum to {sum(shares1):.3f}")   # 1.000 -- still a valid distribution
# takeaway: every exit you add (1) shortens the stay and (2) dilutes every
# other destination's share -- the race is a zero-sum split of probability.
print(f"mean hold: {hold0:.3f}h -> {hold1:.3f}h after adding the rare reboot exit")
# 4.000 -> 3.704: a 0.02 clock you rarely take still trims the stay measurably.
# min of more clocks is always smaller -> no exit is ever truly "free".
# This is the race (min + argmin) that the generator Q books in Chapter 5.

Adding a reboot clock that wins only 7% of races still trims the mean stay from 4.0 to 3.7 hr and shaves probability off every other destination — concrete proof that no exit is “free.” Every rate you add to a state, however rare, makes the state leave sooner and rebalances where it goes.

Chapter 4: Two Clocks in One — Jump Chain vs Holding Time

The race in Chapter 3 told us two things in a single mechanism: when we leave and where we go. It turns out you can ask those two questions separately, with two simpler objects, and lose nothing. This factoring is the cleanest way to understand — and to simulate — a CTMC.

The first object is the embedded jump chain: a plain discrete Markov chain that answers only “where do I go next?” with no notion of time. Its transition probabilities are Pij = qij / λi, each destination’s rate divided by the total leave-rate λi = Σk qik. These are exactly the “where” probabilities from the race — the rates, normalized into a distribution.

The second object is the holding-time clock: for each state, an exponential with rate λi that answers only “how long do I stay?” This is exactly the “when” from the race — the minimum over all the competing clocks, which we proved is exponential with the summed rate.

Here is the claim that earns this chapter its place: knowing the jump chain P and the holding rates λ is exactly the same information as knowing all the rates q_ij. You can recover the rates from them — qij = λi · Pij — and you can recover them from the rates. No information is created or lost in the split; it is a change of coordinates on the same process.

So there are two equivalent ways to simulate a CTMC. The race way (Chapter 3): draw one exponential clock per destination, take the min for the time and the argmin for the destination. The reassemble way (this chapter): first draw the next state from the jump-chain row P_i·, then separately draw the dwell time from exp(λi), then move — alternating “where” and “how long” forever.

These two recipes produce statistically identical sample paths, and the widget below lets you watch the reassemble recipe rebuild the very staircase the race produces. That equivalence is not a coincidence: it is the mathematical fact that the min-and-argmin of independent exponentials factor into an exponential(sum) for the min and an independent categorical(rates/sum) for the argmin. Where and how-long are independent, which is exactly why you may sample them in either order.

Why the factoring is exact, not approximate

The split rests on one fact from Chapter 3’s race: the winning time and the winning identity are statistically independent. We can see this directly from the joint density of “clock i wins at time t,” which was qie−(Σq)t. Factor it:

The first bracket depends only on t (it is the exp(Σq) holding-time density); the second depends only on i (it is the jump-chain probability). A joint density that factors into a function of t times a function of i is the textbook definition of independence. So “when” and “where” carry no information about each other — you lose nothing by drawing them separately.

This is the license for the reassemble recipe. Because the factors are independent, the order of sampling is irrelevant: draw the destination first then the dwell, or the dwell first then the destination — same distribution over trajectories either way. The race draws them jointly (one min, one argmin); the reassemble draws them in sequence. Both are exact samplers of the identical process.

Worked example: state A with q_AB = 1, q_AC = 3

Take a 3-state chain. From A the outgoing rates are qAB = 1 and qAC = 3. First the total leave-rate: λA = qAB + qAC = 1 + 3 = 4.

Jump-chain row from A. PAB = qAB/λA = 1/4 = 0.25 and PAC = qAC/λA = 3/4 = 0.75. They sum to 0.25 + 0.75 = 1.00, a valid distribution, and there is no self-loop: PAA = 0 (a CTMC never “jumps” to where it already is).

Holding time in A. Exponential with rate λA = 4, so mean dwell = 1/4 = 0.25 time units.

Reading it together: from A we typically sit about 0.25 time units, then go to C three times as often as to B. Now reconstruct the race answer as a sanity check: the race from A draws clocks exp(1) and exp(3); their min is exp(4) (mean 0.25, matching the dwell) and P(C wins) = 3/(1+3) = 0.75 (matching P_AC). The two views agree exactly — that is the equivalence, in numbers.

The quiz, worked

State A has qAB = 2, qAC = 2. Then λA = 4, so PAB = 2/4 = 0.5 (a fair coin for direction) and the mean dwell is 1/4 = 0.25 — a faster clock than either single rate (1/2 = 0.5) because both clocks race.

python
# Build the embedded jump matrix P from the generator Q, then simulate by
# alternating: draw next state from P's row, draw dwell from exp(lambda_i).
import numpy as np

def jump_matrix(Q):
    Q = np.asarray(Q, float)
    lam = -np.diag(Q)              # total leave-rate per state = -q_ii
    P = Q / lam[:, None]            # divide each row by its leave-rate
    np.fill_diagonal(P, 0.0)        # no self-jumps: P_ii = 0
    return P, lam

def simulate_reassemble(Q, s0, T):
    P, lam = jump_matrix(Q)
    t, s, path = 0.0, s0, [(0.0, s0)]
    while t < T:
        dwell = np.random.exponential(1.0 / lam[s])   # HOW LONG (clock)
        t += dwell
        s = np.random.choice(len(lam), p=P[s])    # WHERE (jump chain)
        path.append((t, s))
    return path

Q = [[-4, 1, 3], [2, -2, 0], [0, 1, -1]]
P, lam = jump_matrix(Q)
print(P[0])   # [0.   0.25 0.75]  -> P_AB=0.25, P_AC=0.75
print(1/lam[0]) # 0.25  -> mean dwell in A

There is no single library call for “CTMC simulate” in base NumPy; the two-liner above is the idiomatic version. The point is that jump_matrix(Q) — normalize off-diagonals by the row leave-rate, zero the diagonal — is the only conversion you need.

Worked example: a full hand-trajectory on the 3-state chain

Let us walk three jumps by hand on Q = [[−4,1,3],[2,−2,0],[0,1,−1]] (states A,B,C), using the reassemble recipe with prepared random draws. From A: λA = 4, jump row [0, 0.25, 0.75].

Jump 1 (from A). Dwell from exp(4): with U = 0.5, dwell = −ln(0.5)/4 = 0.1733. Clock now 0.1733. Destination with U2 = 0.6: cumulative [B:0.25, C:1.0], and 0.6 > 0.25 → go to C.

Jump 2 (from C). C’s only exit is to B at rate 1, so λC = 1, jump row [0,1,0] (B with certainty). Dwell from exp(1): with U = 0.3, dwell = −ln(0.3)/1 = 1.204. Clock now 0.1733 + 1.204 = 1.377. Destination: B (no choice).

Jump 3 (from B). B exits to A at rate 2, so λB = 2, jump row [1,0,0]. Dwell from exp(2): with U = 0.8, dwell = −ln(0.8)/2 = 0.1116. Clock now 1.377 + 0.1116 = 1.489. Destination: A. Trajectory so far: A →(0.173) C →(1.377) B →(1.489) A — a where-sequence A,C,B,A with dwells 0.173, 1.204, 0.112 stitched on.

From-scratch: confirm reassemble reproduces the race occupancy

python
# Reassemble and race must give the SAME long-run occupancy. Check it.
import numpy as np

Q = np.array([[-4.0, 1.0, 3.0], [2.0, -2.0, 0.0], [0.0, 1.0, -1.0]])
lam = -np.diag(Q)
P = Q / lam[:, None]; np.fill_diagonal(P, 0.0)

def occupancy(use_race, T=200000.0):
    t, s = 0.0, 0; time_in = np.zeros(3)
    while t < T:
        if use_race:                                   # RACE: min/argmin of clocks
            times = np.array([np.random.exponential(1/Q[s,j]) if Q[s,j]>0 else np.inf for j in range(3)])
            dwell, nxt = times.min(), int(times.argmin())
        else:                                          # REASSEMBLE: dwell then jump
            dwell = np.random.exponential(1 / lam[s]); nxt = np.random.choice(3, p=P[s])
        time_in[s] += dwell; t += dwell; s = nxt
    return time_in / time_in.sum()

print("race      :", np.round(occupancy(True), 3))
print("reassemble:", np.round(occupancy(False), 3))
# both print the same occupancy vector -> the two recipes are one process

The two printed vectors agree to two or three decimals, which is the equivalence made empirical: whether you race the clocks jointly or reassemble the dwell and destination separately, the chain spends the same fraction of time in each state. The factorization is not a modeling choice — it is a theorem you can measure.

Why the jump chain has no self-loops

A subtle but important point: the embedded jump chain always has Pii = 0. In continuous time, “jumping to where you already are” is meaningless — a jump is by definition a change of state. The holding-time clock already accounts for “staying put for a while”; the jump chain only records transitions. So a CTMC’s jump chain is a discrete Markov chain with a guaranteed zero diagonal.

Contrast this with a discrete-time Markov chain, whose transition matrix routinely has a positive diagonal Pii > 0 meaning “stay this tick.” That self-loop probability is exactly the artifact of the tick — the “nothing happened” outcome. Removing the tick removes the self-loop: in continuous time, time-in-state lives in the clock, and transitions live in the (loop-free) jump chain. The two responsibilities are cleanly separated.

Comparing the two simulation styles

The reassemble style is what real simulators use, because it needs only two random draws per jump regardless of how many exits a state has, whereas the race draws one clock per exit. For a state with 50 possible destinations, reassemble draws 2 numbers; the naive race draws 50. Same distribution, far cheaper — which is exactly the Gillespie algorithm of Chapter 9.

From-scratch: recover Q from a simulated path

The deepest test of “(P, λ) = Q” is the inverse problem: simulate a long trajectory, then estimate the rates back from it. The maximum-likelihood estimate is beautifully simple — q̂ij = (number of i→j jumps) / (total time spent in i). Counts over time, exactly the rate definition.

python
# Simulate a CTMC, then estimate Q back: q_ij = (i->j jumps) / (time in i).
import numpy as np

Q_true = np.array([[-4.0, 1.0, 3.0], [2.0, -2.0, 0.0], [0.0, 1.0, -1.0]])
lam = -np.diag(Q_true)
P = Q_true / lam[:, None]; np.fill_diagonal(P, 0.0)

t, s = 0.0, 0
time_in = np.zeros(3); jump_ct = np.zeros((3, 3))
while t < 500000:
    dwell = np.random.exponential(1 / lam[s]); nxt = np.random.choice(3, p=P[s])
    time_in[s] += dwell; jump_ct[s, nxt] += 1; t += dwell; s = nxt

Q_hat = jump_ct / time_in[:, None]            # off-diagonal MLE rates
for i in range(3): Q_hat[i, i] = -Q_hat[i].sum()    # fill diagonal
print("recovered Q:")
print(np.round(Q_hat, 2))                          # ~ Q_true
# row 0 ~ [-4, 1, 3]: counts/time recover the generating rates

The estimated Q̂ recovers the true generator to about two decimals: counting i→j transitions and dividing by time-in-i gives back the rates. This closes the loop — rates generate paths (forward), and paths estimate rates (inverse) — and it is the exact M-step that Chapter’s sibling CT-HMM lesson uses to learn a generator from data.

A subtle trap: jump-chain π ≠ CTMC π

One genuine gotcha the split exposes: the stationary distribution of the embedded jump chain is not the same as the CTMC’s time-stationary π. The jump chain counts visits (how often you land in a state); the CTMC π weights by time (how long you stay). A state you visit often but leave instantly (high λ) gets a big jump-chain share but a small time-share.

The conversion is exactly the holding time: πi ∝ νi/λi, where ν is the jump-chain stationary distribution. Divide the visit-frequency by the leave-rate (multiply by the mean dwell) and renormalize to get the time-stationary π. Forgetting this division is a classic CTMC error — it conflates “how often visited” with “how much time spent.”

python
# Show jump-chain pi (visits) differs from CTMC pi (time), and convert.
import numpy as np

Q = np.array([[-4.0, 1.0, 3.0], [2.0, -2.0, 0.0], [0.0, 1.0, -1.0]])
lam = -np.diag(Q)
P = Q / lam[:, None]; np.fill_diagonal(P, 0.0)

# jump-chain stationary nu: left eigenvector of P for eigenvalue 1
vals, vecs = np.linalg.eig(P.T)
nu = np.abs(vecs[:, np.argmin(np.abs(vals - 1))].real); nu /= nu.sum()

# CTMC time-stationary pi: weight visits by mean dwell 1/lambda, renormalize
pi = (nu / lam); pi /= pi.sum()
print("jump-chain nu (visits):", np.round(nu, 3))   # how OFTEN visited
print("CTMC pi (time)        :", np.round(pi, 3))   # how much TIME spent
# nu and pi DIFFER: state 0 is visited often but left fast (lambda=4),
# so its TIME share is smaller than its VISIT share.

The two distributions genuinely differ: state A is entered frequently but abandoned in ~0.25 time units, so its time-share πA is smaller than its visit-share νA. The ν/λ reweighting is the bridge — and a reminder that “where” and “how long” really are two separate questions, even at equilibrium.

Why this split is the right mental model

Stepping back: the jump-chain-plus-clock decomposition is not just a simulation trick, it is the cleanest way to think about a CTMC. When you look at a system, ask the two questions separately. “Where does it tend to go from here?” is the jump chain — pure routing, no time. “How long does it linger before going?” is the clock — pure timing, no routing. Many modeling mistakes come from tangling these two; the split keeps them honest.

It also explains a deep feature of continuous time: you can change all the holding rates (speed the chain up or slow it down) without touching where it goes, by scaling the diagonal while keeping the off-diagonal ratios fixed. The jump chain is invariant to a global time-rescaling; only the dwell distribution stretches. That separation of “route” from “pace” is unique to the continuous-time view and is exactly what the generator Q = Λ(P − I) factorization (Chapter 5) makes precise.

So carry this forward: every CTMC is a discrete routing chain wearing an exponential wristwatch. Chapter 5 stacks both into one matrix Q; Chapter 9’s Gillespie simulator alternates the two draws forever. The whole rest of the lesson is variations on “where” times “how long.”

python
# Time-rescaling: scale all rates by c -> same jump chain, dwells shrink by c.
import numpy as np

Q = np.array([[-4.0, 1.0, 3.0], [2.0, -2.0, 0.0], [0.0, 1.0, -1.0]])

def jump_and_dwell(Q):
    lam = -np.diag(Q)
    P = Q / lam[:, None]; np.fill_diagonal(P, 0.0)
    return P, 1 / lam            # jump chain, mean dwells

P, dwell = jump_and_dwell(Q)
P2, dwell2 = jump_and_dwell(3.0 * Q)   # triple every rate

print("jump chain unchanged:", np.allclose(P, P2))   # True -- WHERE is invariant
print("dwells scaled by 1/3:", np.allclose(dwell2, dwell / 3))  # True -- only PACE changed
print("original dwells:", np.round(dwell, 3))    # [0.25 0.5  1.  ]
print("3x-rate dwells :", np.round(dwell2, 3))   # [0.083 0.167 0.333]
# Routing (jump chain P) is identical; the chain just runs 3x faster.

# Consequence: the stationary distribution is ALSO unchanged by rescaling,
# because pi depends only on the balance of flows, not their absolute speed.
from scipy.linalg import expm
pi_a = expm(Q * 10000)[0]
pi_b = expm(3 * Q * 10000)[0]
print("pi unchanged by 3x rate:", np.allclose(pi_a, pi_b))   # True
print("pi (either) =", np.round(pi_a, 3))                # same stationary vector
# Only the TIME to reach pi changes (3x faster); the destination is identical.
# This is why "route" and "pace" are genuinely independent knobs on a CTMC:
# the off-diagonal RATIOS set where it goes, the diagonal MAGNITUDES set how fast.
print("mean dwell A, 1x vs 3x:", round(dwell[0], 3), "vs", round(dwell2[0], 3))  # 0.25 vs 0.083
print("P[0] row, 1x vs 3x:", np.round(P[0], 3), np.round(P2[0], 3))  # identical routing
# Same routing row, faster clock: the canonical "where vs how-long" separation.
# Knowing (P, lambda) is knowing Q -- and you may rescale lambda without
# touching P, the cleanest proof that the two halves are truly independent.
# Gillespie (Chapter 9) exploits exactly this: draw dwell from lambda, jump via P.
# Two random numbers per jump, independent of how many exits the state has.
# That O(1)-draws-per-jump efficiency is why event-driven simulation scales.

Tripling every rate leaves the jump chain P identical while dividing all dwell times by three — the route is unchanged, only the pace. The stationary distribution is unchanged too, since π depends on the balance of flows, not their absolute speed. This cleanly separates “where it goes” from “how fast it gets there,” the deepest consequence of the where/how-long split and a property no discrete-time chain has.

Chapter 5: The Generator Matrix Q

So far the rates have been scattered around: fail 0.2, repair 1.0, maintenance 0.05. To do anything systematic — compute finite-time answers, find equilibria, simulate — we need one tidy ledger that holds every rate. That ledger is the generator matrix Q, and it is the single object that generates the entire process.

The rule for filling it in has two parts. The off-diagonal entry qij ≥ 0 (for i ≠ j) is the rate of jumping from state i to state j — exactly the clock rates from the race. The diagonal entry is set so the row balances: qii = −Σj≠i qij, the negative of the total leave-rate.

That single rule — every row of Q sums to zero — is the defining property of a generator, and it is just conservation of probability written as arithmetic. The diagonal says “probability leaks OUT of state i at total rate λ_i” (negative, because it is leaving); the off-diagonals say “and here is where every bit of it goes.” Out equals in-redirected; nothing is created or destroyed.

It pays to read Q geometrically. A row tells you everything about leaving a state: the diagonal is the total exit rate (so the holding time is exp(−q_ii)), and the off-diagonals are how that exit splits across destinations (so the jump-chain row is q_ij/(−q_ii)). Both Chapter 2’s clock and Chapter 4’s jump chain are sitting inside a single row of Q.

Units matter and are a frequent source of confusion. Every entry of Q has units of 1/time — they are rates, not probabilities. You cannot read qUP,DOWN = 0.2 as “a 0.2 probability of going down.” It is “a tendency of 0.2 per hour.” The dimensionless probability object (rows summing to 1) is P(t), which Q generates — that is the next two chapters.

One more structural fact for later: because every row sums to zero, the all-ones vector is killed by Q on the right, Q·1 = 0. That means Q always has an eigenvalue of exactly 0, with the all-ones right-eigenvector. This is not a curiosity — it is the mode that does not decay, and in Chapter 7–8 it carries the stationary distribution. The conservation rule and the long-run behavior are the same fact.

It also helps to see Q as Q = Λ(P − I), where Λ is the diagonal of holding rates and P is the embedded jump matrix. The off-diagonal of Λ(P−I) is λiPij = qij (rate = holding-rate times jump-probability), and the diagonal is λi(0 − 1) = −λi. This factorization is the cleanest statement of “a CTMC is a jump chain P timed by an exponential clock Λ” — Chapter 4’s split, written as matrices.

Uniformization: turning Q into a probability matrix (preview)

There is a slick way to convert the rate-world Q back into a probability-world matrix, used heavily in practice. Pick a rate μ at least as big as the fastest leave-rate (μ ≥ maxi λi). Then U = I + Q/μ is an honest stochastic matrix: its rows sum to 1, and all entries are nonnegative because μ is big enough to keep 1 + qii/μ ≥ 0.

This uniformization reinterprets the CTMC as a discrete chain U that ticks at a fixed fast rate μ, where most ticks are self-loops (“nothing happened”). It is the principled version of Chapter 0’s tick — and it gives a numerically stable series for exp(Qt) in Chapter 7: exp(Qt) = e−μt Σk (μt)k/k! · Uk, a Poisson-weighted average of powers of U. We only flag it here; the point is that Q and a probability matrix are two encodings of the same process, related by a single shift-and-scale.

Concretely, for our 2-state Q with μ = 1.0 (the max leave-rate): U = I + Q/1.0 = [[0.8, 0.2],[1.0, 0.0]]. Row UP keeps 80% probability on itself and sends 20% to DOWN; row DOWN sends everything back to UP. Rows sum to 1, entries nonnegative — a genuine stochastic matrix carved out of the generator.

Reading a single row three ways

Every row of Q is a complete description of leaving one state, and it answers three of our earlier questions at once. Take row UP = [−0.2, 0.2]. Chapter 2 (when): the diagonal magnitude 0.2 is the total leave-rate λ, so the holding time is exp(0.2) with mean 5 hr. Chapter 3/4 (where): the off-diagonal 0.2 divided by λ gives the jump-chain probability 0.2/0.2 = 1.0 — UP always goes to DOWN. Chapter 6 (flow): the row is the instantaneous rate of probability leaving (diagonal) and arriving (off-diagonals).

So you never have to “remember” the holding distribution and the jump chain as separate objects — they are both encoded in one row, recoverable by reading the diagonal versus the off-diagonals. Q is the compression of everything in Chapters 2–4 into a single grid.

Worked example: a 3-state generator

Add a MAINTENANCE state to the machine. From UP: fail to DOWN at 0.2, go to MAINT at 0.05. From DOWN: repair to UP at 1.0. From MAINT: return to UP at 0.5. Build Q row by row (order UP, DOWN, MAINT).

Row UP. Off-diagonals 0.2 (to DOWN) and 0.05 (to MAINT); their sum is 0.25, so the diagonal is −0.25. Row: [−0.25, 0.2, 0.05], sum = 0 ✓.

Row DOWN. One off-diagonal 1.0 (to UP); diagonal −1.0. Row: [1.0, −1.0, 0], sum = 0 ✓.

Row MAINT. One off-diagonal 0.5 (to UP); diagonal −0.5. Row: [0.5, 0, −0.5], sum = 0 ✓.

Read row UP back: holding time exp(0.25) (mean 4 hr, matching Chapter 3’s race), and the jump chain from UP is [0, 0.2/0.25, 0.05/0.25] = [0, 0.80, 0.20] — 80% to DOWN, 20% to MAINT, exactly the race result. The whole multi-exit story of Chapter 3 is sitting in one row of this matrix.

Worked example: the canonical UP/DOWN generator

Two states, UP and DOWN. The machine fails (UP→DOWN) at rate a = 0.2/hr and is repaired (DOWN→UP) at rate b = 1.0/hr. Build Q row by row.

Row UP. The only off-diagonal is qUP,DOWN = 0.2. So the diagonal is qUP,UP = −0.2. Row check: −0.2 + 0.2 = 0 ✓.

Row DOWN. The only off-diagonal is qDOWN,UP = 1.0. So the diagonal is qDOWN,DOWN = −1.0. Row check: 1.0 − 1.0 = 0 ✓.

Read it back: the diagonal −0.2 means “UP leaks out at total rate 0.2” (holding time exp(0.2), mean 5 hr); the off-diagonal 0.2 says “all of that leak goes to DOWN.” The diagonal −1.0 means DOWN leaks at rate 1.0 (mean stay 1 hr); the 1.0 says it all returns to UP. This single Q drives every computation in the rest of the lesson.

The quiz, worked

A valid generator row reads [−0.7, qAB, 0.3]. The row must sum to zero: −0.7 + qAB + 0.3 = 0, so qAB = 0.7 − 0.3 = 0.4. And it must be nonnegative (it is an outgoing rate), which 0.4 satisfies. So qAB = 0.4.

python
# Assemble Q from a dict of (i,j)->rate, auto-fill the diagonal, then validate.
import numpy as np

def build_Q(n, rate_dict):
    Q = np.zeros((n, n))
    for (i, j), r in rate_dict.items():
        Q[i, j] = r                  # off-diagonal outgoing rate
    for i in range(n):
        Q[i, i] = -(Q[i].sum() - Q[i, i])  # diagonal = -(row off-diag sum)
    return Q

def validate(Q, tol=1e-9):
    rows_zero = np.allclose(Q.sum(axis=1), 0.0, atol=tol)
    offdiag_ok = (Q - np.diag(np.diag(Q)) >= -tol).all()
    return rows_zero and offdiag_ok

Q = build_Q(2, {(0, 1): 0.2, (1, 0): 1.0})
print(Q)            # [[-0.2  0.2] [ 1.  -1. ]]
print(validate(Q))  # True: rows sum to 0, off-diagonals >= 0

The library one-liner for the row-sum check is just np.allclose(Q.sum(axis=1), 0) — that one assertion is the whole definition of “is this a valid generator.”

The zero eigenvalue, from the row-sum rule

One structural consequence deserves its own derivation because it powers Chapters 7 and 8. Because every row sums to zero, multiplying Q by the all-ones column vector 1 = [1,1,…,1]T gives each entry as that row’s sum — which is zero. So Q·1 = 0 = 0·1.

That is literally the eigenvalue equation Q·v = λ·v with eigenvector v = 1 and eigenvalue λ = 0. Hence a generator always has 0 as an eigenvalue, with the all-ones right-eigenvector. In Chapter 7 this is the mode of exp(Qt) that does not decay (since e0·t = 1); in Chapter 8 its left-eigenvector is the stationary distribution π. Conservation of probability and the existence of equilibrium are the same algebraic fact.

From-scratch: validate, convert, and read a generator

python
# A complete generator toolkit: validate, extract holding rates and jump chain,
# and confirm the all-ones vector is the zero-eigenvalue eigenvector.
import numpy as np

Q = np.array([[-0.25, 0.2, 0.05],
              [ 1.0, -1.0, 0.0],
              [ 0.5, 0.0, -0.5]])

assert np.allclose(Q.sum(axis=1), 0)              # rows sum to 0
offdiag = Q - np.diag(np.diag(Q))
assert (offdiag >= -1e-12).all()                    # off-diagonals nonnegative

lam = -np.diag(Q)                                   # holding rates per state
Pjump = offdiag / lam[:, None]                       # embedded jump chain
print("holding rates :", lam)             # [0.25 1.  0.5]
print("mean dwells   :", 1 / lam)         # [4.  1.  2. ]
print("jump row UP   :", Pjump[0])         # [0.   0.8  0.2]  -> Ch3 race

ones = np.ones(3)
print("Q @ 1 =", Q @ ones)                      # [0. 0. 0.] -> eigenvalue 0
vals = np.linalg.eigvals(Q)
print("eigenvalues:", np.round(vals, 3))      # one is exactly 0, others < 0

The printout ties the chapter together: the holding rates are the diagonal magnitudes (Chapter 2), the jump row is the normalized off-diagonals (Chapters 3–4), and Q@1 = 0 confirms the zero eigenvalue (Chapters 7–8). Four lessons, one matrix, one validation pass.

Worked example: uniformization by hand

Convert the 2-state Q = [[−0.2, 0.2],[1.0, −1.0]] into a stochastic matrix. The max leave-rate is μ = 1.0 (state DOWN). Form U = I + Q/μ.

Row UP: [1, 0] + [−0.2, 0.2]/1.0 = [0.8, 0.2]. Row DOWN: [0, 1] + [1.0, −1.0]/1.0 = [1.0, 0.0]. So U = [[0.8, 0.2],[1.0, 0.0]]. Check: both rows sum to 1, all entries in [0,1] — a genuine transition matrix. The self-loop probability 0.8 on row UP is the “nothing happened this μ-tick” probability; the 0.2 is a real UP→DOWN move.

This is the rigorous version of Chapter 0’s tick: instead of guessing a dt, you let the Poisson clock at rate μ decide the ticks, and most ticks are self-loops. It is also the bridge to Chapter 7’s exact exp(Qt), which is a Poisson-weighted sum of powers of U.

python
# Uniformization: Q -> stochastic U, then recover exp(Qt) as a Poisson sum.
import numpy as np, math
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
mu = -np.diag(Q).max()              # at least the fastest leave-rate
U = np.eye(2) + Q / mu               # the uniformized stochastic matrix
print("U =", U)                          # [[0.8 0.2] [1. 0. ]] rows sum to 1
assert np.allclose(U.sum(axis=1), 1) and (U >= 0).all()

def expm_uniform(Q, t, mu, K=60):
    U = np.eye(len(Q)) + Q / mu
    P = np.zeros_like(U); Uk = np.eye(len(Q))
    for k in range(K):
        wk = math.exp(-mu * t) * (mu * t) ** k / math.factorial(k)  # Poisson weight
        P += wk * Uk; Uk = Uk @ U
    return P
print(np.round(expm_uniform(Q, 3, mu)[0], 4))   # [0.8379 0.1621]
print(np.round(expm(Q * 3)[0], 4))           # [0.8379 0.1621] -> agree

The Poisson-weighted sum of powers of U reproduces the exact exp(Qt) answer [0.8379, 0.1621] — numerically stable because U has nonnegative entries and the weights are an honest probability distribution. This is the third independent route to P(t), after Euler (Ch6) and eigendecomposition (Ch7), all landing on the same verified numbers.

Worked example: round-tripping Q ↔ (P, λ)

Chapter 4 claimed that (P, λ) carries exactly the same information as Q. Verify the round trip on the 3-state generator. Forward (Q → P, λ): holding rates are the diagonal magnitudes λ = [0.25, 1.0, 0.5]; jump rows are off-diagonals over λ, giving P = [[0, 0.8, 0.2],[1, 0, 0],[1, 0, 0]].

Backward (P, λ → Q): rebuild each off-diagonal as qij = λi·Pij. Row UP: 0.25×0.8 = 0.2 (to DOWN), 0.25×0.2 = 0.05 (to MAINT); diagonal −0.25. That is exactly the original row [−0.25, 0.2, 0.05]. No information was lost in the split — it is a reversible change of coordinates, as promised.

python
# Round-trip: Q -> (P, lambda) -> Q must return the original generator.
import numpy as np

Q = np.array([[-0.25, 0.2, 0.05],
              [ 1.0, -1.0, 0.0],
              [ 0.5, 0.0, -0.5]])

def to_jump(Q):
    lam = -np.diag(Q)
    P = (Q - np.diag(np.diag(Q))) / lam[:, None]   # off-diag / leave-rate
    return P, lam

def from_jump(P, lam):
    Q = P * lam[:, None]               # q_ij = lambda_i * P_ij
    for i in range(len(lam)): Q[i, i] = -lam[i]   # restore the diagonal
    return Q

P, lam = to_jump(Q)
print("jump P :", np.round(P, 3).tolist())
print("lambda :", lam)                       # [0.25 1.   0.5 ]
print("round-trip OK:", np.allclose(from_jump(P, lam), Q))  # True

The reconstructed Q is identical to the original to machine precision — the “where” (jump chain P) and the “when” (holding rates λ) together hold every bit of information in the generator. This is Chapter 4’s factoring, made invertible in code.

Worked example: validating a malformed generator

Practice spotting invalid generators, because mistakes here are common. Candidate A: [[−0.3, 0.2, 0.1],[0.5, −0.5, 0],[0, 1, −1]] — row sums are 0, 0, 0 ✓ and off-diagonals nonnegative ✓, so it is valid. Candidate B: [[−0.3, 0.2, 0.2], …] — row 0 sums to 0.1 ≠ 0, so it is not a generator (probability is being created). Candidate C: [[0.3, −0.2, −0.1], …] — the diagonal is positive and off-diagonals negative, the signs inverted; invalid (rates cannot be negative, diagonals cannot be positive).

The two checks — rows sum to zero, off-diagonals ≥ 0 — catch every malformed generator. They are the rate-world analog of “rows sum to 1, entries in [0,1]” for a probability matrix, and forgetting them is the single most common CTMC bug.

python
# A strict generator validator with informative failure messages.
import numpy as np

def check_generator(Q):
    Q = np.asarray(Q, float)
    problems = []
    if not np.allclose(Q.sum(axis=1), 0, atol=1e-9):
        problems.append(f"rows must sum to 0, got {Q.sum(axis=1)}")
    off = Q - np.diag(np.diag(Q))
    if (off < -1e-12).any():
        problems.append("off-diagonals must be >= 0 (they are rates)")
    if (np.diag(Q) > 1e-12).any():
        problems.append("diagonals must be <= 0 (negative leave-rate)")
    return problems or ["valid generator"]

print(check_generator([[-0.3,0.2,0.1],[0.5,-0.5,0],[0,1,-1]]))  # valid
print(check_generator([[-0.3,0.2,0.2],[0.5,-0.5,0],[0,1,-1]]))  # row 0 sums to 0.1
print(check_generator([[ 0.3,-0.2,-0.1],[0.5,-0.5,0],[0,1,-1]]))  # sign errors

The validator names exactly what is wrong: a non-zero row sum (probability not conserved) or sign errors (negative rates / positive diagonal). Running this assertion before any computation catches the bugs that otherwise surface as nonsensical probabilities downstream.

Worked example: from a labeled diagram to Q

In practice you start from a hand-drawn state diagram with labeled arrows, not a matrix. Suppose: HEALTHY→INFECTED at rate 0.3, INFECTED→RECOVERED at rate 0.5, INFECTED→DEAD at rate 0.1, RECOVERED→HEALTHY at rate 0.2 (immunity wanes), DEAD absorbing. Translate the arrows directly into off-diagonal entries, then fill diagonals to balance each row.

Row HEALTHY: one arrow out (0.3 to INFECTED), diagonal −0.3. Row INFECTED: two arrows out (0.5 + 0.1 = 0.6), diagonal −0.6. Row RECOVERED: one arrow (0.2 to HEALTHY), diagonal −0.2. Row DEAD: no arrows, diagonal 0 (absorbing). Every arrow becomes one positive off-diagonal; every row’s diagonal is set to make it sum to zero. That mechanical translation — arrows to off-diagonals, balance the diagonal — is the entire skill of building a generator from a model.

python
# Build Q from a labeled edge list (a state diagram), auto-balancing diagonals.
import numpy as np

states = ["HEALTHY", "INFECTED", "RECOVERED", "DEAD"]
idx = {s: i for i, s in enumerate(states)}
edges = [("HEALTHY", "INFECTED", 0.3),
         ("INFECTED", "RECOVERED", 0.5),
         ("INFECTED", "DEAD", 0.1),
         ("RECOVERED", "HEALTHY", 0.2)]

Q = np.zeros((4, 4))
for src, dst, rate in edges:
    Q[idx[src], idx[dst]] = rate        # each arrow -> one off-diagonal
for i in range(4):
    Q[i, i] = -Q[i].sum()                # balance each row to zero
print(Q)
print("rows sum to 0:", np.allclose(Q.sum(1), 0))   # True
print("DEAD absorbing:", np.allclose(Q[3], 0))     # True (no out-arrows)

The edge list translates one-to-one into Q: every labeled arrow is a positive off-diagonal, and the diagonal auto-balances to keep rows summing to zero. DEAD has no out-arrows, so its row is all zeros — an absorbing trap (Chapter 10). This is the workflow for every real model: draw the diagram, list the arrows, let the diagonals fall out.

Q is everything the rest of the lesson needs

Pause to appreciate what we now hold. The single matrix Q determines the holding time of every state (its diagonal), the destination probabilities from every state (its off-diagonals), the finite-time transition probabilities (via exp(Qt), Chapter 7), the long-run distribution (via πQ=0, Chapter 8), and the recipe for simulation (Chapter 9). Nothing else is needed; Q is the complete specification of a CTMC.

This is the payoff of the “rates not probabilities” choice. A discrete chain needs a different transition matrix for every time horizon you care about; a CTMC needs only Q, and every horizon’s answer is exp(Qt) for the appropriate t. One matrix, parameterized by time, instead of a new matrix per tick. That compression is exactly what we bought by removing the clock in Chapter 0.

python
# One Q, every horizon: P(t) for any t is just exp(Q*t). No new matrices.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
for t in [0.5, 1, 3, 10]:
    print(f"P(UP up after {t:4}h) = {expm(Q*t)[0,0]:.4f}")
# 0.5h: 0.9067   1h: 0.8612   3h: 0.8379   10h: 0.8334 (-> pi)
# Every horizon from the SAME Q -- the generator is the whole model.
print("limit (t->inf):", np.round(expm(Q * 1000)[0, 0], 4))   # 0.8333 = pi_UP
# short horizons -> near 1 (just started UP); long horizons -> pi (forgot start)
# The discrete-chain alternative would need a SEPARATE matrix for each of these
# horizons; the CTMC reuses the one Q via exp(Q*t). That is the compression.
print("P(2.5h):", np.round(expm(Q * 2.5)[0], 4))   # any t you want, one Q
# Q is the COMPLETE model: holding times, jump chain, P(t), pi -- all from here.
# Read the diagonal for dwell, off-diagonals for routing, exp(Qt) for any horizon.
# Validate it with one line: assert np.allclose(Q.sum(axis=1), 0).

The single generator Q produces the answer for every time horizon through exp(Qt) — 0.5 hr, 3 hr, or the long-run limit, all from one matrix. That is the compression the generator buys, and the rest of the lesson is learning to compute and read it.

Chapter 6: Kolmogorov’s Forward & Backward Equations

We have the instantaneous rates (Q). Now we want the finite-time answer: “starting UP, what is the chance I am DOWN after 3 hours?” That quantity is the transition function Pij(t) — the probability of being in state j at time t given you started in state i. The job of this chapter is to connect the rates Q to the finite-time P(t).

The connection comes from asking how P changes over one tiny sliver dt. Over dt, the chain either stays put or makes one jump, and the jump probabilities are exactly Chapter 1’s linear rule: P(t+dt) ≈ P(t)·(I + Q·dt). The I is “stay,” the Q·dt is the rate-weighted “move.”

Rearrange and let dt→0. (P(t+dt) − P(t))/dt ≈ P(t)·Q, which in the limit is a differential equation: dP/dt = P·Q, with P(0) = I. This is the Kolmogorov forward equation. It is the matrix version of Chapter 1’s scalar seed dx/dt = q·x — same shape, more dimensions.

There is a twin. If you instead condition the tiny step on the first instant rather than the last, you get dP/dt = Q·P, the Kolmogorov backward equation. Both equations have the same solution P(t); they differ only in which end you hold fixed while differentiating. Forward fixes the start and varies the end-time (good for propagating an initial distribution forward); backward fixes the end and varies the start (good for hitting/absorption questions, Chapter 10).

The picture to hold is probability as a fluid sloshing between states. The current distribution is a set of bar heights summing to 1. Each small step, mass flows along every edge in amount (source probability)·(rate)·dt. Inflow to one state is outflow from another, so the total height stays exactly 1 — conservation, every step.

We can solve this numerically right now with the crudest possible method — Euler integration: take the distribution, add P·Q·dt, repeat. It is exactly the “tiny tick” we mocked in Chapter 0, and it works as an approximation, but it has error that grows with dt. Chapter 7 will replace it with the exact closed-form solution, no stepping at all. Euler here is the bridge that makes the exact form unsurprising.

Deriving the forward equation from Chapman–Kolmogorov

The honest derivation starts from a self-consistency rule called Chapman–Kolmogorov: getting from i to j in time t + dt means getting to some intermediate state k in time t, then from k to j in the last dt. In matrix form that is just P(t + dt) = P(t)·P(dt) — transition over a long interval factors through any split point.

Now use Chapter 1’s linear rule for the short factor: over a tiny dt, P(dt) ≈ I + Q·dt (stay, plus rate-weighted move). Substitute: P(t+dt) ≈ P(t)·(I + Q·dt) = P(t) + P(t)·Q·dt.

Subtract P(t) and divide by dt: (P(t+dt) − P(t))/dt ≈ P(t)·Q. Take dt→0 and the left side becomes the derivative by definition: dP/dt = P·Q, with P(0) = I. That is the Kolmogorov forward equation, derived in three lines from “split the interval” plus “short steps move at rate Q.”

The backward equation comes from splitting at the other end: P(t+dt) = P(dt)·P(t) = (I + Q·dt)·P(t), giving dP/dt = Q·P. Same short factor, placed on the left instead of the right. Both are exact in the limit, and because matrix multiplication here commutes with the shared solution, both yield P(t) = exp(Qt).

Worked example: one Euler step from fully UP

Use Q = [[−0.2, 0.2], [1.0, −1.0]], start fully UP so the row vector is P = [1, 0], and take one step of dt = 0.1.

Compute P·Q. Row-vector times matrix: [1, 0]·Q. First entry: 1·(−0.2) + 0·(1.0) = −0.2. Second entry: 1·(0.2) + 0·(−1.0) = 0.2. So P·Q = [−0.2, 0.2].

Scale by dt. P·Q·dt = [−0.2, 0.2]×0.1 = [−0.02, 0.02]. This is the change vector dP.

Add to P. Pnew = [1, 0] + [−0.02, 0.02] = [0.98, 0.02].

Read it: after 0.1 hr, there is a 98% chance still UP and a 2% chance now DOWN. And the conservation check: 0.98 + 0.02 = 1.000 exactly — mass was moved, not created. (For comparison, the exact answer from Chapter 7 is 0.9803 / 0.0197; Euler’s 0.98 / 0.02 is close, and would be closer with a smaller dt.)

Worked example: three Euler steps, watching error accumulate

Keep stepping from P = [1, 0] with dt = 0.1 to see the error compound. After step 1 (above): [0.98, 0.02]. Step 2: compute P·Q = [0.98, 0.02]·Q. First entry 0.98·(−0.2) + 0.02·(1.0) = −0.196 + 0.02 = −0.176; second 0.98·(0.2) + 0.02·(−1.0) = 0.196 − 0.02 = 0.176. Scale by 0.1: [−0.0176, 0.0176]. Add: [0.9624, 0.0376].

Step 3: P·Q = [0.9624, 0.0376]·Q. First 0.9624·(−0.2) + 0.0376·(1.0) = −0.19248 + 0.0376 = −0.15488; scale and add gives [0.94691, 0.05309]. After 0.3 hr Euler says ~0.9469 UP. The exact exp(0.3Q) gives ~0.9474 — the gap is ~0.0005 and grows roughly linearly in the number of steps for a fixed dt. Halve dt (double the steps) and the per-step error quarters but you take twice as many, so total error halves: Euler is first-order accurate.

This is exactly why Chapter 7’s closed form matters: exp(Qt) has zero discretization error at any t, for one matrix operation, instead of a tower of accumulating Euler steps.

The quiz, worked

Start from P(0) = I and take one tiny step of the forward equation. P(dt) ≈ P(0) + P(0)·Q·dt = I + Q·dt. So after a tiny dt, P ≈ I + Q·dt — probabilities move at exactly the rates in Q. This is the first-order term of the matrix exponential we build next.

python
# Euler-integrate dP/dt = P@Q from P = I; row-sums stay ~1 (conservation).
import numpy as np

def euler_Pt(Q, t, steps):
    Q = np.asarray(Q, float)
    dt = t / steps
    P = np.eye(Q.shape[0])      # P(0) = I
    for _ in range(steps):
        P = P + P @ Q * dt       # forward step: P += P@Q*dt
    return P

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
print(euler_Pt(Q, 0.1, 1))     # [[0.98 0.02] [0.1 0.9]] one coarse step
print(euler_Pt(Q, 3.0, 3000)[0])  # ~[0.838 0.162] -> matches exp(3Q) row 0

The library one-liner that makes Euler obsolete: scipy.linalg.expm(Q*t) — the exact P(t), no stepping. That is Chapter 7.

From-scratch: Euler vs the exact answer as dt shrinks

python
# Show Euler is first-order: total error ~ halves when dt halves.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
t = 3.0
exact = expm(Q * t)                       # the truth, no stepping

def euler_Pt(Q, t, steps):
    dt = t / steps; P = np.eye(len(Q))
    for _ in range(steps): P = P + P @ Q * dt
    return P

for steps in [3, 30, 300, 3000]:
    approx = euler_Pt(Q, t, steps)
    err = np.abs(approx - exact).max()
    print(f"steps={steps:5d} (dt={t/steps:.4f}): max err={err:.6f}")
# steps=3   (dt=1.0000): max err~0.05      (coarse: visibly off)
# steps=30  (dt=0.1000): max err~0.005     (10x finer -> ~10x smaller)
# steps=300 (dt=0.0100): max err~0.0005    (first-order: error ~ dt)
# steps=3000(dt=0.0010): max err~0.00005

# Conservation holds at every step: rows of the iterate sum to ~1.
print(euler_Pt(Q, t, 300).sum(axis=1))   # [1. 1.] -> mass conserved

Each tenfold refinement shrinks the error roughly tenfold — the signature of a first-order method, where total error scales like dt. Chapter 7 replaces this whole ladder with one exact call; here the value is seeing why the exact form is worth wanting.

Two readings: distribution row-vector vs full matrix

The flow equation comes in two practical sizes. If you only care about one starting distribution, propagate a row vector: dp/dt = p·Q, p(0) = your start. That is cheap — one vector evolving. If you want the answer from every start at once, propagate the full matrix dP/dt = P·Q, P(0) = I, whose row i is the distribution started from state i.

They are the same equation; the matrix version just runs n row-vector versions in parallel (one per identity row). When you finally reach exp(Qt) in Chapter 7, the matrix gives you the complete transition table and the row-vector gives you p(0)·exp(Qt), the evolved distribution. Choose the vector when you have a known start, the matrix when you want the whole story.

From-scratch: a better integrator (RK4) closes the gap faster

Euler is first-order; a fourth-order Runge–Kutta (RK4) step uses four rate evaluations to cancel more error terms, converging like dt⁴ instead of dt. It is still “stepping,” but dramatically more accurate per step — useful when you cannot use expm (e.g. time-varying Q).

python
# Euler vs RK4 for dP/dt = P@Q: RK4 nails it in far fewer steps.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
t, exact = 3.0, expm(Q * 3.0)

def step_euler(P, dt): return P + P @ Q * dt
def step_rk4(P, dt):
    k1 = P @ Q
    k2 = (P + dt / 2 * k1) @ Q
    k3 = (P + dt / 2 * k2) @ Q
    k4 = (P + dt * k3) @ Q
    return P + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

for steps in [6, 30]:
    dt = t / steps
    Pe = Pr = np.eye(2)
    for _ in range(steps): Pe = step_euler(Pe, dt); Pr = step_rk4(Pr, dt)
    print(f"steps={steps}: euler err={np.abs(Pe-exact).max():.2e}  rk4 err={np.abs(Pr-exact).max():.2e}")
# steps=6:  euler err~3e-02  rk4 err~2e-05   (RK4 ~1000x better at same cost class)
# steps=30: euler err~5e-03  rk4 err~5e-08

At just 6 steps RK4 is already a thousand times more accurate than Euler, because it cancels the leading error terms. Both are still beaten by expm (zero discretization error), but RK4 is the tool of choice when Q itself varies with time and no closed form exists.

Worked example: the probability-flow bookkeeping of one step

Let us track the actual mass moving in one Euler step, edge by edge, to see conservation as plumbing. Start P = [0.6, 0.4] (60% UP, 40% DOWN), dt = 0.1. The UP→DOWN flow is (source prob)·(rate)·dt = 0.6·0.2·0.1 = 0.012. The DOWN→UP flow is 0.4·1.0·0.1 = 0.040.

Update UP: it loses 0.012 (out to DOWN) and gains 0.040 (in from DOWN), net +0.028 → 0.6 + 0.028 = 0.628. Update DOWN: gains 0.012, loses 0.040, net −0.028 → 0.4 − 0.028 = 0.372. Total: 0.628 + 0.372 = 1.000 — every unit that left one bar arrived in the other. That is why row-sums stay exactly 1: the off-diagonal flows are equal-and-opposite ledger entries.

This edge-by-edge view is precisely what the widget animates: the arrow widths are these prob·rate·dt flows, and the bar heights are the running totals. Conservation is not imposed — it falls out of the row-sum-zero structure of Q (Chapter 5), because what Q subtracts on the diagonal it adds back across the off-diagonals.

python
# One Euler step as explicit edge flows; confirm mass is conserved.
import numpy as np

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
P = np.array([0.6, 0.4]); dt = 0.1

flow_ud = P[0] * Q[0, 1] * dt      # UP -> DOWN mass
flow_du = P[1] * Q[1, 0] * dt      # DOWN -> UP mass
print(f"UP->DOWN flow = {flow_ud:.3f}   DOWN->UP flow = {flow_du:.3f}")

P_new = np.array([P[0] - flow_ud + flow_du,    # UP: -out +in
                  P[1] - flow_du + flow_ud])   # DOWN: -out +in
print("P_new =", np.round(P_new, 3), " sum =", P_new.sum())   # [0.628 0.372] sum 1.0
print("matches P + P@Q*dt:", np.allclose(P_new, P + P @ Q * dt))  # True

The explicit flow bookkeeping reproduces P + P@Q*dt exactly, and the sum stays 1.000. The matrix form P@Q is just these edge flows added up per state — the diagonal collects the outflows, the off-diagonals deliver the inflows.

Worked example: propagating a known starting distribution

Often you do not want the whole matrix — just “given I start 70% UP / 30% DOWN, where am I after a short time?” Propagate the row vector p = [0.7, 0.3] through dp/dt = p·Q with dt = 0.1. Compute p·Q: first entry 0.7·(−0.2) + 0.3·(1.0) = −0.14 + 0.3 = 0.16; second 0.7·(0.2) + 0.3·(−1.0) = 0.14 − 0.3 = −0.16.

Scale and add: pnew = [0.7, 0.3] + [0.16, −0.16]×0.1 = [0.716, 0.284]. The UP share rose because the strong repair rate is pulling the 30% DOWN mass back faster than the weak fail rate is leaking UP mass away. Net flow is toward UP whenever the current down-share exceeds its equilibrium — the distribution is sliding toward π = [0.833, 0.167].

python
# Propagate a distribution forward two ways: row-vector Euler vs p0 @ expm.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
p0 = np.array([0.7, 0.3])

def propagate_euler(p, T, steps):
    dt = T / steps
    for _ in range(steps): p = p + p @ Q * dt   # one row-vector evolves
    return p

print("one coarse step :", np.round(propagate_euler(p0.copy(), 0.1, 1), 4))   # [0.716 0.284]
print("t=3 (fine Euler):", np.round(propagate_euler(p0.copy(), 3.0, 3000), 4))
print("t=3 (exact expm):", np.round(p0 @ expm(Q * 3), 4))   # sliding toward [0.833 0.167]
# one coarse step : [0.716 0.284]   (matches the hand calc)
# t=3 fine vs exact agree -> Euler converges to p0 @ expm(Q t)

The single coarse step reproduces the hand calculation [0.716, 0.284], and the fine Euler run agrees with the exact p0 @ expm(Q*3) — both slide the distribution toward π. Propagating one vector is far cheaper than the full matrix when you only care about a single known start.

The backward equation in action: a hitting probability

The forward equation propagates a distribution forward; the backward equation answers “starting from each state, what is the probability of hitting a target before time t?” Both share the solution exp(Qt), but the backward view fixes the end condition (the target) and varies the start — ideal for “will I reach state j in time?” questions.

Concretely: for our UP/DOWN chain, “starting UP, probability of being DOWN at t=3” is read off column DOWN of exp(3Q) at row UP — the same matrix, now read by which end state you ask about. Forward and backward are bookend readings of one exponential; the code below shows they agree.

python
# Forward (fix start, vary end) vs backward (fix end, vary start): same expm.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
P3 = expm(Q * 3)            # solves BOTH dP/dt=PQ and dP/dt=QP

# forward: start UP=[1,0], read distribution over end states
print("forward, start UP  :", np.round(np.array([1., 0.]) @ P3, 4))  # [0.8379 0.1621]

# backward: fix end = DOWN (column 1), read prob from each start
print("backward, end DOWN :", np.round(P3[:, 1], 4))   # [0.1621 0.1896]
# row UP of forward == entry [UP,DOWN] of P3 == [UP] of backward column DOWN
print("both agree on P(UP->DOWN,3):", round(P3[0, 1], 4))   # 0.1621

# verify backward equation numerically: d/dt expm = Q @ expm = expm @ Q
h = 1e-6
deriv = (expm(Q * (3 + h)) - expm(Q * 3)) / h
print("PQ == QP == dP/dt:", np.allclose(deriv, P3 @ Q, atol=1e-4) and np.allclose(deriv, Q @ P3, atol=1e-4))  # True

The forward read (start UP) and the backward read (end DOWN) agree on P(UP→DOWN, 3) = 0.1621, and the numerical derivative confirms dP/dt = PQ = QP — the two equations really do share one solution. Pick forward to push a known start ahead, backward to ask about a fixed target from every start.

Why this ODE is the heart of the lesson

The equation dP/dt = PQ is the hinge of the entire subject, so it is worth stating what it buys us. It converts the instantaneous description (rates, the matrix Q) into the finite-time description (probabilities, the matrix P(t)). Everything we want — “chance of being DOWN in 3 hours,” “long-run availability” — is a question about P(t), and this ODE is the bridge from the Q we can write down to the P(t) we actually care about.

Euler integration solves it approximately by taking the “tiny tick” literally, and its first-order error is the price of stepping. Chapter 7 will solve the same ODE exactly, recognizing that a linear matrix ODE has a matrix-exponential solution — the same way the scalar dx/dt = qx has solution eqt. The shape of the equation guarantees the shape of the answer; we are one chapter from the closed form.

So hold the picture: Q is the rate of change operator, P(t) is the accumulated answer, and dP/dt = PQ is the law connecting them. Euler is the crude numerical bridge; exp(Qt) is the exact one. The whole arc from Chapter 5 to Chapter 8 is “write Q, solve this ODE, read off the answer.”

python
# Euler with too-large dt can produce NEGATIVE probabilities -- a stability bug.
import numpy as np

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
for dt in [0.5, 1.5, 2.5]:
    P = np.array([1.0, 0.0])
    for _ in range(10): P = P + P @ Q * dt
    print(f"dt={dt}: P={np.round(P,3)}  valid={bool((P>=0).all())}")
# dt=0.5: P=[0.833 0.167] valid=True   (stable)
# dt=2.5: P=[ ... negative ... ] valid=False  (Euler overshoots past 0!)
# The fix: keep dt < 1/max(leave-rate), or just use expm (always valid).

# expm has NO such constraint -- valid at any t, no step size to tune.
from scipy.linalg import expm
print("expm(Q*25) row 0:", np.round(expm(Q * 25)[0], 4))   # [0.8333 0.1667] valid
# expm never goes negative and never needs a stability check -- it IS the
# exact limit of infinitely-many infinitesimal Euler steps. Use it.
P_unstable = np.array([1.0, 0.0])
for _ in range(10): P_unstable = P_unstable + P_unstable @ Q * 2.5
print("Euler dt=2.5 gives invalid P:", P_unstable)   # negative entries -> garbage
print("any negative? ", bool((P_unstable < 0).any()))   # True -> not a distribution
# expm has no such failure mode -- prefer it whenever a closed form exists.
# Reserve Euler/RK4 for time-VARYING Q where no matrix exponential applies.

With too large a dt, Euler can drive a probability negative — the explicit step overshoots when dt exceeds the fastest state’s time-constant 1/λ. This stability constraint is another reason to prefer the exact exp(Qt), which is a valid stochastic matrix for every t with no step-size to tune. Euler is a teaching bridge, not a production tool.

Chapter 7: P(t) = exp(Qt) — The Matrix Exponential

Chapter 6 left us with dP/dt = PQ, P(0) = I, solved approximately by Euler stepping. Now we solve it exactly. The answer is the most important formula in the subject, and it falls straight out of the scalar analogy we planted in Chapter 1.

Recall the scalar growth law: dx/dt = q·x with x(0) = x0 solves to x(t) = eqt·x0. The matrix equation has the identical shape, so the solution has the identical form: P(t) = exp(Qt), the matrix exponential, where P(0) = I plays the role of x0.

What does “e raised to a matrix” even mean? By definition it is the same power series as the scalar one, with matrix powers: exp(Qt) = I + Qt + (Qt)²/2! + (Qt)³/3! + …. That series always converges, and it is exactly the limit of Chapter 6’s Euler stepping as dt→0 (each factor (I + Q·dt) multiplied infinitely often becomes exp(Qt)).

For a small chain you do not sum the series by hand — you use eigendecomposition. Write Q = V·diag(λ)·V−1 where λ are eigenvalues and V their eigenvectors. Then the powers telescope and exp(Qt) = V·diag(eλt)·V−1. The matrix exponential becomes ordinary scalar exponentials on the eigenvalues — one easy per-mode calculation.

Two facts about the eigenvalues organize everything. First, because rows of Q sum to zero, one eigenvalue is always exactly 0 (the non-decaying mode, carrying the long-run answer). Second, every other eigenvalue of a generator has a negative real part, so its eλt decays. As t grows, the decaying modes die and only the λ=0 mode survives — which is precisely the chain forgetting its start and approaching equilibrium (Chapter 8).

So P(t) is a sum of a constant piece (the λ=0 mode → the stationary distribution) plus exponentially decaying corrections (the other modes). Each entry P_ij(t) starts at the identity value at t=0 and bends smoothly toward its equilibrium value as the corrections fade. Four eigen-numbers describe all the curves.

Why diagonalizing works: powers telescope

The reason eigendecomposition tames the matrix exponential is that it diagonalizes every power at once. If Q = V·D·V−1 with D diagonal, then Q² = V·D·V−1·V·D·V−1 = V·D²·V−1 — the inner V−1V cancels. By induction Qk = V·Dk·V−1.

Substitute into the power series: exp(Qt) = Σk (Qt)k/k! = V·(Σk (Dt)k/k!)·V−1 = V·exp(Dt)·V−1. And because D is diagonal, exp(Dt) is just the scalar exponentials eλt on the diagonal. The whole infinite matrix series collapses to “exponentiate each eigenvalue.”

Deriving the 2-state closed form

For Q = [[−a, a], [b, −b]], find the eigen-pieces. Eigenvalues: 0 (rows sum to zero) and −(a+b) (from the trace). The right eigenvector for 0 is the all-ones [1, 1]T; the left eigenvector for 0 is the stationary [b, a]/(a+b). Assembling V·diag(1, e−(a+b)t)·V−1 and reading the (UP,UP) entry gives the formula we use.

Read the two pieces: the constant b/(a+b) is the λ=0 mode (the equilibrium, Chapter 8); the term a/(a+b)·e−(a+b)t is the λ=−(a+b) mode (the transient that decays). At t = 0 they sum to 1 (you start UP); as t→∞ the transient dies and only the equilibrium remains. Every other entry of P(t) is the same constant-plus-decay shape with different coefficients.

Worked example: UP/DOWN with a = 0.2, b = 1.0 at t = 3

Generator Q = [[−0.2, 0.2], [1.0, −1.0]]. Find its eigenvalues. One is 0 (rows sum to zero). The trace of Q is the sum of eigenvalues: tr(Q) = −0.2 + (−1.0) = −1.2, and since the eigenvalues are 0 and λ2, we get λ2 = −1.2 = −(a+b).

The closed form for the 2-state chain (which you can derive from the eigendecomposition) is PUP→UP(t) = b/(a+b) + a/(a+b)·e−(a+b)t. Plug numbers: b/(a+b) = 1.0/1.2 = 0.8333 and a/(a+b) = 0.2/1.2 = 0.1667. So PUP→UP(t) = 0.8333 + 0.1667·e−1.2t.

Sanity at t = 0. 0.8333 + 0.1667·e0 = 0.8333 + 0.1667 = 1.0000 ✓ (you start UP, so P_UP→UP(0) = 1).

Evaluate at t = 3. e−1.2×3 = e−3.6 = 0.02732. Then PUP→UP(3) = 0.8333 + 0.1667×0.02732 = 0.8333 + 0.00456 = 0.8379.

So P(DOWN at t=3 | start UP) = 1 − 0.8379 = 0.1621. This is the exact answer to the question we have been chasing since Chapter 0, and it matches scipy.linalg.expm(Q*3) to four decimals. Notice the constant term 0.8333 is exactly the equilibrium availability we will derive in Chapter 8.

The quiz, worked

Why is one eigenvalue of Q always exactly 0? Because every row of Q sums to zero, the matrix-vector product Q·1 = 0 for the all-ones vector 1. That is the definition of an eigenvector with eigenvalue 0. This is the mode that does not decay, and it carries the stationary distribution.

python
# P(t) three ways: eigendecomposition, the power series, and scipy.expm.
import numpy as np
from scipy.linalg import expm

def matrix_exp_eig(Q, t):
    vals, V = np.linalg.eig(Q)            # eigenvalues + right eigenvectors
    D = np.diag(np.exp(vals * t))      # exponentiate the EIGENVALUES
    return (V @ D @ np.linalg.inv(V)).real  # recompose

def matrix_exp_series(Q, t, terms=60):
    Qt = Q * t; P = np.eye(Q.shape[0]); term = np.eye(Q.shape[0])
    for k in range(1, terms):
        term = term @ Qt / k             # (Qt)^k / k!
        P = P + term
    return P

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
for P in (matrix_exp_eig(Q, 3), matrix_exp_series(Q, 3), expm(Q * 3)):
    print(np.round(P[0], 4))   # all -> [0.8379 0.1621], rows sum to 1

The library one-liner: P = scipy.linalg.expm(Q*t). It uses a numerically robust scaling-and-squaring method — in practice you always call this, and the eigendecomposition is for understanding why it works.

From-scratch: decompose P(t) into its modes

To make the “constant plus decay” structure tangible, compute P(t) at several times and watch the transient term shrink while the constant term holds.

python
# Track the two modes of P(t) separately: equilibrium (lambda=0) + transient.
import numpy as np
from scipy.linalg import expm

a, b = 0.2, 1.0
Q = np.array([[-a, a], [b, -b]])
eq = b / (a + b)                       # equilibrium P(UP->UP) = 0.8333

for t in [0.0, 1.0, 3.0, 10.0]:
    transient = (a / (a + b)) * np.exp(-(a + b) * t)   # decaying part
    closed = eq + transient                                     # closed-form P(UP->UP)
    truth = expm(Q * t)[0, 0]                                # scipy ground truth
    print(f"t={t:5.1f}: eq={eq:.4f} + transient={transient:.4f} = {closed:.4f}  (expm {truth:.4f})")
# t=  0.0: eq=0.8333 + transient=0.1667 = 1.0000  (expm 1.0000)  start UP
# t=  3.0: eq=0.8333 + transient=0.0046 = 0.8379  (expm 0.8379)  <- the verified number
# t= 10.0: eq=0.8333 + transient=0.0000 = 0.8333  (expm 0.8333)  transient gone

# Eigenvalues confirm the decay rate: 0 (kept) and -(a+b) (dies).
print("eigenvalues:", np.round(np.linalg.eigvals(Q), 3))   # [0.  -1.2]

The table is the whole chapter in one printout: at every t the answer is the fixed equilibrium 0.8333 plus a transient that decays at rate a+b = 1.2. By t = 10 the transient has vanished and P(UP→UP) sits exactly at the stationary value — the bridge into Chapter 8.

How scipy actually does it: scaling and squaring

Eigendecomposition is great for understanding but fragile numerically (it fails for non-diagonalizable Q and is sensitive to near-equal eigenvalues). Production code uses scaling and squaring, which exploits the identity exp(Qt) = (exp(Qt/2m))2m. Pick m large enough that Qt/2m is tiny, approximate exp of that tiny matrix by a short Taylor (or Padé) series where it is very accurate, then square the result m times to undo the scaling.

The intuition: the power series converges fast only for small matrices, so shrink the matrix, exponentiate cheaply, and recover the full exponent by repeated squaring (each squaring doubles the effective t). It is the same “halve the problem, then recombine” trick as fast exponentiation of numbers, lifted to matrices.

python
# Scaling-and-squaring from scratch: shrink, Taylor, then square back.
import numpy as np
from scipy.linalg import expm

def expm_scale_square(M, m=10, terms=8):
    S = M / (2 ** m)                 # 1) scale down so S is tiny
    E = np.eye(M.shape[0]); term = np.eye(M.shape[0])
    for k in range(1, terms):
        term = term @ S / k          # 2) short Taylor: accurate for tiny S
        E = E + term
    for _ in range(m):
        E = E @ E                  # 3) square m times -> exp(M)
    return E

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
mine = expm_scale_square(Q * 3)
print(np.round(mine[0], 4))            # [0.8379 0.1621]
print(np.allclose(mine, expm(Q * 3)))  # True -> matches scipy

Our hand-rolled scaling-and-squaring matches scipy to machine precision and produces the verified [0.8379, 0.1621]. This is the algorithm behind the one-liner — robust even when Q cannot be diagonalized (repeated eigenvalues), which is exactly when the eigen-method breaks.

Worked example: the full eigendecomposition by hand

Let us assemble exp(3Q) the eigen way, entry by entry, so the recipe is concrete. Eigenvalues are λ1 = 0 and λ2 = −1.2. Right eigenvectors: for λ1 = 0, solve Qv = 0 → v1 = [1, 1]T (all-ones). For λ2 = −1.2, solve (Q + 1.2I)v = 0 → v2 = [1, −6]T (since (−0.2+1.2)·1 + 0.2·(−6) = 1 − 1.2 = −0.2... scaled; the eigenvector direction is what matters).

Stack them as columns of V, build D = diag(e0·3, e−1.2·3) = diag(1, 0.02732), and form exp(3Q) = V·D·V−1. Reading the (UP,UP) entry reproduces 0.8333·1 + 0.1667·0.02732 = 0.8379 — the constant-plus-decay structure, now seen as “weight on the λ=0 mode plus weight on the decaying λ=−1.2 mode.” The eigenvectors are the modes; the eλt are how each mode fades.

Reading all four entries of P(t)

The full transition matrix at t = 3 is P(3) = [[0.8379, 0.1621], [0.8104, 0.1896]] (rows = start UP, start DOWN). Row UP: starting up, 83.79% up / 16.21% down after 3 hr. Row DOWN: starting down, 81.04% up / 18.96% down after 3 hr — the high repair rate has already pulled most of the down-starters back up.

Notice the two rows are already close (0.8379 vs 0.8104), and they converge to the same row [0.8333, 0.1667] as t→∞ — that common row is the stationary π. The whole matrix is collapsing toward a rank-one matrix whose every row is π, which is the geometric meaning of “the chain forgets its start.” Each row obeys the same constant-plus-decay law with different transient weights.

python
# The full P(t) matrix and its collapse toward the rank-one pi-matrix.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
for t in [3, 10, 50]:
    P = expm(Q * t)
    print(f"t={t}:")
    print(np.round(P, 4))
    print("  row gap:", round(np.abs(P[0] - P[1]).max(), 5))   # rows converge
# t=3:  rows [0.8379 0.1621] / [0.8104 0.1896], gap 0.0275
# t=50: rows [0.8333 0.1667] / [0.8333 0.1667], gap ~0  -> rank-one pi-matrix

The “row gap” shrinks to zero as t grows: both rows become π. This is exp(Qt) literally forgetting which state you started in — the transient modes die and only the shared equilibrium row remains.

Worked example: starting DOWN instead

The closed form for the other diagonal is PDOWN→DOWN(t) = a/(a+b) + b/(a+b)·e−(a+b)t — same structure, swapped coefficients. Plug a = 0.2, b = 1.0 at t = 3: constant a/(a+b) = 0.2/1.2 = 0.1667; transient coefficient b/(a+b) = 1.0/1.2 = 0.8333; e−3.6 = 0.02732.

So PDOWN→DOWN(3) = 0.1667 + 0.8333×0.02732 = 0.1667 + 0.02277 = 0.1895, and therefore PDOWN→UP(3) = 1 − 0.1895 = 0.8105. This matches the row-DOWN of the matrix we printed ([0.8104, 0.1896]) to the rounding. The high repair rate means a machine that starts DOWN is already 81% likely to be UP after 3 hr — it bounces back fast, exactly as the transient’s large coefficient 0.8333 says.

From-scratch: how many series terms do you need?

python
# Watch the power series Sum (Qt)^k/k! converge to expm; count terms needed.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
t = 3.0; truth = expm(Q * t)
Qt = Q * t
P = np.eye(2); term = np.eye(2)
for k in range(1, 25):
    term = term @ Qt / k          # (Qt)^k / k!
    P = P + term
    err = np.abs(P - truth).max()
    if k in (2, 5, 10, 15):
        print(f"after {k:2d} terms: P[0]={np.round(P[0],4)}  max err={err:.2e}")
# after  2 terms: P[0]=[0.83  0.17 ]  max err=2e-01   (crude)
# after  5 terms: P[0]=[0.8377 0.1623] max err=4e-03
# after 15 terms: P[0]=[0.8379 0.1621] max err=2e-12  (converged)

The raw power series converges to the verified [0.8379, 0.1621] in about 15 terms here, but for a Q with large entries it needs many more — which is exactly why production code uses scaling-and-squaring instead of summing the series directly. Understanding the series, then trusting expm: the chapter in one workflow.

Why every transient decays: eigenvalues have non-positive real part

We claimed the non-zero eigenvalues of a generator have negative real parts, so their eλt decay. The reason is the Gershgorin circle theorem: every eigenvalue lies within a disk centered at a diagonal entry qii with radius equal to the off-diagonal row sum Σj≠i|qij|. For a generator, that radius is exactly −qii = λi (rows sum to zero), so the disk is centered at −λi with radius λi — it touches the imaginary axis at 0 and lies entirely in the left half-plane.

So no eigenvalue can have positive real part: every mode either stays put (the λ=0 equilibrium) or decays (Re λ < 0). This is the algebraic guarantee that exp(Qt) never blows up and always settles toward π — a CTMC cannot have a runaway growing mode, because conservation forbids it.

python
# Confirm: a generator's eigenvalues sit in the closed left half-plane,
# with exactly one at 0. Test on a random valid 4-state generator.
import numpy as np

def random_generator(n):
    Q = np.random.rand(n, n)        # nonnegative off-diagonals
    np.fill_diagonal(Q, 0.0)
    for i in range(n): Q[i, i] = -Q[i].sum()   # diagonal balances the row
    return Q

for _ in range(3):
    Q = random_generator(4)
    vals = np.linalg.eigvals(Q)
    max_re = max(v.real for v in vals)
    n_zero = sum(abs(v) < 1e-9 for v in vals)
    print(f"max Re(eig)={max_re:+.6f}  zeros={n_zero}")
# max Re(eig)=+0.000000  zeros=1  (every time: one 0, rest in left half-plane)

Across random generators the largest eigenvalue real part is always 0 (the equilibrium mode) and everything else is strictly negative — the decay guarantee, verified. This is why P(t) always converges and never explodes: the generator’s spectrum cannot escape the left half-plane.

Complex eigenvalues: oscillating approach

A chain with three or more states can have complex eigenvalues (in conjugate pairs). Their real part still controls decay, but the imaginary part adds oscillation: e(α + iω)t = eαt(cos ωt + i sin ωt). The eαt (with α < 0) shrinks the amplitude while ω makes the approach to equilibrium spiral rather than settle monotonically.

Physically this is a chain with a preferred direction of circulation — like the one-way ring A→B→C→A — where probability sloshes around the loop with a damped wobble before settling on π. The 2-state chains in this lesson have purely real eigenvalues (0 and −(a+b)), so they approach π smoothly with no overshoot. But knowing the complex case exists explains why some real systems ring before they settle: the generator’s eigenvalues have an imaginary part.

The reassuring constant across all cases: real parts are never positive, so the amplitude never grows. A CTMC can spiral toward equilibrium, but it can never spiral away — conservation, encoded in the row-sum-zero structure of Q, forbids any growing mode. Equilibrium is always reached; only the style of approach (monotone vs oscillating) varies.

python
# exp(Qt) is a stochastic matrix at every t: rows sum to 1, entries in [0,1].
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
for t in [0.5, 3.0, 50.0]:
    P = expm(Q * t)
    ok_rows = np.allclose(P.sum(1), 1)        # conservation
    ok_nonneg = (P >= -1e-12).all()           # valid probabilities
    print(f"t={t:4.1f}: rows-sum-1={ok_rows} nonneg={ok_nonneg}  P[0]={np.round(P[0],4)}")
# Every t: rows sum to 1, all entries >= 0 -- exp(Qt) is always a transition
# matrix, which entrywise exponentiation of Q would VIOLATE (the misconception).
entrywise = np.exp(Q * 3)            # the WRONG thing
print("entrywise exp row sums (broken):", np.round(entrywise.sum(1), 3))  # NOT 1

# And the true expm vs entrywise differ wildly -- they are not the same object.
print("true expm row 0     :", np.round(expm(Q * 3)[0], 4))   # [0.8379 0.1621]
print("entrywise exp row 0 :", np.round(entrywise[0], 4))       # [0.5488 1.8221] nonsense
# entrywise[0,1] = e^0.6 = 1.82 > 1 -- not even a probability. Never do this.
# The matrix exponential mixes the states through its cross terms (Qt)^k;
# exponentiating entries independently throws away exactly that coupling.
diff = np.abs(expm(Q * 3) - entrywise).max()
print("max difference true-vs-entrywise:", round(diff, 3))   # huge -> totally different
# Rule: matrix exponential = power series of (Qt); element-wise exp is a
# different (and meaningless-for-CTMCs) operation. They share only the name "exp".
# In NumPy: scipy.linalg.expm(Q*t) is right; np.exp(Q*t) is the WRONG element-wise one.
# The hand method (eigendecompose, exponentiate eigenvalues, recompose) agrees with expm.
# For repeated eigenvalues use expm directly -- it never needs diagonalizability.
# Bottom line: understand via eigenvalues, compute via scipy.linalg.expm.
# The eigenvalues also reveal the decay rates (transients) and the kept 0-mode (pi).

At every t the true exp(Qt) is a valid stochastic matrix (rows sum to 1, entries nonnegative), while entrywise np.exp(Q*3) produces rows that do not sum to 1 and even contains entries above 1 — concrete proof that the matrix exponential is the power series, not element-wise exponentiation. The conservation that the eigenvalue-0 mode guarantees shows up as “rows always sum to 1.”

When eigendecomposition fails

Not every generator is diagonalizable. A chain with a repeated eigenvalue and too few independent eigenvectors (a defective matrix) cannot be written V·D·V−1. The matrix exponential still exists — the power series always converges — but you need the Jordan form, where exp(Qt) picks up polynomial-in-t factors like t·eλt alongside the pure exponentials. Scaling-and-squaring sidesteps this entirely, which is why it, not eigendecomposition, is the library default. For the 2-state chains in this lesson the eigenvalues are distinct, so the clean eigen-formula always applies.

Chapter 8: The Stationary Distribution π

Run the up/down machine for a long time and the question “is it up right now?” stops depending on how it started. It just hovers at some fixed availability — for our rates, about 83% up. That long-run fraction of time in each state is the stationary distribution π, the resting heartbeat of the system.

We saw π already, hiding in Chapter 7: as t→∞ the decaying eigen-modes vanish and only the λ=0 mode survives, leaving every row of P(t) equal to the same vector — that vector is π. So we could find π by taking t large in exp(Qt), but there is a far cheaper way: solve one balance equation.

The defining equation is πQ = 0 with Σπ = 1. The first part says “the net rate of probability flow into every state is zero” (a row-vector times Q being zero is the steady-state condition); the second part says π is a genuine distribution. Together they pin down a unique π for a well-behaved chain.

For a 2-state chain there is an even more transparent form: detailed balance of flows. The probability flow UP→DOWN is πUP·a and the flow DOWN→UP is πDOWN·b. At equilibrium these must be equal — whatever leaves UP for DOWN must return — so πUP·a = πDOWN·b. That one equation plus normalization gives π instantly.

The result πUP = b/(a+b) has a memorable reading: you spend more time in states you return to faster. A big repair rate b pulls you back to UP quickly and keeps you there; a big fail rate a drags you to DOWN. The fraction of time UP is the return-rate’s share of the total rate.

One conceptual guard: stationary does not mean frozen. The chain keeps jumping forever — UP→DOWN→UP→… without end. What is constant is the probability of finding it in each state. Like a busy office whose per-room headcount is steady though people keep walking around, πQ=0 means net flow is zero, not that motion stops.

Where πQ = 0 comes from

The defining equation is not arbitrary — it is just “the distribution stops changing.” Let p(t) be the row-vector distribution; Chapter 6 says it evolves by dp/dt = p·Q. Stationary means the distribution is constant in time, so its derivative is zero: dp/dt = 0. Substituting gives π·Q = 0 directly. The balance equation is the steady-state of the flow equation, full stop.

This πQ = 0 is the left null-vector of Q, the partner of the all-ones right null-vector from Chapter 5. Q has a zero eigenvalue; its right-eigenvector is 1 (conservation), its left-eigenvector is π (equilibrium). Two sides of the same λ=0 mode.

Global balance vs detailed balance

There are two strengths of balance, and the difference matters. Global balance (πQ = 0) only requires that total inflow to each state equals total outflow — flows may circulate. Detailed balance requires the stronger pairwise condition πiqij = πjqji for every edge: each link carries equal traffic both ways.

Every 2-state chain satisfies detailed balance automatically (there is only one edge, and global = detailed). But a 3-state chain with a one-way cycle A→B→C→A does not: probability circulates around the loop, so πAqAB ≠ πBqBA (indeed qBA = 0). Yet global balance still holds — total in equals total out at each state — so a unique π exists. Detailed balance is sufficient for stationarity but not necessary; chains with it are called reversible.

Worked example: availability of the UP/DOWN machine

Rates a = 0.2 (fail) and b = 1.0 (repair). Balance the flows: πUP·a = πDOWN·b, i.e. 0.2·πUP = 1.0·πDOWN.

Use normalization πUP + πDOWN = 1, so πDOWN = 1 − πUP. Substitute: 0.2·πUP = 1.0·(1 − πUP).

Expand: 0.2·πUP = 1 − πUP. Collect: 0.2·πUP + πUP = 1, i.e. 1.2·πUP = 1.

Solve: πUP = 1/1.2 = 0.8333, and πDOWN = 1 − 0.8333 = 0.1667. The machine is up 83.33% of the time — exactly the constant term b/(a+b) = 1.0/1.2 we saw inside P(t) as t→∞. The general formula confirms it: πUP = b/(a+b) = 1.0/1.2 = 0.8333.

Verify πQ = 0. [0.8333, 0.1667]·Q: first entry 0.8333·(−0.2) + 0.1667·(1.0) = −0.1667 + 0.1667 = 0; second entry 0.8333·(0.2) + 0.1667·(−1.0) = 0.1667 − 0.1667 = 0. So πQ = [0, 0] ✓.

Worked example: a 3-state cycle, global but not detailed balance

Consider a one-way ring: A→B, B→C, C→A, each at rate 1, with no reverse edges. By symmetry the stationary distribution is π = [1/3, 1/3, 1/3]. Verify global balance at A: inflow is from C at rate πC·1 = 1/3; outflow is to B at rate πA·1 = 1/3. Equal — global balance holds, so π is stationary.

But detailed balance fails outright: the A→B flow is πAqAB = 1/3, while the B→A flow is πBqBA = (1/3)·0 = 0 — not equal. Probability genuinely circulates around the ring; there is a net current. This chain is not reversible, yet it still has a perfectly good unique stationary distribution. The lesson: πQ = 0 (global) is the real requirement; detailed balance is a bonus that some chains have and many do not.

Why π is unique (when the chain is irreducible)

Chapter 7 showed P(t) tends to a matrix whose every row is the same vector. That sameness is exactly uniqueness: no matter where you start, you land on one π. Algebraically, an irreducible generator has a simple zero eigenvalue (multiplicity one), so its left null-space is one-dimensional — a single direction, pinned to a probability vector by the normalization Σπ = 1. If the chain is reducible (Chapter 10), the zero eigenvalue repeats and there are multiple stationary distributions, one per island.

The quiz, worked

For a 2-state chain with up-rate (fail) a and down-rate (repair) b, the fraction of time UP at equilibrium is πUP = b/(a+b) — the return rate over the total. You spend more time in the state you come back to faster.

python
# Solve pi Q = 0 with sum(pi)=1 by appending a normalization row, then
# cross-check against the lambda=0 left-eigenvector of Q.
import numpy as np

def stationary(Q):
    Q = np.asarray(Q, float); n = Q.shape[0]
    A = np.vstack([Q.T, np.ones(n)])   # [Q^T ; 1...1]  (pi Q = 0 -> Q^T pi^T = 0)
    b = np.append(np.zeros(n), 1.0)   # [0...0 , 1] (normalization)
    pi, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
    return pi

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
print(stationary(Q))   # [0.8333 0.1667]

from scipy.linalg import expm
print(expm(Q * 100)[0])  # [0.8333 0.1667] -> long-time row of P(t) equals pi

The library one-liner: solve the same augmented system with np.linalg.lstsq as above, or read off the left null-space of Q via scipy.linalg.null_space(Q.T) and normalize. Both give the same π.

How fast does it converge? The spectral gap

Chapter 7 told us P(t) is equilibrium plus decaying transients eλt for the negative eigenvalues. The slowest transient — the eigenvalue closest to 0 — sets the convergence speed. Its magnitude is the spectral gap. For our chain the only nonzero eigenvalue is −(a+b) = −1.2, so deviations from π shrink like e−1.2t.

A useful summary number is the relaxation time 1/|λ2| = 1/1.2 = 0.833 hr — the time for the gap to π to shrink by a factor of e. After a few relaxation times, the chain has “forgotten” its start. Bigger total rate a+b means faster mixing; a chain with tiny rates lingers near its initial distribution for a long time before settling onto π.

From-scratch: solve, verify, and time the convergence

python
# Three checks: solve pi, confirm pi Q = 0, and measure relaxation to pi.
import numpy as np
from scipy.linalg import expm

a, b = 0.2, 1.0
Q = np.array([[-a, a], [b, -b]])
pi = np.array([b, a]) / (a + b)         # [0.8333, 0.1667]
print("pi      =", np.round(pi, 4))
print("pi @ Q  =", np.round(pi @ Q, 6))    # [0. 0.] -> balance confirmed

# Start far from pi and watch the gap decay like exp(-(a+b) t).
p0 = np.array([0.0, 1.0])                  # start fully DOWN
for t in [0.0, 0.83, 1.67, 5.0]:
    pt = p0 @ expm(Q * t)               # distribution at time t
    gap = np.abs(pt - pi).max()
    print(f"t={t:5.2f}: p={np.round(pt,4)}  gap to pi={gap:.4f}")
# t= 0.00: p=[0.    1.   ]  gap to pi=0.8333  (starts far)
# t= 0.83: p=[0.52  0.48 ]  gap to pi=0.31    (~1 relaxation time: gap/e)
# t= 5.00: p=[0.832 0.168]  gap to pi=0.0014  (essentially settled)

eigs = np.linalg.eigvals(Q)
print("relaxation time =", round(1 / abs(sorted(eigs)[0]), 3))  # 0.833 hr

The gap-to-π column shrinks by about a factor of e every 0.83 hr — exactly the relaxation time the spectral gap predicts. After ~5 hr (six relaxation times) the chain is within 0.0014 of equilibrium regardless of starting fully DOWN. Convergence speed is an eigenvalue, not a mystery.

Worked example: detailed balance in two numbers

Verify reversibility of the 2-state chain directly. Detailed balance asks πUP·a = πDOWN·b. Left: 0.8333×0.2 = 0.16667. Right: 0.1667×1.0 = 0.16667. Equal — the UP→DOWN traffic exactly matches the DOWN→UP traffic, so the 2-state chain is reversible. (Every 2-state chain is; there is only one edge to balance.) This pairwise balance is stronger than the global πQ=0, and when it holds it gives an instant shortcut to π: just set the cross-flows equal and normalize.

From-scratch: three independent routes to π

python
# pi via: (1) closed form, (2) power iteration on exp(Q dt), (3) Monte Carlo.
import numpy as np, math, random
from scipy.linalg import expm

a, b = 0.2, 1.0
Q = np.array([[-a, a], [b, -b]])

# (1) closed form
print("closed form  :", np.round(np.array([b, a]) / (a + b), 4))   # [0.8333 0.1667]

# (2) power iteration: repeatedly push a distribution through a short step
A = expm(Q * 0.5)              # one-step matrix for dt=0.5
p = np.array([1.0, 0.0])
for _ in range(200): p = p @ A   # converges to the left eigenvector
print("power iter   :", np.round(p, 4))         # [0.8333 0.1667]

# (3) Monte Carlo: time-weighted occupancy from a long Gillespie run
t, s, time_in = 0.0, 0, [0.0, 0.0]
while t < 200000:
    lam = a if s == 0 else b
    dwell = -math.log(random.random()) / lam
    time_in[s] += dwell; t += dwell; s = 1 - s   # 2-state: flip
total = sum(time_in)
print("monte carlo  :", np.round([x / total for x in time_in], 3))  # ~[0.833 0.167]

All three land on [0.8333, 0.1667]: the closed form, the power iteration (which is just “run P(dt) until it forgets the start”), and the time-weighted Gillespie occupancy. The agreement is the punchline of the chapter — π is simultaneously an algebraic null-vector, the limit of P(t), and the long-run fraction of time spent in each state.

Worked example: solving π for the 3-state machine by hand

Take the 3-state generator from Chapter 5: Q = [[−0.25, 0.2, 0.05],[1.0, −1.0, 0],[0.5, 0, −0.5]] (UP, DOWN, MAINT). Write πQ = 0 column by column. Column UP: −0.25πU + 1.0πD + 0.5πM = 0. Column DOWN: 0.2πU − 1.0πD = 0 → πD = 0.2πU. Column MAINT: 0.05πU − 0.5πM = 0 → πM = 0.1πU.

Now normalize: πU + 0.2πU + 0.1πU = 1 → 1.3πU = 1 → πU = 0.7692. Then πD = 0.2×0.7692 = 0.1538 and πM = 0.1×0.7692 = 0.0769. So π = [0.769, 0.154, 0.077] — the machine is up ~77% of the time once a maintenance state competes for it, down from 83% in the 2-state model. Adding an exit from UP necessarily lowers UP’s share.

Verify column UP as a check: −0.25×0.7692 + 1.0×0.1538 + 0.5×0.0769 = −0.1923 + 0.1538 + 0.0385 = 0 ✓. The balance holds, confirming the hand solve.

python
# 3-state stationary: solve pi Q = 0 with normalization, verify the hand answer.
import numpy as np

Q = np.array([[-0.25, 0.2, 0.05],
              [ 1.0, -1.0, 0.0],
              [ 0.5, 0.0, -0.5]])
n = 3
A = np.vstack([Q.T, np.ones(n)])      # pi Q = 0  AND  sum pi = 1
b = np.append(np.zeros(n), 1.0)
pi = np.linalg.lstsq(A, b, rcond=None)[0]
print("pi      =", np.round(pi, 4))      # [0.7692 0.1538 0.0769]
print("pi @ Q  =", np.round(pi @ Q, 6))  # [0 0 0] -> balance confirmed
print("UP share vs 2-state 0.8333:", round(pi[0], 3))   # 0.769 < 0.833

The solver confirms [0.769, 0.154, 0.077] and πQ = 0. The drop in UP-share from 0.833 to 0.769 is the quantitative cost of the maintenance exit — every new way to leave UP eats into the time spent there, exactly as the balance equations dictate.

Sensitivity: how π moves when a rate moves

A practical question: if you double the repair rate b, how much availability do you gain? With πUP = b/(a+b), push b from 1.0 to 2.0 (keeping a = 0.2): availability rises from 1.0/1.2 = 0.833 to 2.0/2.2 = 0.909 — a 7.6-point gain. Halve a instead (fail less often, 0.2→0.1): availability rises to 1.0/1.1 = 0.909 too. Reducing failures and speeding repairs are interchangeable levers here, both worth ~7.6 points.

The formula makes the trade explicit: availability depends only on the ratio b/a (rewrite b/(a+b) = 1/(1 + a/b)). Whether you halve the fail rate or double the repair rate, you double b/a and land on the same availability. That is the kind of design insight a stationary distribution hands you for free — one formula, every what-if.

python
# Availability sensitivity: vary fail rate a and repair rate b, see b/a drive it.
def availability(a, b): return b / (a + b)

base = availability(0.2, 1.0)
print(f"baseline (a=0.2,b=1.0): {base:.4f}")             # 0.8333
print(f"double repair (b=2.0) : {availability(0.2,2.0):.4f}")  # 0.9091
print(f"halve fail   (a=0.1)  : {availability(0.1,1.0):.4f}")  # 0.9091  (same!)

# Both moves double the ratio b/a -> identical availability.
for a, b in [(0.2, 1.0), (0.1, 1.0), (0.2, 2.0), (0.05, 1.0)]:
    print(f"a={a}, b={b}: ratio b/a={b/a:5.1f}  availability={availability(a,b):.4f}")
# ratio 5.0 -> 0.833; ratio 10.0 -> 0.909; ratio 20.0 -> 0.952
# availability is a function of b/a alone -> the design knob is the RATIO.

The table shows availability tracking the ratio b/a exactly: ratio 5 gives 0.833, ratio 10 gives 0.909, ratio 20 gives 0.952. Reliability engineering lives in this one curve — spend your budget wherever it moves b/a the most, because that, and only that, sets the steady-state uptime.

From-scratch: start-independence, four random starts

The defining promise of an irreducible chain is that any start converges to the same π. Let us demonstrate it the most convincing way: drop four wildly different initial distributions and watch them all funnel to [0.8333, 0.1667].

python
# Four different starts, one destination: start-independence of pi.
import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
starts = [[1.0, 0.0], [0.0, 1.0], [0.5, 0.5], [0.9, 0.1]]

for t in [1, 5, 20]:
    Pt = expm(Q * t)
    rows = [np.round(np.array(s) @ Pt, 3) for s in starts]
    spread = max(r[0] for r in rows) - min(r[0] for r in rows)
    print(f"t={t:2d}: UP-shares {[r[0] for r in rows]}  spread={spread:.3f}")
# t= 1: UP-shares [0.834, 0.747, 0.79, 0.825]  spread=0.087  (still differ)
# t= 5: UP-shares [0.834, 0.83, 0.832, 0.833]   spread=0.004  (converging)
# t=20: UP-shares [0.833, 0.833, 0.833, 0.833]  spread=0.000  (all -> pi)

The “spread” between the four starts collapses from 0.087 at t=1 to 0.000 by t=20: every initial distribution forgets itself and lands on the same πUP = 0.833. That collapse is the chain mixing — the transient modes dying at rate a+b = 1.2 — and it is what makes π a property of the chain, not of where you happened to start.

python
# Detailed-balance test: is the chain reversible? (pi_i q_ij == pi_j q_ji)
import numpy as np

def reversible(Q, pi):
    n = len(pi)
    for i in range(n):
        for j in range(n):
            if i != j and not np.isclose(pi[i] * Q[i, j], pi[j] * Q[j, i]):
                return False
    return True

# 2-state chain: reversible (only one edge to balance)
Q2 = np.array([[-0.2, 0.2], [1.0, -1.0]])
pi2 = np.array([1.0, 0.2]); pi2 /= pi2.sum()
print("2-state reversible:", reversible(Q2, pi2))    # True

# 3-state one-way ring A->B->C->A: NOT reversible (circulating current)
Q3 = np.array([[-1., 1., 0.], [0., -1., 1.], [1., 0., -1.]])
pi3 = np.array([1, 1, 1]) / 3
print("3-ring reversible :", reversible(Q3, pi3))    # False (but pi3 still stationary)
print("3-ring pi Q = 0   :", np.allclose(pi3 @ Q3, 0))   # True -- global balance holds
# net A->B current = pi_A*q_AB - pi_B*q_BA = 1/3 - 0 = 1/3 != 0 -> irreversible
print("net A->B current  :", round(pi3[0] * Q3[0, 1] - pi3[1] * Q3[1, 0], 3))   # 0.333
# A nonzero current means probability circulates A->B->C->A forever, even
# though each state's occupancy (1/3) is constant. Stationary, not static.
print("reversible chains have ZERO net current on every edge (detailed balance).")
# so: pi Q = 0 always; detailed balance (zero current) is the stronger bonus.

The 2-state chain passes the reversibility test; the one-way ring fails it (probability circulates, so no edge balances pairwise) yet still satisfies global balance πQ = 0. This is the chapter’s key distinction in code: detailed balance is a special property, but πQ = 0 is the universal requirement for stationarity.

Chapter 9: Showcase — The Living CTMC (Gillespie)

This is the payoff. One generator Q, four windows onto it, all wired together. Move a rate and watch the holding-time histograms stretch, the P(t) curves rebend, individual sample lives jitter, and the long-run π split shift — the whole theory breathing at once.

The engine driving the live sample path is the Gillespie algorithm — the exact, event-driven simulator we promised back in Chapter 0. It never uses a fixed tick. Instead it leaps directly from one jump to the next, spending compute only where the state actually changes.

Each Gillespie step is exactly Chapter 4’s reassemble recipe, two random draws. First, the when: from the current state i, the total rate is λi = −qii, so draw a holding time −ln(U)/λi from a uniform U. Advance the clock by that amount. Second, the where: draw the next state from the jump-chain row qij/λi. Move there and repeat.

Watch the three agreements form as you let it run. The per-state holding-time histograms fill in and match the exponential(λ_i) density. The empirical fraction of time spent in each state creeps toward the exact π from Chapter 8. And the P(t) curves you can read at any cursor time are the exact exp(Qt) from Chapter 7. Theory and simulation, converging before your eyes.

Try to break it. Toggle to 3 states for a richer chain. Crank the repair rate b high and watch π pile onto UP while the DOWN dwell bars shrink. Make a rate tiny and watch its dwell histogram stretch far to the right. Every knob moves all four panels coherently because they are all faces of the single matrix Q.

There is no quiz here — the simulator is the test. If you can predict, before you drag a slider, which way each of the four panels will move, you understand CTMCs. That is the whole lesson, made tactile.

Worked example: one Gillespie step by hand (3-state)

Start in state A with rates qAB = 1, qAC = 3, so λA = 4.

Draw the when. Sample a uniform U = 0.5. Holding time = −ln(U)/λA = −ln(0.5)/4 = 0.6931/4 = 0.1733. Advance the clock to t = 0.1733.

Draw the where. The jump-chain row from A is [PAB, PAC] = [1/4, 3/4] = [0.25, 0.75], with cumulative thresholds [B: 0.25, C: 1.0]. Sample a second uniform U2 = 0.6. Since 0.6 > 0.25, we fall past B into the C interval — jump to C.

Result: at t = 0.1733 we move A→C. Now we are in C; repeat with C’s rates. Two random numbers per jump — one exponential draw for the time, one uniform draw for the destination — exactly the race of Chapter 3, sampled.

python
# The full Gillespie loop: leap event-to-event, no fixed tick anywhere.
import numpy as np

def gillespie(Q, s0, T):
    Q = np.asarray(Q, float); n = Q.shape[0]
    lam = -np.diag(Q)
    Pjump = Q / lam[:, None]; np.fill_diagonal(Pjump, 0.0)
    t, s = 0.0, s0; times, states = [0.0], [s0]
    while t < T:
        hold = np.random.exponential(1.0 / lam[s])   # WHEN: -ln(U)/lambda
        t += hold
        s = np.random.choice(n, p=Pjump[s])     # WHERE: q_ij / lambda_i
        times.append(t); states.append(s)
    return np.array(times), np.array(states)

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
times, states = gillespie(Q, 0, 100000)
# time-weighted occupancy converges to pi = [0.833, 0.167]:
dur = np.diff(times); occ = np.array([dur[states[:-1]==k].sum() for k in (0,1)])
print(occ / occ.sum())   # ~[0.833 0.167], matching stationary(Q)

The library one-liner: many ecosystems ship Gillespie directly (e.g. the SSA solvers in gillespy2 or scipy-based chemical-kinetics packages), but the loop above is the entire algorithm — the same engine that powers chemical-kinetics and queueing simulators.

Chapter 10: Edge Cases & When It Breaks

The clean theory assumed you always eventually leave every state and the chain mixes to one π. Real models break these assumptions, and knowing the failure modes is how you decide whether to trust a model at all.

The first failure mode is an absorbing state: a state with no outgoing rates, qii = 0 and all qij = 0. Its total leave-rate is 0, so the holding time is infinite — once entered, never left. A “crashed beyond repair” state is absorbing; every sample path eventually funnels in and stops changing, and the column of P(t) for that state tends to 1.

The second is reducibility: the state space splits into clusters that cannot reach each other. Then there is no single π — the long-run distribution depends on which island you started in. A chain must be irreducible (all states communicate) for the unique-π promise of Chapter 8 to hold.

The third is explosion: if rates grow without bound as you move through states, holding times shrink toward zero so fast that infinitely many jumps fit into a finite time. A finite-state CTMC can never explode (its rates are bounded), but infinite-state models (some queues, branching processes) can, and then P(t) is not even well-defined past the explosion time.

For absorbing chains the interesting question changes from “what is π?” (trivially, all mass on the absorbing state) to “how long until I get absorbed?” That is the expected time to absorption, and it has a clean linear-system answer. Collect the transient states (the non-absorbing ones) and let QT be the sub-generator among them. The expected hitting times τ solve −QT·τ = 1 (the all-ones vector).

For a simple chain on a single path, this collapses to something you can see: the expected time to absorption is just the sum of the mean holding times along the path. Sit in each transient state for its average dwell, then move on, until you fall into the trap. We will work exactly that case below.

Deriving the absorption-time system −Q_T·τ = 1

Why does −QT·τ = 1 give the expected time to absorption? Condition on the first jump (a “first-step” argument). Starting in transient state i, you wait a mean holding time 1/λi, then jump to some state j with probability Pij = qij/λi. From j you still expect τj more time (zero if j is absorbing).

Multiply through by λi: λiτi = 1 + Σj qijτj. Move the sum left: λiτi − Σj qijτj = 1. Since −λi = qii (the transient diagonal), the left side is exactly −(QT·τ)i. Stacking over all transient i gives the matrix equation −QT·τ = 1 — derived from one first-step recursion, no magic.

For a single-path chain this matrix system collapses to the sum-of-dwells shortcut, because each transient state has exactly one exit and the recursion unrolls into τi = 1/λi + τnext. The general solve handles branching paths, where you cannot just add dwells because the route is uncertain.

Worked example: time to failure of WORK→DEG→FAILED

A reliability chain: WORK degrades to DEG at rate λWORK = 0.5, DEG fails to FAILED at rate λDEG = 0.3, and FAILED is absorbing (all its rates are 0). What is the expected time from WORK to FAILED?

Because the path is a single line WORK→DEG→FAILED, the expected absorption time is the sum of the mean holding times in the transient states. Mean stay in WORK: 1/λWORK = 1/0.5 = 2.0 time units. Mean stay in DEG: 1/λDEG = 1/0.3 = 3.3333 time units.

Add them: 2.0 + 3.3333 = 5.3333 time units to failure, on average. That is the answer the linear system −QT·τ = 1 would produce too — here it simplifies to a sum because each transient state has exactly one exit.

And the stationary distribution is degenerate: all probability eventually sits on FAILED, so π = [0, 0, 1] (in order WORK, DEG, FAILED). There is no “equilibrium availability” to speak of — the only steady state is death.

The quiz, worked

A state i has a generator row that is entirely zero (qii = 0 and all qij = 0). Its total leave-rate is 0, so the holding time is infinite: once entered, it is never left — the state is absorbing. Its column of P(t) tends to 1, and if it is reachable, π puts all its mass there.

python
# Detect absorbing states, and solve expected time-to-absorption: -Q_T tau = 1.
import numpy as np

def is_absorbing(Q, i):
    return np.allclose(Q[i], 0.0)          # a fully-zero row = no exits

def time_to_absorption(Q, transient):
    Q = np.asarray(Q, float)
    QT = Q[np.ix_(transient, transient)]     # sub-generator on transient states
    ones = np.ones(len(transient))
    return np.linalg.solve(-QT, ones)        # expected hitting times tau

# WORK=0, DEG=1, FAILED=2 (absorbing)
Q = np.array([[-0.5, 0.5, 0.0],
              [ 0.0,-0.3, 0.3],
              [ 0.0, 0.0, 0.0]])
tau = time_to_absorption(Q, [0, 1])
print(tau[0])          # 5.3333 -> expected time from WORK to FAILED
print(is_absorbing(Q, 2))  # True

The library one-liner is the solve itself: tau = np.linalg.solve(-Q_T, np.ones(k)), where Q_T is the transient block of Q. Everything else is bookkeeping about which states are transient.

Worked example: a reducible chain has no single π

Take four states split into two islands that never communicate: {A, B} with A&rlrhar;B at rates 1 each, and {C, D} with C⇌D at rates 1 and 3. There is no edge between the islands. Within island 1, balance gives π = [0.5, 0.5]; within island 2, π = [3/4, 1/4] (you spend more time in C, which you return to faster).

Now the long-run distribution depends entirely on where you start: begin in A and you converge to [0.5, 0.5, 0, 0]; begin in C and you converge to [0, 0, 0.75, 0.25]. There is no single π for the whole chain — the “forget your start” promise of Chapter 8 fails because the chain cannot mix across the islands. This is reducibility, and it is why the irreducibility check is not optional.

Worked example: variance of time-to-absorption

For the single-path WORK→DEG→FAILED chain, the absorption time is a sum of two independent exponential dwells. Independent variances add: Var = 1/λWORK² + 1/λDEG² = 1/0.25 + 1/0.09 = 4 + 11.11 = 15.11. The standard deviation is √15.11 = 3.89 time units — against a mean of 5.33. So “5.33 give or take ~3.9”: time-to-failure is highly variable, a crucial caveat reliability engineers care about, not just the mean.

From-scratch: absorption time three ways

python
# Mean absorption time: linear-solve, sum-of-dwells, and Monte Carlo agree.
import numpy as np, math, random

Q = np.array([[-0.5, 0.5, 0.0],
              [ 0.0,-0.3, 0.3],
              [ 0.0, 0.0, 0.0]])      # WORK, DEG, FAILED(absorbing)

# 1) linear system -Q_T tau = 1 on transient block {WORK, DEG}
QT = Q[:2, :2]
tau = np.linalg.solve(-QT, np.ones(2))
print("linear solve  :", round(tau[0], 4))    # 5.3333

# 2) sum of mean dwells (works only on a single path)
print("sum of dwells :", round(1 / 0.5 + 1 / 0.3, 4))   # 5.3333

# 3) Monte Carlo: simulate many WORK->FAILED first-passage times
def absorb_time():
    t, s = 0.0, 0
    while s != 2:
        lam = -Q[s, s]
        t += -math.log(random.random()) / lam
        s += 1                       # single path: always advance
    return t
samples = [absorb_time() for _ in range(100000)]
mc = sum(samples) / len(samples)
print("monte carlo   :", round(mc, 3))           # ~5.33
var = sum((x - mc) ** 2 for x in samples) / len(samples)
print("variance      :", round(var, 2))          # ~15.1 -> std ~3.9

All three routes land on 5.33, and the Monte Carlo variance ~15.1 matches the hand calculation — the mean is reliable but the spread is large, exactly the reliability caveat above. The linear solve is the only one that survives when the path branches.

Worked example: branching absorption (which trap?)

When there are two absorbing states, the new question is “which one do I hit?” Suppose from a transient state S you go to trap GOOD at rate 1 and trap BAD at rate 3, with no other exits. This is just Chapter 3’s race between two absorbing destinations: P(hit BAD) = 3/(1+3) = 0.75 and P(hit GOOD) = 1/4 = 0.25. The absorption probabilities are the jump-chain probabilities; the absorption time is the dwell 1/(1+3) = 0.25.

With intermediate transient states the hitting probabilities solve a linear system too: B = (−QT)−1·R, where R holds the rates from each transient state into each absorbing state. The expected-time solve −QT·τ = 1 and the hitting-probability solve −QT·B = R share the same matrix −QT — one factorization answers both “how long” and “where to.”

From-scratch: hitting probabilities and a reducible chain

python
# Absorption probabilities B = (-Q_T)^-1 R, and a reducible-chain demo.
import numpy as np
from scipy.linalg import expm

# States: S=0 (transient), GOOD=1, BAD=2 (both absorbing)
Q = np.array([[-4.0, 1.0, 3.0],
              [ 0.0, 0.0, 0.0],
              [ 0.0, 0.0, 0.0]])
QT = Q[:1, :1]                 # transient block (just S)
R = Q[:1, 1:]                  # rates into the absorbing states
B = np.linalg.solve(-QT, R)    # hitting probabilities
print("P(hit GOOD, BAD) =", np.round(B[0], 3))   # [0.25 0.75]

# Reducible chain: two islands, long-run answer depends on the start.
Qr = np.array([[-1., 1., 0., 0.],
               [ 1.,-1., 0., 0.],     # island {A,B}
               [ 0., 0.,-1., 1.],
               [ 0., 0., 3.,-3.]])    # island {C,D}
print("from A:", np.round(([1,0,0,0] @ expm(Qr * 50)), 3))  # [0.5 0.5 0 0]
print("from C:", np.round(([0,0,1,0] @ expm(Qr * 50)), 3))  # [0 0 0.75 0.25]
# Two different limits -> NO single pi. Irreducibility is required.

The first solve gives the 25/75 split into the two traps; the reducible run shows the long-run distribution genuinely depends on the starting island ([0.5,0.5,0,0] from A vs [0,0,0.75,0.25] from C). Both confirm the chapter’s warnings numerically: absorbing chains funnel to traps, and reducible chains have no single π.

Explosion: when infinitely many jumps fit in finite time

The third failure mode — explosion — cannot happen in a finite chain but is real in infinite-state models. Imagine a chain where state n has leave-rate λn = 2n: each state is held half as long as the previous. The expected time to take infinitely many jumps is Σ 1/2n, a convergent geometric series — so the chain takes infinitely many jumps in finite total time and “explodes” past the boundary. After the explosion time, P(t) is not even defined.

The guard is a divergence condition: as long as Σ 1/λn diverges (holding times do not shrink too fast), the chain is non-explosive and everything in this lesson applies. Finite chains have bounded rates, so the sum trivially diverges — they are always safe. Explosion is a hazard only for unbounded-rate models like certain birth processes, and it is why “finite-state” was a quiet assumption throughout.

python
# Survival of a transient population over time: mass drains into the trap.
import numpy as np
from scipy.linalg import expm

# WORK=0, DEG=1, FAILED=2 (absorbing)
Q = np.array([[-0.5, 0.5, 0.0],
              [ 0.0,-0.3, 0.3],
              [ 0.0, 0.0, 0.0]])
p0 = np.array([1.0, 0.0, 0.0])         # start in WORK
print(f"{'t':>4} {'WORK':>7} {'DEG':>7} {'FAILED':>7}")
for t in [0, 2, 5, 10, 30]:
    p = p0 @ expm(Q * t)
    print(f"{t:>4} {p[0]:>7.3f} {p[1]:>7.3f} {p[2]:>7.3f}")
#   t   WORK    DEG  FAILED
#   0  1.000  0.000  0.000
#   5  0.082  0.213  0.705   (draining into the trap)
#  30  0.000  0.000  1.000   (all mass absorbed -> pi=[0,0,1])

The FAILED column climbs monotonically to 1 while WORK and DEG drain to 0 — the mass funnels into the absorbing state, confirming π = [0, 0, 1]. There is no equilibrium between the transient states; the only steady state is the trap, which is why the interesting question became “how long until absorption,” answered above.

Worked example: a branching reliability chain

Make the chain branch so the sum-of-dwells shortcut fails and the linear system earns its keep. From WORK (rate out 0.5) you go to DEG with probability 0.8 or straight to a SPARE state with probability 0.2; DEG fails at 0.3; SPARE returns to WORK at 0.4. FAILED is the absorbing trap. Now the route is uncertain — you might loop WORK→SPARE→WORK several times before degrading — so you cannot just add two dwells.

The linear system −QT·τ = 1 over the three transient states {WORK, DEG, SPARE} handles the branching automatically: each τi is the mean dwell plus the probability-weighted expected times from wherever you jump. Solving it (below) gives a mean time to FAILED that accounts for every possible looping path — impossible to get by hand-summing dwells.

python
# Branching chain: WORK can loop via SPARE. Only the linear solve is correct.
import numpy as np

# states: WORK=0, DEG=1, SPARE=2, FAILED=3 (absorbing)
# WORK -> DEG at 0.4, WORK -> SPARE at 0.1 (so 0.8/0.2 split, total leave 0.5)
Q = np.array([[-0.5, 0.4, 0.1, 0.0],
              [ 0.0,-0.3, 0.0, 0.3],
              [ 0.4, 0.0,-0.4, 0.0],
              [ 0.0, 0.0, 0.0, 0.0]])
transient = [0, 1, 2]
QT = Q[np.ix_(transient, transient)]
tau = np.linalg.solve(-QT, np.ones(len(transient)))
print("mean time to FAILED from each transient state:")
print("  WORK :", round(tau[0], 3))   # accounts for WORK<->SPARE looping
print("  DEG  :", round(tau[1], 3))   # 1/0.3 = 3.333 (one exit)
print("  SPARE:", round(tau[2], 3))   # dwell + back to WORK
# tau[WORK] > the naive 5.33 because the SPARE loop delays degradation.

The solve gives a WORK-to-FAILED time larger than the single-path 5.33, because the SPARE loop sends the system back to WORK and postpones failure — a delay the sum-of-dwells shortcut simply cannot represent. Whenever the route branches or loops, −QT·τ = 1 is the only correct tool.

The fundamental matrix (−Q_T)⁻¹

The matrix N = (−QT)−1 that appears in both the absorption-time solve and the hitting-probability solve has a beautiful meaning of its own: entry Nij is the expected total time spent in transient state j before absorption, starting from i. The expected absorption time τi = Σj Nij is then just “sum the time spent in every transient state” — which is why τ = N·1 (the row sums of N).

This is the continuous-time analog of the discrete absorbing-chain fundamental matrix. It packages every absorption statistic: row sums give expected times, products with the exit-rate matrix give hitting probabilities, and second moments (for variance) come from N applied twice. One inverse, computed once, answers the entire family of “what happens before the trap” questions — the single most useful object for any absorbing CTMC.

The practical takeaway: when you face an absorbing chain, compute N = (−QT)−1 first, then read off whatever you need. Expected time to failure, probability of each failure mode, time spent in each degraded state — all are simple reads of one matrix, no resimulation required.

python
# The fundamental matrix N answers the whole absorbing-chain family at once.
import numpy as np

# WORK=0, DEG=1, FAILED=2 (absorbing)
Q = np.array([[-0.5, 0.5, 0.0],
              [ 0.0,-0.3, 0.3],
              [ 0.0, 0.0, 0.0]])
QT = Q[:2, :2]
N = np.linalg.inv(-QT)            # N[i,j] = expected time in j starting from i
print("fundamental matrix N:")
print(np.round(N, 3))                   # [[2.    3.333][0.    3.333]]
print("time in WORK from WORK :", round(N[0, 0], 3))   # 2.0  (1/0.5)
print("time in DEG  from WORK :", round(N[0, 1], 3))   # 3.333 (1/0.3)
print("expected absorption (row sum):", round(N[0].sum(), 3))  # 5.333
# One inverse -> time-in-each-state AND total time to absorption.

# Second moments for VARIANCE come from N too: Var(tau) = (2N - I) N 1 - (N1)^2
ones = np.ones(2)
tau = N @ ones
var = ((2 * N - np.eye(2)) @ tau) - tau ** 2
print("Var(absorption time) from WORK:", round(var[0], 2))   # ~15.11 -> std ~3.9
print("std dev:", round(var[0] ** 0.5, 2), "vs mean", round(tau[0], 2))   # 3.89 vs 5.33
# N also gives hitting probabilities: B = N @ R for absorbing-state rates R.
# One inverse (-Q_T)^-1 powers mean, variance, AND which-trap probabilities.
# Compute N once, then read off every "what happens before the trap" statistic.
# This is the continuous-time analog of the discrete absorbing-chain fundamental matrix.
# No resimulation needed -- the inverse already contains every answer.
# For a single-path chain it collapses to the simple sum of mean dwells.
# spread (3.89) is comparable to the mean (5.33) -- time-to-failure is noisy
# Quoting only the mean MTTF hides this: many failures land far from 5.33h.
print("coefficient of variation:", round(var[0] ** 0.5 / tau[0], 2))   # 0.73 -- substantial spread

The fundamental matrix N directly gives the time spent in each transient state (2 hr in WORK, 3.33 hr in DEG starting from WORK), its row sum is the total expected absorption time 5.33, and a second-moment formula on N yields the variance ~15.11 — matching our hand calculation. One matrix inverse, every absorption statistic including the spread.

Chapter 11: Connections — Where CTMCs Live

You now hold the continuous-time core: rates, exponential clocks, the generator Q, P(t)=exp(Qt), the stationary π, and the Gillespie simulator. Step back and see the family — the same machinery reappears, lightly modified, across a surprising range of models.

Add a fixed tick and you are back to a discrete Markov chain — the tick-based cousin we started by criticizing. The relationship runs both ways: sampling a CTMC at regular intervals dt gives a discrete chain whose one-step matrix is exp(Q·dt).

Add noisy observations of a hidden state and you get a Hidden Markov Model (HMM). The underlying dynamics are still a Markov chain; you simply never see the state directly, only emitted symbols, and must infer it with Forward/Viterbi. The continuous-time version (CT-HMM) keeps the generator Q and the gap-aware exp(Q·dt) we built here.

Restrict the jumps to plus-or-minus one on a count and you get a birth-death process — the backbone of queueing theory. States are “number in system,” arrivals push up, departures pull down, and Q is tridiagonal. This is the M/M/1 queue, worked below.

Add actions and rewards and you get a continuous-time Markov Decision Process — the same rate dynamics, now controllable, the foundation for continuous-time reinforcement learning. And bolt a filtering layer on top — estimating the hidden state from a stream of observations — and you are in the territory of the Kalman filter and the general Bayes filter, or its partially-observed control cousin the POMDP.

The thread through all of it is the generator-and-clocks viewpoint. Reliability engineering is up/down CTMCs with more states; chemical kinetics is a CTMC on molecule counts; epidemic models are CTMCs on compartment populations. Q, holding times, and π are a unifying lens, not a side quest.

Worked example: the M/M/1 queue from a CTMC balance

Model a queue as a CTMC on the count 0, 1, 2, … of jobs in system. Customers arrive at rate λ (rate i→i+1) and are served at rate μ (rate i→i−1). The generator Q is tridiagonal. Take λ = 2/min and μ = 3/min.

Define the utilization ρ = λ/μ = 2/3 = 0.6667. Balance between adjacent states (flow up = flow down across the boundary between i and i+1) gives πi+1·μ = πi·λ, i.e. πi+1 = ρ·πi. So the distribution is geometric: πi = ρi·π0.

Normalize: Σi≥0 ρi·π0 = 1. The geometric series sums to π0/(1−ρ) = 1 (valid because ρ < 1), so π0 = 1 − ρ = 1 − 0.6667 = 0.3333.

So the queue is empty one-third of the time. And π1 = ρ·π0 = 0.6667×0.3333 = 0.2222, π2 = ρ²·π0 = 0.4444×0.3333 = 0.1481, and so on — a real M/M/1 result derived straight from a CTMC balance equation, with no new machinery beyond Chapter 8.

Performance metrics from the geometric π

The geometric π hands us the headline queue metrics for free. The average number in system is L = Σi i·πi; for a geometric distribution this sums to L = ρ/(1−ρ) = 0.6667/0.3333 = 2.0 jobs on average. The mean keeps climbing as ρ→1 — at 90% utilization, L = 0.9/0.1 = 9, and it blows up to infinity as the server saturates.

Little’s law then converts a count into a wait: L = λ·W, where W is the average time a job spends in the system. Solve W = L/λ = 2.0/2 = 1.0 min. So at λ=2, μ=3 a job spends on average one minute in the system, two-thirds of which is waiting and one-third in service. These are the numbers a capacity planner actually quotes — all flowing from one tridiagonal generator and Chapter 8’s balance.

The deep lesson: M/M/1 is not a separate theory to memorize. It is literally a CTMC on the integers with ±1 jumps, and every queueing formula is a stationary-distribution calculation in disguise. The generator-and-balance viewpoint you built in Chapter 8 is the whole of elementary queueing theory.

python
# Birth-death (M/M/1) generator + the closed-form geometric stationary check.
import numpy as np

def birth_death_Q(lam, mu, N):
    Q = np.zeros((N + 1, N + 1))
    for i in range(N + 1):
        if i < N: Q[i, i + 1] = lam   # arrival: i -> i+1
        if i > 0: Q[i, i - 1] = mu   # service: i -> i-1
        Q[i, i] = -(Q[i].sum() - Q[i, i])  # diagonal balances the row
    return Q

lam, mu = 2.0, 3.0; rho = lam / mu
pi0 = 1 - rho
print(rho, pi0)                  # 0.6667 0.3333
print([(1 - rho) * rho ** i for i in range(4)])  # [0.333 0.222 0.148 0.099]

The library one-liner ties back to Chapter 8: stationary(birth_death_Q(2,3,50)) reproduces the geometric π numerically — the same πQ=0 solver, applied to a tridiagonal Q.

From-scratch: the full M/M/1 dashboard

python
# M/M/1 end to end: build Q, solve pi, derive L and W, compare to closed forms.
import numpy as np

def birth_death_Q(lam, mu, N):
    Q = np.zeros((N + 1, N + 1))
    for i in range(N + 1):
        if i < N: Q[i, i + 1] = lam
        if i > 0: Q[i, i - 1] = mu
        Q[i, i] = -(Q[i].sum() - Q[i, i])
    return Q

def stationary(Q):
    n = Q.shape[0]
    A = np.vstack([Q.T, np.ones(n)]); b = np.append(np.zeros(n), 1.0)
    return np.linalg.lstsq(A, b, rcond=None)[0]

lam, mu, N = 2.0, 3.0, 80
rho = lam / mu
pi = stationary(birth_death_Q(lam, mu, N))
print("pi[0..3]      :", np.round(pi[:4], 4))      # [0.3333 0.2222 0.1481 0.0988]
print("empty prob    :", round(pi[0], 4), "vs 1-rho", round(1 - rho, 4))

L = sum(i * pi[i] for i in range(len(pi)))
print("L (num)       :", round(L, 3), "vs rho/(1-rho)", round(rho / (1 - rho), 3))  # 2.0
print("W = L/lambda  :", round(L / lam, 3), "min")   # 1.0  (Little's law)

The numerical π reproduces the geometric closed form, L = 2.0 matches ρ/(1−ρ), and Little’s law turns it into a one-minute average wait — the whole queue characterized by the same πQ=0 solver from Chapter 8 applied to a tridiagonal Q. One generator, every metric.

python
# Utilization sweep: how L and W blow up as the server saturates (rho->1).
mu = 3.0
for lam in [1.0, 2.0, 2.7, 2.97]:
    rho = lam / mu
    L = rho / (1 - rho)          # avg number in system
    W = L / lam                    # avg time in system (Little)
    print(f"rho={rho:.2f}: L={L:6.2f} jobs  W={W:6.2f} min")
# rho=0.33: L=  0.50  W=  0.50   (lightly loaded)
# rho=0.67: L=  2.00  W=  1.00   (our worked case)
# rho=0.90: L=  9.00  W=  3.33   (queue building fast)
# rho=0.99: L= 99.00  W= 33.33   (near saturation -> blows up)

The 1/(1−ρ) blow-up as ρ→1 is the single most important fact in queueing: a server at 99% utilization has a hundred jobs waiting on average. It is why capacity planners keep utilization well below 1, and it comes straight out of the geometric π — a CTMC stationary distribution wearing a queueing hat.

The family at a glance

Each row is the CTMC plus one ingredient. Read down the “what you add” column and the whole landscape of stochastic modeling organizes itself around the generator-and-clocks core you just built. Nothing in that column replaces Q; each only layers something on top.

From-scratch: sampling a CTMC gives a discrete Markov chain

The first family relationship — “add a tick and you get a discrete chain” — is checkable in a few lines. Sample the CTMC at a fixed interval dt; the resulting state sequence is a discrete Markov chain whose one-step transition matrix is exactly exp(Q·dt).

python
# Sample a CTMC every dt; empirical one-step matrix matches exp(Q*dt).
import numpy as np, math, random
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [1.0, -1.0]])
lam = -np.diag(Q); dt = 1.0

# run one long Gillespie path, then read its state on a dt grid
t, s, jumps = 0.0, 0, [(0.0, 0)]
while t < 400000:
    t += -math.log(random.random()) / lam[s]; s = 1 - s; jumps.append((t, s))

def state_at(jumps, tq):
    lo, hi = 0, len(jumps) - 1
    while lo < hi:
        mid = (lo + hi + 1) // 2
        if jumps[mid][0] <= tq: lo = mid
        else: hi = mid - 1
    return jumps[lo][1]

grid = [state_at(jumps, k * dt) for k in range(int(jumps[-1][0] / dt))]
counts = np.zeros((2, 2))
for k in range(len(grid) - 1): counts[grid[k], grid[k + 1]] += 1
emp = counts / counts.sum(axis=1, keepdims=True)
print("empirical 1-step:", np.round(emp[0], 3))    # [0.819 0.181]
print("exp(Q*dt) row 0 :", np.round(expm(Q * dt)[0], 3))  # [0.819 0.181]

The empirically measured one-step matrix from the sampled path matches exp(Q·dt) = [0.819, 0.181] — sampling a CTMC at fixed intervals really is a discrete Markov chain, with the transition matrix being the matrix exponential of Chapter 7. The two worlds are one shift apart.

One more domain: reliability block diagrams

A concrete payoff outside ML: a redundant system with two servers, each failing at rate λ and repaired at rate μ, is a CTMC on the states {both up, one down, both down}. The generator Q books all the fail/repair rates; the stationary π gives steady-state availability; the time-to-absorption (Chapter 10) gives mean-time-to-system-failure if “both down” is catastrophic. Reliability engineering is, almost entirely, the CTMC machinery of this lesson applied to up/down state spaces — the same Q, π, and absorption-time tools, no new theory.

The same is true across the sciences: chemical kinetics is a CTMC on molecule counts (Gillespie was invented for exactly this); epidemic SIR-type models are CTMCs on compartment populations; population genetics tracks allele counts as birth-death chains. In every case the workflow is identical — write the generator Q, then ask exp(Qt) for transients, πQ=0 for equilibrium, or Gillespie for sample paths. You did not learn a niche tool; you learned the common language of stochastic dynamics in continuous time.

That is the note to end on. The discrete tick we threw away in Chapter 0 was never the point — rates were. Hold the generator-and-clocks picture and the whole zoo above reads as one idea with different state spaces bolted on. What you can simulate from Q, you understand.

If you want to go further, the natural next steps each add exactly one of the ingredients in the table: the HMM lesson adds emissions, the MDP lesson adds actions and rewards, and the Kalman and Bayes filter lessons add continuous state and noisy measurement. Every one of them reuses the rates, holding times, and balance you built here — you are no longer starting from zero.

python
# The redundant-system reliability CTMC: 2 servers, fail rate l, repair rate m.
import numpy as np
from scipy.linalg import expm

l, m = 0.1, 1.0   # per-server fail rate, repair rate
# states: 0 = both up, 1 = one down, 2 = both down
Q = np.array([[-2*l, 2*l, 0],
              [ m, -(m+l), l],
              [ 0, 2*m, -2*m]])
pi = expm(Q * 10000)[0]          # long-run distribution
print("availability (not both-down):", round(1 - pi[2], 5))
# ~0.99 -- redundancy turns two 91%-available servers into a 99% system.

The redundant-system generator gives a system availability near 99% from two individually-modest servers — redundancy quantified, straight from a 3-state CTMC. This is the reliability payoff in code: the same Q, π, and exp(Qt) tools, applied to an up/down state space.