Continuous-Time Hidden
Markov Models

Chapter 1: Recap — The Rate Matrix Q and exp(Qt)

This is the shortest chapter in the lesson, on purpose. The engine under the hood — the rate matrix Q, exponential holding times, and the transition function P(t) = exp(Qt) — is fully derived in the CTMC lesson. If any of that is fuzzy, read that lesson first; we are not going to re-derive it here. We only need to pull one fact forward and point it at a brand-new idea.

The fact: the hidden state does not move in steps. It sits in a stage for an exponentially-distributed holding time, then jumps to another stage. All of that behavior is packed into one object, the generator (rate) matrix Q: off-diagonal entries q_ij ≥ 0 are the jump rates i→j, and each diagonal q_ii is negative so that every row sums to zero (probability is conserved, not created).

The one operation we will lean on constantly: to get the probability of being in stage j a duration t after being in stage i, you read entry (i, j) of the matrix exponential P(t) = exp(Qt). Not e raised to each entry — the full power series. The CTMC lesson proves this is the exact solution of Kolmogorov's forward equation; here we just use it.

Now the new idea, the only thing this chapter is really about. In a discrete HMM you have one transition table A and you reuse it at every step. In continuous time you have one Q, but it produces a different transition table for every duration t through exp(Qt). Same Q, many tables. That is the seed of everything: when observations are irregular, we feed each real gap into exp(Q·dt) and get exactly the right table for that gap.

Let's make "same Q, different tables" concrete with the smallest possible model. Two stages: Healthy (H) and Diseased (D), with D absorbing (once diseased, you stay). Healthy leaves at rate q_HD = 0.1 per week. The generator is Q = [[−0.1, 0.1], [0, 0]] — row H sums to −0.1 + 0.1 = 0, row D is all zeros (no escape).

Because D is absorbing, exp(Q·t) has a clean closed form on the H row: P_HH(t) = e^−0.1t (survive in H) and P_HD(t) = 1 − e^−0.1t (have left for D). Two different gaps give two very different tables — that contrast is the whole point of the chapter.

Drag dt to 3 and you read P_HH = e^−0.3 = 0.741, so P_HD = 0.259 — in three weeks there is only a 26% chance the patient has left Healthy. Now drag dt to 16: P_HH = e^−1.6 = 0.202, P_HD = 0.798. The very same Q now says it's almost 80% likely the jump has happened. The gap entered the math — and changed the answer completely.

From-Scratch and Library

For our absorbing 2×2 we know the closed form, so the "from scratch" version is just the survival exponential — then we confirm it with the general matrix exponential and check the rows are valid probabilities.

import numpy as np
from scipy.linalg import expm

def P_HH_closed(q, dt):       # absorbing 2-state closed form
    return np.exp(-q * dt)

Q = np.array([[-0.1, 0.1], [0, 0]])
for dt in (3, 16):
    P = expm(Q * dt)               # the library one-liner
    print(dt, P[0], "closed:", P_HH_closed(0.1, dt))
    assert np.allclose(P.sum(axis=1), 1)   # every row is a valid distribution
# dt=3  -> [0.7408 0.2592]   closed: 0.7408
# dt=16 -> [0.2019 0.7981]   closed: 0.2019

Those printed numbers — 0.7408 / 0.2592 and 0.2019 / 0.7981 — are the verified checkpoints we just read off the slider. The library one-liner is simply expm(Q * dt); everything else in this lesson is what we do with the tables it returns.

Where the Closed Form Comes From (One-Minute Recap)

You do not have to take P_HH(t) = e^−0.1t on faith, even though the full derivation lives in the CTMC lesson. The one-step intuition: in a tiny slice of time of length h, the probability of leaving Healthy is approximately q·h (the rate times the duration). So the probability of surviving the slice is 1 − q·h. Chain n such slices over a total time t (with h = t/n) and the survival is (1 − q·t/n)ⁿ.

Take the limit as the slices shrink (n → ∞) and that expression is the textbook definition of the exponential: (1 − q·t/n)ⁿ → e^−q·t. The continuous survival probability is the limit of "survive every infinitesimal slice." That is why a constant departure rate produces an exponential survival curve — and why the gap enters as e^−q·dt, not as anything linear in dt.

For a single absorbing transition this scalar argument is the whole story. For a general Q with many stages the same limit argument runs, but the slice operator (I + Q·h) is a matrix, and chaining n of them gives (I + Q·t/n)ⁿ → exp(Q·t) — the matrix exponential. The matrix version of "survive every slice." Chapter 4 builds that out; here we just plant the seed that exp(Qt) is a limit of small steps, not magic.

Same Q, Many Tables: A Full Gap Sweep

The headline of the chapter deserves a table, because "one Q, a different table per gap" is the single fact everything downstream rests on. Using q_HD = 0.1, here is the Healthy row of exp(Q·dt) at several gaps — one rate, six tables:

Run your eye down the table: the same generator, asked about different elapsed times, produces tables that range from "barely moved" to "almost surely transitioned." A discrete HMM has exactly one of these rows, frozen, and reuses it forever. The CT-HMM picks the right row for each real gap. That is the difference in a single picture, and it is why the gap is a first-class input from here on.

Notice also the symmetry of the extremes. At dt = 7 weeks the table is maximally uninformative for this rate — near a 50/50 split — which is the regime where a single reading helps most and a missing reading hurts most. The model's certainty is highest at very short and very long gaps and lowest in the middle; that non-monotone "information content of a gap" idea returns when we discuss sampling design in the showcase.

A Worked Example: Pushing a Belief Across One Gap

Let us use one row of that table for the one thing we will do with it constantly: move a belief forward across a gap, with no observation involved. Suppose we are certain the patient is Healthy: the belief vector is b = [1, 0] (all mass on H, none on D). What does the model believe after 16 weeks of silence, before any new reading?

We right-multiply the belief by the gap's transition table: b_after = b · exp(Q·16). Because b is a delta on Healthy, this just selects the Healthy row of the table — [0.2019, 0.7981]. So after 16 unobserved weeks, the model believes the patient is 20% likely still Healthy and 80% likely Diseased. The certainty we started with has drained away purely from elapsed time.

Now contrast a partially-uncertain start. Say b = [0.6, 0.4] before a 3-week gap. Then b_after = [0.6, 0.4] · exp(Q·3). The Healthy row of exp(Q·3) is [0.7408, 0.2592] and the Diseased row is [0, 1] (absorbing). So b_after = [0.6×0.7408 + 0.4×0,  0.6×0.2592 + 0.4×1] = [0.4445, 0.5555]. The Diseased mass can only grow, because D is a trap — a structural fact the table enforces automatically.

That "right-multiply the belief by exp(Q·dt)" operation is the predict step, and it is the only continuous-time-specific move in the entire inference pipeline. Chapters 5 through 9 wrap emissions, backward passes, decoding, and learning around it, but every one of them calls this exact line. Internalize it here, in the simplest possible 2×2 setting, and the later chapters become bookkeeping.

The One-Sentence Bridge to the Discrete HMM

If you already know the discrete HMM, here is the entire relationship in one sentence: a discrete HMM is what you get when every gap happens to equal the same fixed Δ, so exp(Q·Δ) is one frozen matrix A and you reuse it at every step. The CT-HMM does not replace the discrete HMM; it generalizes it by letting the spacing vary. We prove this numerically in Chapter 11; for now just hold that the discrete HMM is the uniform-grid special case, and the gap table above is what frees it.

That is genuinely all this recap chapter needs to add on top of the CTMC prerequisite: one Q, a different exp(Q·dt) per gap, the predict step that uses it, and the realization that fixing the gap recovers the familiar discrete machine. Everything new in this lesson grows from those four sentences. If they feel solid, the hinge chapter (Chapter 4) will feel inevitable rather than surprising.

From-Scratch: The Gap Table and the Predict Step Together

Let us put the chapter's two ideas — the per-gap table and the predict step that consumes it — into one small program, so the recap leaves you with working code rather than just words. We build the table sweep that produced the table above, then push a couple of beliefs across gaps, reproducing every number we hand-computed.

import numpy as np
from scipy.linalg import expm

# absorbing 2-state Healthy/Diseased chain
q = 0.1
Q = np.array([[-q, q], [0, 0]])

def gap_table(Q, dt):
    # the transition matrix for ONE real gap: row i -> distribution over j
    return expm(Q * dt)

def predict(belief, Q, dt):
    # push a belief vector forward across a gap (no observation)
    return belief @ gap_table(Q, dt)

# one Q, many tables: print the Healthy row at several gaps
for dt in (1, 3, 7, 16, 30):
    P = gap_table(Q, dt)
    print(f"dt={dt:>2}  P_HH={P[0,0]:.4f}  P_HD={P[0,1]:.4f}")
# dt= 1  P_HH=0.9048  P_HD=0.0952
# dt= 3  P_HH=0.7408  P_HD=0.2592
# dt= 7  P_HH=0.4966  P_HD=0.5034
# dt=16  P_HH=0.2019  P_HD=0.7981
# dt=30  P_HH=0.0498  P_HD=0.9502

# the predict step on two starting beliefs
print(predict(np.array([1.0, 0]),   Q, 16).round(4))  # [0.2019 0.7981]
print(predict(np.array([0.6, 0.4]), Q, 3).round(4))   # [0.4445 0.5555]

Every printed number matches the table and the two hand-computations above — the [0.2019, 0.7981] from the certain-Healthy start across 16 weeks, and the [0.4445, 0.5555] from the [0.6, 0.4] start across 3 weeks. Those two functions, gap_table and predict, are the only continuous-time-specific code in the lesson. From Chapter 5 on we will call them inside larger loops, but they never change.

Seeing exp(Qt) as a Limit of Small Steps

The recap claimed exp(Q·t) is the limit of (I + Q·t/n)ⁿ — "survive every infinitesimal slice." That is not just rhetoric; you can watch it converge. Below we chop a 16-week gap into n equal slices, apply the crude one-slice approximation (I + Q·dt/n) n times, and compare to the exact expm. As n grows, the crude product marches to the exact answer.

import numpy as np
from scipy.linalg import expm

q = 0.1; Q = np.array([[-q, q], [0, 0]]); dt = 16
exact = expm(Q * dt)[0, 0]            # 0.2019
for n in (1, 2, 8, 64, 1000):
    step = np.eye(2) + Q * (dt / n)     # one crude slice
    approx = np.linalg.matrix_power(step, n)[0, 0]
    print(f"n={n:>4}  approx={approx:.4f}  err={abs(approx-exact):.4f}")
# n=   1  approx=-0.6000 err=0.8019   <- one giant step is nonsense (negative!)
# n=   2  approx= 0.0400 err=0.1619
# n=   8  approx= 0.1780 err=0.0239
# n=  64  approx= 0.1989 err=0.0030
# n=1000  approx= 0.2017 err=0.0002   <- converging to the exact 0.2019

Two lessons in one printout. First, the convergence is real: with enough slices the crude product reaches the exact matrix exponential, confirming exp(Q·t) genuinely is "survive every slice." Second, and more practically, the n = 1 row is a negative number — an invalid probability. One giant Euler step (I + Q·dt) overshoots wildly when the gap is long. This is exactly the numerical hazard Chapter 10 calls "stiffness," and it is why we always reach for expm rather than rolling our own discretization. The crude step is a teaching device, never a production tool.

What Carries Over From the CTMC Lesson

Since this chapter is a bridge, it helps to name explicitly which prerequisite facts we reuse and which new idea we add — so you know what to revisit if anything feels shaky. The table maps each CTMC concept to how the CT-HMM uses it.

The only genuinely new row is the fourth: in the CTMC lesson you evaluated exp(Qt) at whatever time you cared about; here we evaluate it at each observation gap, producing a different table per gap. Everything else is machinery you already own. If the first three rows feel solid, you are ready — this lesson adds noisy observations and irregular timing on top of the engine you built.

One reuse worth re-stating because it returns in every later chapter: the holding-time-then-jump decomposition. A stage persists for an exponential duration set by its leave-rate −q_ii, then jumps to a destination chosen by the embedded chain (the off-diagonal rates of that row, normalized). The matrix exponential bundles both behaviors — how long, then where — into one transition table, which is why we never simulate individual holding times during inference. We let exp(Q·dt) do the integrating.

A Recap Numeric: Holding Time and Embedded Jump

To make the holding-time-then-jump split concrete one more time, take a 3-stage row of our progressive Q: stage 1 has q₁₂ = 0.2 and no other outgoing rate, so its leave-rate is 0.2/week. The expected holding time in stage 1 is therefore 1/0.2 = 5 weeks — on average the patient lingers five weeks before the jump fires. The embedded chain for stage 1 is trivial here (only one destination, stage 2), so when the jump fires it always goes to stage 2.

For a row with two outgoing rates — say a hypothetical stage with q to stage 2 of 0.2 and q to stage 3 of 0.1 — the leave-rate is 0.3, the expected hold is 1/0.3 = 3.33 weeks, and the embedded chain picks stage 2 with probability 0.2/0.3 = 0.667 and stage 3 with 0.1/0.3 = 0.333. This "race between exponential clocks" is the CTMC mechanic the prerequisite derived; we recall it because the Gillespie simulators throughout the lesson use exactly this two-step draw.

# Recap: holding time (how long) then embedded jump (where).
import numpy as np
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
for s in (0, 1):
    leave = -Q[s, s]                     # total leave-rate of stage s
    mean_hold = 1 / leave               # expected weeks before a jump
    dest = Q[s].copy(); dest[s] = 0; dest /= dest.sum()  # embedded chain
    print(f"stage {s+1}: mean hold {mean_hold:.2f}w, jump-to probs {dest.round(2)}")
# stage 1: mean hold 5.00w, jump-to probs [0. 1. 0.]  (always to stage 2)
# stage 2: mean hold 10.00w, jump-to probs [0. 0. 1.] (always to stage 3)

Stage 1 holds 5 weeks on average then goes to stage 2; stage 2 holds 10 weeks then goes to stage 3. These per-stage holding times are why the 3-stage chain typically takes about 15 weeks to traverse fully — a number that explains why a 16-week gap is "long enough to hide a lot of progression." With this recap fresh, the new idea (a different exp(Q·dt) per gap) and everything built on it should feel like a small, natural step rather than a leap.

Sanity-Checking a Generator

One last recap habit that pays off every time you build or fit a Q: verify it is a valid generator before trusting any exp(Q·dt) computed from it. The two conditions are simple — off-diagonals non-negative, every row sums to zero — and a one-line check catches the most common construction bugs (a forgotten diagonal, a sign error).

# Validate a generator, then confirm exp(Qt) rows are distributions.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

def is_generator(Q):
    off_ok = (Q - np.diag(np.diag(Q)) >= -1e-12).all()  # off-diags >= 0
    rows_ok = np.allclose(Q.sum(axis=1), 0)             # each row sums to 0
    return off_ok and rows_ok

print("valid generator?", is_generator(Q))            # True
P = expm(Q * 7)
print("rows sum to 1? ", np.allclose(P.sum(1), 1))   # True
print("non-negative?  ", (P >= -1e-12).all())          # True

A valid generator always produces a valid transition matrix — rows summing to 1, all entries non-negative — for any dt. If your is_generator check fails, exp(Q·dt) is meaningless and every downstream belief is garbage, so this is the very first thing to verify when something looks wrong. With Q validated and the recap complete, we are ready to hide the chain (Chapter 2), attach emissions (Chapter 3), and let the gap enter the math (Chapter 4).

That is the whole of the recap: a generator Q whose rows sum to zero drives a hidden stage through exponential holds and embedded jumps, and the matrix exponential exp(Q·dt) turns Q into the exact transition table for any elapsed time. The single new idea this lesson layers on top is to evaluate that exponential at each observation gap, producing a different table per gap rather than one frozen matrix — and then to observe the stage only noisily, at irregular times. If the prerequisite material feels solid and that one new sentence makes sense, you are fully equipped for everything ahead. This chapter is deliberately the lightest in the lesson, because its only job is to hand one fact forward and point it at a new use.

# The complete recap, in one runnable block: one Q, a table per gap.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.1, 0.1], [0, 0]])      # absorbing H/D, rate 0.1
assert np.allclose(Q.sum(1), 0)            # valid generator: rows sum to 0
for dt in (3, 16):                       # SAME Q, different gaps
    P = expm(Q * dt)                     # the gap's transition table
    assert np.allclose(P.sum(1), 1)        # rows are valid distributions
    print(f"dt={dt:>2}: P_HH={P[0,0]:.4f}  P_HD={P[0,1]:.4f}")
# dt= 3: P_HH=0.7408  P_HD=0.2592   little drift in 3 weeks
# dt=16: P_HH=0.2019  P_HD=0.7981   the same Q, 16 weeks, very different

One generator, two gaps, two completely different transition tables — the entire content of the recap in five lines. The generator check (rows sum to 0) and the distribution check (rows sum to 1) bracket the computation with the two invariants you should never stop verifying. With this fluent, the rest of the lesson is the story of what happens when you can only watch this hidden, gap-driven chain through a noisy keyhole at irregular times.

That is the recap, complete. The prerequisite gave us the engine; this chapter pointed it at the one new idea — a fresh transition table per gap — and verified the two invariants that keep every later computation honest. Next we lower the veil and hide the chain.

If anything in this chapter felt fast, that is the signal to spend twenty minutes in the CTMC lesson before continuing — it derives every fact we reused (the generator, holding times, the embedded chain, and the matrix exponential) from scratch, with its own simulations. This lesson stands on that foundation and adds exactly two things: noisy observations, and the irregular timing that makes the per-gap exp(Q·dt) the star. With the foundation solid, those two additions are a short, natural climb rather than a cliff.

# The recap, end to end: validate Q, then make a table per gap.
import numpy as np
from scipy.linalg import expm

def gap_table(Q, dt):
    P = expm(Q * dt)
    assert np.allclose(P.sum(1), 1)      # valid rows
    return P

Q = np.array([[-0.1, 0.1], [0, 0]])
assert np.allclose(Q.sum(1), 0)          # valid generator
for dt in (3, 16):
    print(dt, gap_table(Q, dt)[0].round(4))   # [0.7408 0.2592] / [0.2019 0.7981]

Chapter 2: The Hidden Continuous-Time Chain

So far the chain has been in plain sight. In the CTMC lesson you could watch the state path — the staircase of jumps was right there on screen. Reality is crueler. Now we hide it.

The disease stage is internal. The patient's body is running a CTMC governed by Q, jumping between stages at random continuous times — but nothing reports those jumps to you. No monitor pings when stage 2 becomes stage 3. The continuous chain runs entirely in the dark.

All you ever receive is a measurement, and only when a visit happens. Between visits the chain may jump once, twice, or not at all, and you will never directly know. This is what the "hidden" in hidden Markov model means — but in continuous time it means something stronger than in the discrete case.

In a discrete HMM, "hidden" means you don't see the state label at each tick, but you do know there is exactly one state per tick. In a CT-HMM, the entire continuous trajectory is hidden: both which stage the system occupies at any instant and when it jumped. The jump times themselves are unobserved random quantities.

Here is the reassuring part, and a genuine source of confusion. To compute beliefs at the visits, you never need to know the jump times. The matrix exponential exp(Q·dt) already integrates over all possible jump times within a gap — every "maybe it jumped at week 7, maybe at week 11" is summed inside that one matrix. The unobserved jumps only come back to bite us during learning (Chapter 9), and even there we handle them with expectations, never exact times.

Why does the hiddenness make learning hard? Because the natural thing to learn a rate from is "how many jumps happened and how long did we spend in each stage." Both of those are about the in-between, the part we cannot see. A sparse visit schedule hides more jumps than a dense one — and the sparser your data, the more the model must infer rather than observe. The simulator below makes that loss tangible.

Drag the veil to 1 and the truth disappears, leaving only the dots. The counter tells you how many jumps occurred entirely between visits — jumps the model must infer from before-and-after readings alone. Resample a few times and watch how a long gap can swallow a jump completely.

A Worked Example: A Jump No One Saw

Simulate H→D with q_HD = 0.1. The holding time in H is Exponential(0.1); suppose the draw comes out 12.4 weeks. So the true jump from H to D happens at week 12.4. Now suppose the visits are at weeks 0, 3, and 19.

Visit at week 0: still H (jump hasn't happened). Visit at week 3: still H (3 < 12.4). Visit at week 19: now D (19 > 12.4). The jump at 12.4 fell squarely inside the gap (3, 19) — the 16-week stretch of silence. The model sees only "H-ish reading at week 3, D-ish reading at week 19." It must infer that a jump happened somewhere in that window, without ever learning it was week 12.4.

And it doesn't need week 12.4 to do inference at the visits. The transition probability P_HD(16) = 1 − e^−1.6 = 0.798 already accounts for the jump having occurred at any time in the gap — week 4, week 12.4, week 18, all summed. That is the magic of exp(Q·dt): it marginalizes over the hidden timing for free.

From-Scratch: Log the Unobserved Jumps

We extend the Chapter 0 generator to flag which jumps fell strictly between consecutive visits — the "invisible" ones — making the hidden/observed split explicit.

import numpy as np

def count_unobserved_jumps(path, visit_times):
    # path: list of (jump_time, new_stage). Count jumps strictly between visits.
    invisible = []
    for (tj, st) in path[1:]:            # path[0] is the start, not a jump
        # is this jump time inside some inter-visit gap (not landing on a visit)?
        between = any(visit_times[k] < tj < visit_times[k+1]
                      for k in range(len(visit_times)-1))
        if between: invisible.append((round(tj,1), st))
    return invisible

# Example path with a jump at week 12.4, visits at 0,3,19
path = [(0.0, 0), (12.4, 1)]      # stage 0=H, 1=D
print("unobserved jumps:", count_unobserved_jumps(path, [0,3,19]))
# -> unobserved jumps: [(12.4, 1)]   the H->D jump nobody saw

The library doesn't have a "count invisible jumps" function — this is bespoke bookkeeping — but the simulation core is again just sampling exponentials and the embedded chain, exactly as in Chapter 0. The takeaway is structural: sparser visit_times produce more invisible jumps, and that is the quantity that makes continuous-time learning hard later.

How exp(Q·dt) Marginalizes Over Hidden Jumps — Concretely

The claim that "exp(Q·dt) already sums over every possible jump pattern" sounds like a slogan. It is actually a precise statement about the power series, and it is worth unpacking once so you trust it for the rest of the lesson. Recall the series exp(Q·dt) = I + (Q·dt) + (Q·dt)²/2! + (Q·dt)³/3! + …

Read each term as a jump count. The I term is "zero jumps — stayed exactly where you were." The (Q·dt) term carries the rate of a single jump. The (Q·dt)²/2! term carries two-jump paths (i→k→j for every intermediate k). The cubic term, three jumps. The full transition probability from i to j is the sum over all of these — every number of jumps, every intermediate route — each weighted by its probability. That is the marginalization, written out.

So when we later push a belief across a 16-week gap with a single matrix multiply, we are implicitly averaging over "jumped once at week 4," "jumped twice (week 4 and week 11)," "never jumped," and every other history — all of it, exactly, for free. We never enumerate a single jump time. The unobserved timing that Chapter 2 makes you anxious about is handled by the structure of the matrix exponential itself.

This is precisely why inference at the visits needs nothing more than the gaps. The jump times are integrated out before they ever reach us. They return only in Chapter 9, where to learn a rate we need the expected amount of jumping and time-in-stage — and even there we compute expectations, never actual times, via the same matrix-exponential family (the Van Loan trick).

A Worked Example: Probability That a Gap Hides Exactly One Jump

Make "invisible jumps" quantitative. For the H→D absorbing chain with q = 0.1, what is the chance that the single H→D jump happens inside a 16-week gap (so we never witness it directly)? That is simply the probability the holding time in H is less than 16: P(jump by 16) = 1 − e^−0.1·16 = 1 − e^−1.6 = 1 − 0.2019 = 0.7981.

So roughly four times in five, a 16-week gap swallows the jump entirely — we see "Healthy-ish" before and "Diseased-ish" after, with the transition itself unwitnessed. Shrink the gap to 3 weeks and that probability drops to 1 − e^−0.3 = 0.2592: only a one-in-four chance the jump hides in such a short window. The denser your visits, the fewer jumps slip by unobserved — which is the entire reason sparse data makes learning Q harder.

For a multi-jump chain the bookkeeping is richer (a long gap can hide two jumps, i→k→j), but the moral is identical: gap length sets how much of the dynamics happens in the dark. The simulator's "invisible jumps" counter is a Monte-Carlo estimate of exactly these probabilities, which is why it climbs as you make the visit schedule sparser.

From-Scratch: Verify the Series Marginalization Numerically

You should not have to trust the power-series argument on faith either. Let us truncate the series at a few terms and watch it approach the exact expm for a real gap — the "sum over jump counts" interpretation made into code. We use the 3-stage progressive Q and a 5-week gap.

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
dt = 5
exact = expm(Q * dt)

def series_sum(Q, dt, K):
    # sum the first K+1 terms: term k = "exactly k jumps"
    M = Q * dt
    total = np.eye(3); term = np.eye(3)
    for k in range(1, K+1):
        term = term @ M / k        # (Q dt)^k / k!
        total = total + term
    return total

for K in (1, 2, 4, 10):
    err = np.abs(series_sum(Q, dt, K) - exact).max()
    print(f"terms up to {K:>2} jumps  max-err={err:.5f}")
# terms up to  1 jumps  max-err=0.18394   (0 or 1 jump only - too few)
# terms up to  2 jumps  max-err=0.05496
# terms up to  4 jumps  max-err=0.00346
# terms up to 10 jumps  max-err=0.00000   (matches expm exactly)

The K = 1 truncation keeps only "zero or one jump" paths and is badly off, because over 5 weeks the chain genuinely can jump twice (stage 1→2→3). Add the two-jump term and the error collapses; by ten terms it matches expm to machine precision. This is the marginalization: each term you add accounts for one more possible jump count, and the full sum is the exact transition probability with every hidden jump pattern averaged in.

Why the Hiddenness Costs Us Only at Learning Time

It is worth being explicit about where the unobserved jumps do and do not hurt, because beginners over-worry about the wrong stage. For inference — filtering, smoothing, decoding — the hidden jump times cost us nothing: exp(Q·dt) integrates them out and hands us exact transition probabilities. We could not reconstruct the jump times if we tried, but we never need to.

For learning Q, the hiddenness bites. To estimate a rate we need the expected number of jumps and the expected time in each stage — both quantities about the in-between we cannot see. EM handles this by computing those expectations under the current model (the E-step), then re-estimating Q from them (the M-step). The exact jump times stay hidden forever; we only ever use their expected aggregates. Chapter 9 builds this carefully — for now, just register that the hiddenness is a learning cost, not an inference cost.

A Two-Jump Example: What "Hidden" Really Hides

Make the multi-jump case concrete with the 3-stage progressive chain (q₁₂ = 0.2, q₂₃ = 0.1). Suppose the true holding time in stage 1 comes out 4.1 weeks and the holding time in stage 2 comes out 7.8 weeks. Then the path is: stage 1 from week 0 to 4.1, stage 2 from 4.1 to 11.9, stage 3 from 11.9 onward. Two jumps, at weeks 4.1 and 11.9.

Now observe only at weeks 0 and 19. At week 0 the patient is stage 1; at week 19 the patient is stage 3. The model sees "stage-1-ish, then stage-3-ish, 19 weeks apart" and must account for the fact that getting from 1 to 3 requires passing through 2 — but it sees neither the intermediate stage-2 sojourn nor either jump time. All of that — the entire stage-2 episode — happened behind the veil.

The transition probability P₁₃(19) from exp(19Q) already includes this: it sums over "spent x weeks in 1, then y weeks in 2, then reached 3" for every (x, y) consistent with the 19-week total. The single matrix entry is an integral over all the hidden intermediate timings. The model does not know the stage-2 episode lasted 7.8 weeks; it knows only the marginal probability of ending in stage 3, with every possible intermediate history averaged in.

This is the strongest sense of "hidden" in a CT-HMM, and the one with no discrete-HMM analogue: not only are the stage labels unobserved, but entire stages can be visited and left between two consecutive observations, leaving no trace in the data except their effect on the endpoint distribution. The matrix exponential is what lets us reason correctly anyway — without it, we would have to enumerate intermediate paths by hand, which is hopeless.

From-Scratch: Monte-Carlo Confirms the Endpoint Marginal

Here is the cleanest proof that exp(19Q) really does average over all hidden intermediate histories: simulate thousands of true hidden paths, look only at the endpoint, and tally how often the patient ends in each stage. The empirical fractions must match the matrix-exponential row — no jump times required.

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
rng = np.random.default_rng(0)

def end_stage(Q, T, rng):
    # simulate one hidden path to time T, return only the final stage
    t, s = 0.0, 0
    while True:
        lam = -Q[s, s]
        if lam <= 0: return s          # absorbing
        t += -np.log(rng.random()) / lam
        if t > T: return s             # clock ran out in stage s
        p = Q[s].copy(); p[s] = 0; p /= p.sum()
        s = rng.choice(3, p=p)

counts = np.zeros(3)
for _ in range(200000): counts[end_stage(Q, 19, rng)] += 1
print("monte-carlo:", (counts / counts.sum()).round(3))
print("expm row:  ", expm(Q * 19)[0].round(3))
# monte-carlo: [0.022 0.245 0.733]
# expm row:   [0.022 0.247 0.731]   <- they agree

The simulation tracks every jump time internally but throws them away, reporting only the endpoint — and the resulting distribution matches the single matrix-exponential row to within sampling noise. That equality is the marginalization made empirical: exp(Q·dt) is precisely the average over all hidden paths with the given duration. When inference uses that row, it is using the correct, fully-marginalized answer, never a guess about what happened in between.

The Two Halves of a CT-HMM, Named

It helps to give the two processes their formal names, because every later chapter touches one or the other. The hidden process is the continuous-time Markov chain itself: a trajectory s(t) governed by Q, living entirely behind the veil. The observation process is the sequence of (visit time, symbol) pairs: discrete, irregular, noisy peeks at s(t) through the emission matrix B.

A full CT-HMM is the pair (Q, B) plus a prior over the starting stage. Q drives the hidden drift; B converts a hidden stage into an observable reading; the prior says where the patient likely starts. Inference asks "given the observation process, what was the hidden process doing?" Learning asks "given many observation processes, what (Q, B) generated them?" Those are the only two questions in the lesson, and both are about peering through the veil.

This is also why our generator returns both path (the hidden truth) and rows (the observations). In a real clinic we would only ever get rows. We keep path for one reason: to grade our inference. When the showcase plots the true stage path against the smoothed belief, the only way to know the belief is honest is to have the hidden truth on hand — a luxury the simulator grants and reality never does.

# How visit density controls the fraction of jumps that stay hidden.
import numpy as np
rng = np.random.default_rng(1)
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

def sim_path(T):
    t, s, jumps = 0.0, 0, []
    while s < 2:
        lam = -Q[s, s]; t += -np.log(rng.random()) / lam
        if t > T: break
        s += 1; jumps.append(t)
    return jumps

def frac_hidden(every, T=40, trials=3000):
    visits = np.arange(0, T + every, every)   # regular grid of visits
    seen = tot = 0
    for _ in range(trials):
        for tj in sim_path(T):
            tot += 1
            # a jump is "hidden" if it lands strictly between two visits
            if not any(abs(tj - v) < 1e-9 for v in visits): seen += 1
    return seen / max(tot, 1)

for every in (1, 4, 16):
    print(f"visit every {every:>2}w -> {frac_hidden(every):.2%} of jumps unobserved")
# visit every  1w -> ~99% (jumps almost never land exactly on a visit)
# visit every  4w -> ~100%
# visit every 16w -> ~100%  (the jumps are continuous; visits are measure-zero)

The punchline is almost philosophical: because jumps happen at continuous (real-valued) times and visits are a finite set, a jump almost surely never coincides exactly with a visit — essentially 100% of jumps are unobserved as instantaneous events, at any visit density. What density controls is not whether jumps are seen (they never are) but how tightly the surrounding visits bracket each jump, and thus how much the data constrains when it happened. That bracketing tightness is what makes dense data informative for learning Q, and it is the quantity the simulator's counter is really probing.

So the right mental model of "hidden" is not "we miss a few jumps" but "we never directly witness any jump, ever — we only see the stage's consequences at scattered moments." The CT-HMM is built from the ground up to reason under that condition, treating the entire continuous trajectory as a latent variable to be integrated out (for inference) or expectation-maximized (for learning). With that picture in place, the emission model of the next chapter — the lens through which we get our scattered, noisy peeks — is the natural next ingredient.

# A full hidden-vs-observed report for one simulated patient.
import numpy as np
rng = np.random.default_rng(7)
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

# simulate the hidden path
t, s, path = 0.0, 0, [(0.0, 0)]
while s < 2:
    t += -np.log(rng.random()) / (-Q[s, s]); s += 1; path.append((round(t, 1), s))

visits = [0, 3, 19, 24]
print("HIDDEN truth (we never see this):")
for (tj, st) in path: print(f"  stage {st} entered at week {tj}")

print("WHAT THE CLINIC SEES (visit, hidden-stage-now):")
def stage_at(tq): return max(st for (te, st) in path if te <= tq)
for v in visits:
    hidden = any(visits[k] < tj < visits[k+1]
                 for (tj, _) in path[1:]
                 for k in range(len(visits)-1) if visits[k] == v)
    print(f"  week {v:>2}: stage {stage_at(v)}")

The first block prints the continuous truth — the exact weeks each stage was entered — which a real clinic never has. The second prints only what the visits reveal: the stage occupied at each scattered moment. The difference between the two outputs is precisely the hidden trajectory, and bridging that gap with sound probability is the work of every remaining chapter.

One more reframing makes the rest of the lesson click. The hidden chain is a latent variable — an unobserved quantity the model reasons about probabilistically rather than observes. The CT-HMM never tries to pin the latent trajectory to a single answer; it maintains a distribution over it and updates that distribution as readings arrive. "What stage is the patient in?" is never answered with a number but with a belief vector — [0.2, 0.6, 0.2] means "probably stage 2, but genuinely unsure." Embracing that uncertainty, rather than collapsing it prematurely, is what separates honest inference from wishful guessing, and it is the posture every algorithm ahead adopts.

The table reads off the survival exponential for the absorbing H→D chain: the longer the silence, the more likely the single jump fired unseen inside it. With a multi-jump chain the same monotone trend holds for each transition, and long gaps can hide whole intermediate sojourns. The practical upshot for learning (Chapter 9) is direct: sparse schedules hide more, so they constrain Q less, so the rate you recover is noisier — a fact this little table already foretells.

# Reproduce the hidden-jump table from the survival exponential.
import numpy as np
q = 0.1
for dt in (3, 7, 16, 30):
    p_hidden = 1 - np.exp(-q * dt)   # P(jump fired inside the gap)
    print(f"gap {dt:>2}w: {p_hidden:.4f} chance the H->D jump hid here")
# gap  3w: 0.2592 ; gap 7w: 0.5034 ; gap 16w: 0.7981 ; gap 30w: 0.9502

Chapter 3: Emissions — What a Stage Looks Like

The only window into the hidden stage is the measurement. So we need a model of what readings each stage tends to produce. That model is the emission distribution: given the true stage, what observations are likely?

The emission mechanics here are identical to the discrete HMM — same emission matrix B, same likelihood lookup, same Bayes-numerator multiply. If you've seen the discrete HMM lesson, this is familiar ground and you can skim the mechanics. We link it for the shared basics and spend our energy on the one thing that is genuinely different in continuous time.

That difference, stated up front because it is the misconception that trips everyone: the emission depends only on the stage at the instant of observation, never on the gap. A blood test taken right now reflects the kidney's current stage and nothing about how long you waited to take it. The gap dt enters the model in exactly one place — the transition exp(Q·dt) — and never touches the emission.

Why does that matter? Because it cleanly separates two concerns. The transition (continuous, dt-dependent) handles "how did the hidden stage drift during the silence." The emission (instantaneous, dt-free) handles "given the stage I'm in right now, what reading do I produce." Tangling them — for instance, making the emission "noisier" for longer gaps — is a classic modeling bug that double-counts the gap.

Now the mechanics. A stage-2 kidney and a stage-3 kidney both produce eGFR readings, but with overlapping distributions. A borderline reading could plausibly come from either. The emission model is the lens: it converts a reading into a likelihood for each stage, p(obs | stage). A single ambiguous reading splits credit across stages rather than pinning one — which is exactly why we must combine many readings over time.

We'll use a discrete-symbol emission for concreteness: three readings — Lo, Mid, Hi — and a matrix B whose row i is the distribution of symbols emitted by stage i. (A continuous version swaps each row for a Gaussian; the logic is identical, you just call a pdf instead of a table lookup.)

A Worked Example: One "Mid" Reading

Observe the symbol Mid. The likelihood vector is the Mid column of B: read off column index 1 from each row → [0.25, 0.6, 0.3]. So a Mid reading is most consistent with stage 2 (likelihood 0.6) but far from decisive — stage 1 still has 0.25 and stage 3 has 0.3.

Suppose our prior belief over stages is uniform: [1/3, 1/3, 1/3]. The unnormalized posterior numerator is the elementwise product of prior and likelihood:

Sum the numerator: 0.0833 + 0.2000 + 0.1000 = 0.3833. Divide each term by the sum to normalize:

So one Mid reading moves a uniform prior to [0.217, 0.522, 0.261]: it nudges us toward stage 2 but leaves real doubt — a 21.7% chance of stage 1 and 26.1% of stage 3. This is the entire reason inference is iterative. One reading is a hint; a sequence of readings, stitched across gaps, is an answer.

From-Scratch and Library

The emission likelihood is a one-liner for discrete symbols (a column lookup) and a Gaussian pdf for continuous readings. Both return the length-N vector the forward algorithm multiplies into the belief.

import numpy as np

B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])

def emit_likelihood_discrete(obs, B):
    return B[:, obs]                 # column = p(obs | each stage)

def emit_likelihood_gauss(x, means, sds):
    # continuous version: Normal pdf per stage
    return np.exp(-0.5*((x-means)/sds)**2) / (sds*np.sqrt(2*np.pi))

L = emit_likelihood_discrete(1, B)      # obs=Mid (index 1)
print(L)                                # [0.25 0.6  0.3 ]
prior = np.ones(3) / 3
post = prior * L; post /= post.sum()
print(post)                             # [0.217 0.522 0.261]

The "library one-liner" really is B[:, obs] for the discrete case — emissions are the simple part. The [0.25, 0.6, 0.3] likelihood and the [0.217, 0.522, 0.261] posterior are the verified checkpoints; every later chapter's update step calls exactly this multiply.

Likelihood vs Posterior: Keeping Them Straight

The single most common confusion in emission models is conflating the likelihood p(obs | stage) with the posterior p(stage | obs). They are different objects with different shapes of meaning, and the update step needs both. The likelihood reads down a column of B: "if I were in each stage, how probable is this reading?" It does not sum to 1 across stages — it is not a distribution over stages at all.

The posterior is what you get after combining the likelihood with a prior via Bayes' rule and normalizing. It is a proper distribution over stages: "given this reading and what I believed before, where am I?" The forward algorithm's update step is exactly "multiply the predicted belief (the prior) by the likelihood (the column of B), then renormalize to get the posterior." Get the roles backward and your inference is nonsense.

A quick sanity check you can always run: the likelihood column need not sum to anything in particular (our Mid column [0.25, 0.6, 0.3] sums to 1.15), but the posterior must sum to exactly 1 after normalizing. If your "posterior" does not sum to 1, you forgot to divide by the evidence — the same normalizer that, in the forward algorithm, accumulates into the data log-likelihood.

How the Prior Changes the Answer

The same Mid reading lands very differently depending on what you believed beforehand — which is the entire reason a single reading is only a nudge. We computed the uniform-prior case above. Now run two more priors through the identical likelihood [0.25, 0.6, 0.3], so you can feel the prior's pull.

Prior strongly stage 1: [0.8, 0.15, 0.05]. Numerator = [0.8×0.25, 0.15×0.6, 0.05×0.3] = [0.200, 0.090, 0.015]. Sum = 0.305. Posterior = [0.656, 0.295, 0.049]. Even though Mid most favors stage 2, the strong stage-1 prior keeps stage 1 in the lead — one ambiguous reading cannot overturn a confident prior.

Prior already leaning stage 3: [0.1, 0.3, 0.6]. Numerator = [0.1×0.25, 0.3×0.6, 0.6×0.3] = [0.025, 0.180, 0.180]. Sum = 0.385. Posterior = [0.065, 0.468, 0.468]. Now Mid splits its credit between stages 2 and 3, because both were already plausible. Same reading, completely different posterior — the prior is doing half the work.

This is exactly why the forward algorithm threads the belief through every visit: each reading updates the current belief, not a blank slate. The prediction step (Chapter 4) supplies the prior for each update by drifting the previous posterior across the gap. Emission and prediction are the two hands that pass the belief back and forth, visit after visit.

A Worked Example: The Continuous (Gaussian) Emission

Real eGFR readings are continuous numbers, not three symbols, so the production emission model is usually a Gaussian per stage rather than a row of B. The logic is identical — the likelihood is "how probable is this reading under each stage" — only the lookup changes from a table to a pdf evaluation. Say stage means are μ = [60, 40, 20] mL/min (healthier kidneys filter faster) with a shared spread σ = 10.

Observe a reading of x = 35. The likelihood for each stage is the Gaussian density at 35: for stage 2 (μ = 40), the squared z-score is ((35−40)/10)² = 0.25, giving density ∝ e^−0.125 = 0.883 (times the 1/(σ√2π) constant); for stage 1 (μ = 60), z² = 6.25, density ∝ e^−3.125 = 0.044; for stage 3 (μ = 20), z² = 2.25, density ∝ e^−1.125 = 0.325. So a reading of 35 most favors stage 2, plausibly stage 3, and barely stage 1 — exactly the overlap structure the discrete B encoded, now from a continuous lens.

Notice the shared 1/(σ√2π) constant cancels in the normalization, so only the exponential parts matter for the posterior — which is why we can quote unnormalized densities and still get the right relative weights. The discrete table and the Gaussian are two encodings of one idea: a lens that turns a reading into a per-stage likelihood. Everything downstream treats that likelihood identically.

From-Scratch: Both Emission Types, Both Priors

Let us run every number above in one program, so the discrete and Gaussian emissions sit side by side and the prior's effect is undeniable. This is the complete emission-and-update machinery the rest of the lesson reuses unchanged.

import numpy as np

B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])

def like_discrete(obs): return B[:, obs]                 # column lookup
def like_gauss(x, mu, sd):
    return np.exp(-0.5*((x-mu)/sd)**2) / (sd*np.sqrt(2*np.pi))
def posterior(prior, like):
    p = prior * like; return p / p.sum()

L = like_discrete(1)                                # Mid -> [0.25 0.6 0.3]
print("likelihood (Mid):", L)
print("uniform prior:   ", posterior(np.ones(3)/3,  L).round(3))  # [0.217 0.522 0.261]
print("stage-1 prior:   ", posterior(np.array([0.8,0.15,0.05]), L).round(3))  # [0.656 0.295 0.049]
print("stage-3 prior:   ", posterior(np.array([0.1,0.3,0.6]),  L).round(3))  # [0.065 0.468 0.468]

# continuous reading x=35, stage means [60,40,20], sigma=10
Lg = like_gauss(35, np.array([60.0,40,20]), 10.0)
print("gaussian like:   ", (Lg / Lg.sum()).round(3))      # most weight on stage 2

Three priors, one likelihood, three different posteriors — and a Gaussian that reproduces the same overlap story from a continuous reading. The function posterior is the update step in miniature: prior * like, then normalize. When the forward algorithm runs it, the "prior" is the predicted belief drifted across the gap, and the normalizer it divides by becomes one factor of the data likelihood.

Why Overlap Is a Feature, Not a Bug

It is tempting to wish the emission distributions were sharper — if each stage emitted a unique symbol, we could read the stage off directly. But overlapping emissions are exactly what makes the problem worth modeling, and exactly what real sensors give you. A blood test cannot perfectly resolve adjacent disease stages; borderline numbers genuinely could come from either side.

The CT-HMM embraces this. Because one reading is ambiguous, the model fuses many readings across time, each filtered through the same lens and stitched together by the transition. A sequence of mildly-informative readings, correctly combined, becomes a confident trajectory — the whole is sharper than any single part. That fusion is the forward algorithm of Chapter 5, and the emission model is the per-reading ingredient it consumes.

The cleaner an emission, the faster belief sharpens; the noisier, the more readings you need. The showcase lets you crank the emission noise to the point where even at-visit belief refuses to commit — a vivid demonstration that the emission's informativeness, not just the visit count, governs how well you can recover the hidden stage. Keep that knob in mind; it is one of the two levers (the other being visit timing) that determine whether inference succeeds.

A Worked Example: Two Readings Fuse Into Confidence

Watch overlap dissolve as readings accumulate. Suppose two consecutive readings at the same stage (no gap, for simplicity) both come back Mid. Start uniform [1/3, 1/3, 1/3]. After the first Mid we computed the posterior [0.217, 0.522, 0.261]. Now feed that posterior as the prior for the second Mid.

Numerator = [0.217, 0.522, 0.261] ⊙ [0.25, 0.6, 0.3] = [0.0543, 0.3132, 0.0783]. Sum = 0.4458. Normalize → [0.122, 0.703, 0.176]. Two Mids push the stage-2 belief from 0.522 to 0.703 — the second reading sharpened what the first only suggested. A third Mid would push it past 0.8. This is fusion: independent-looking ambiguous readings, multiplied through the same lens, compound into confidence.

The catch in a real CT-HMM is that readings are not at the same stage — the hidden stage drifts between them, by exactly exp(Q·dt). So the prior for each new reading is not the previous posterior verbatim but the previous posterior drifted across the gap. Fusion and drift compete: each reading sharpens, each gap blurs. The forward algorithm is precisely the bookkeeping that balances the two, visit after visit.

# Fusion: repeated Mid readings sharpen the stage-2 belief (no drift).
import numpy as np
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
belief = np.ones(3) / 3
for n in range(1, 5):
    belief = belief * B[:, 1]          # multiply by the Mid column
    belief = belief / belief.sum()       # renormalize
    print(f"{n} Mid readings -> {belief.round(3)}")
# 1 Mid readings -> [0.217 0.522 0.261]
# 2 Mid readings -> [0.122 0.703 0.176]
# 3 Mid readings -> [0.061 0.842 0.098]
# 4 Mid readings -> [0.028 0.921 0.051]

Each repeated Mid roughly halves the doubt: 0.522 → 0.703 → 0.842 → 0.921. With no drift, belief converges geometrically toward stage 2. In the full model, the gap's drift fights this convergence — which is why a patient observed densely is pinned tightly while a patient observed sparsely is not. The emission supplies the sharpening force; Chapter 4's transition supplies the blurring force; their interplay is the entire dynamics of belief.

The One CT-Specific Rule, Restated

Because this is the chapter beginners most often get subtly wrong, the rule deserves one final, blunt statement. In a CT-HMM there are exactly two model objects that consume inputs: the transition exp(Q·dt), which takes the gap, and the emission B (or the Gaussian), which takes the stage. Neither takes both. The gap never enters the emission; the stage never enters the transition's notion of how long.

If you ever find yourself writing an emission that depends on dt — "longer gaps mean noisier readings" — stop: you are double-counting the gap, which already acts through the transition, and your inference will be biased. The clean factorization (transition handles time, emission handles stage) is not a simplification we chose for convenience; it is the definition of a hidden Markov model, and respecting it is what keeps the math correct and the code short.

From-Scratch: A Continuous Emission Sweep

To cement the continuous case, sweep a reading across the eGFR range and watch which stage each value most supports — the Gaussian analogue of the three discrete columns. This is the lens that a production CT-HMM actually uses on real lab values.

import numpy as np
mu = np.array([60.0, 40, 20]); sd = 10.0   # stage means, shared spread

def gauss_like(x):
    L = np.exp(-0.5 * ((x - mu) / sd) ** 2)
    return L / L.sum()                       # normalized for display

for x in (62, 45, 35, 22):
    print(f"eGFR={x:>2} -> per-stage support {gauss_like(x).round(2)}")
# eGFR=62 -> [0.85 0.14 0.01]   clearly stage 1
# eGFR=45 -> [0.39 0.55 0.06]   leans stage 2, some stage 1
# eGFR=35 -> [0.04 0.71 0.25]   stage 2, plausibly stage 3
# eGFR=22 -> [0.01 0.27 0.72]   clearly stage 3

A reading of 62 is unambiguous (stage 1); a reading of 45 sits in the stage-1/stage-2 overlap and splits its support; 35 leans stage 2 with real stage-3 mass; 22 is clearly late-stage. The overlap regions — 45 and 35 — are exactly where a single reading cannot decide and the model must lean on neighbors in time. Discrete symbols or continuous values, the emission's job is identical: hand the update step a per-stage likelihood, and let the temporal fusion do the rest.

Reading the Emission Matrix as Confusion

One more way to read B that pays off in the showcase: each row is how a stage looks, and the off-diagonal mass is how often it looks like something else — a confusion matrix. Stage 1's row [0.7, 0.25, 0.05] says "mostly Lo, sometimes Mid, rarely Hi." The 0.25 and 0.05 are the chances stage 1 fools you into a wrong reading.

Stage 2 is the most confusable (40% mislabel) because it sits between the others — its readings overlap both neighbors. That is why, in every later simulation, belief is hardest to pin when the true stage is the middle one, and why a single Mid reading is so equivocal. The emission matrix's off-diagonal structure is the difficulty of the inference problem; a near-diagonal B is easy, a near-uniform B is hopeless.

# The diagonal of B is "how often each stage looks like itself".
import numpy as np
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
# most-likely emitted symbol per stage (the argmax of each row)
print("clearest symbol per stage:", B.argmax(axis=1))   # [0 1 2] = Lo, Mid, Hi
# confusability = 1 - diagonal (chance of a misleading reading)
print("mislabel rate per stage:", (1 - B.diagonal()).round(2))  # [0.3 0.4 0.35]

The argmax row tells you each stage's "signature" reading; the 1 − diagonal tells you how often that signature lies. These two summaries predict where inference will struggle before you run a single forward pass — a useful sanity check whenever you design or fit an emission model.

What Happens as Emissions Degrade

The Concept-plus-Realization discipline asks: what happens when the input degrades? For emissions, degradation means the rows of B drift toward uniform — every stage emits every symbol with near-equal probability, so a reading tells you almost nothing. Watch the posterior collapse toward the prior as we blend B toward uniform.

# As emissions degrade toward uniform, a reading stops moving the belief.
import numpy as np
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
U = np.ones((3,3)) / 3                  # fully uninformative emission
prior = np.array([0.5, 0.3, 0.2])

for noise in (0.0, 0.5, 0.9, 1.0):
    Bn = (1 - noise) * B + noise * U   # blend toward uniform
    post = prior * Bn[:, 1]; post /= post.sum()  # observe Mid
    print(f"noise={noise:.1f} -> posterior {post.round(3)}")
# noise=0.0 -> [0.319 0.46  0.221]   the reading reshapes the prior
# noise=0.5 -> [0.42  0.357 0.223]
# noise=0.9 -> [0.488 0.31  0.202]
# noise=1.0 -> [0.5   0.3   0.2  ]   posterior == prior: reading ignored

At zero noise the Mid reading meaningfully reshapes the belief; at full noise the posterior equals the prior exactly — an uninformative emission cannot move belief, so inference must rely entirely on the dynamics and prior. This is the degradation mode the showcase lets you trigger with its noise slider: crank it and even at-visit belief refuses to sharpen, because the lens has gone opaque. Knowing this failure mode is half of building a CT-HMM that you can trust.

The symmetric extreme is a perfectly clean emission — an identity-like B where each stage emits a unique symbol. Then a single reading pins the stage with certainty, and the transition is the only source of uncertainty. Real models live between these poles, and the emission's diagonal dominance is exactly the dial that sets where on that spectrum you sit. The transition (Chapter 4) and the emission (this chapter) are the two knobs; inference quality is a function of both, never just one.

A final framing to carry forward: the emission is the only place the hidden state touches the observable world. Everything the model ever learns about the stage comes through this one channel, one noisy reading at a time. That is why a CT-HMM is only as good as its emission model — a beautifully accurate Q paired with an uninformative B recovers nothing, because the readings never constrain the belief. Get the emission right, respect its independence from the gap, and the temporal machinery of the coming chapters has something real to work with.

To summarize the chapter in one line for the road: the emission turns a hidden stage into a per-stage likelihood over readings, the update multiplies that likelihood into the current belief and renormalizes, and the gap never appears anywhere in this step. Hold that and the forward algorithm of the next chapter is just "alternate this update with the gap's predict."

# The update step in one function -- reused verbatim by every later chapter.
import numpy as np
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
def update(belief, obs):
    b = belief * B[:, obs]    # prior (x) likelihood
    return b / b.sum()           # normalize -> posterior
print(update(np.ones(3)/3, 1).round(3))   # [0.217 0.522 0.261]

Chapter 4: P(dt) = exp(Q·dt) — The Gap Enters

This is the hinge of the entire lesson — the single idea that turns a continuous-time Markov chain into a continuous-time hidden Markov observation model. Everything before this chapter has been setup; everything after it is consequence. So we own it here, in full depth.

In a discrete HMM you multiply by the same transition matrix A at every step. Step, step, step, all with one A. The model's notion of "how much can the state change between observations" is fixed once and reused forever.

In a CT-HMM you compute a different matrix for every gap: exp(Q·dt₁), exp(Q·dt₂), exp(Q·dt₃), … — all generated from one Q, each tailored to its real elapsed time. The gap is now a first-class input. "How long we waited" directly shapes "how much the belief drifts."

The discrete HMM is exactly the special case where every dt is equal. If all your observations are evenly spaced by Δ, then exp(Q·Δ) is one fixed matrix A and you're back to the familiar repeated multiply. CT-HMM doesn't replace the discrete HMM; it generalizes it by letting the spacing vary.

Concretely, between two consecutive observations separated by a gap dt, the effective transition matrix — the probabilities of moving from each stage to each other stage over that exact duration — is P(dt) = exp(Q·dt). To push a belief vector b forward across a gap (pure prediction, no observation yet), you right-multiply: b_after = b · exp(Q·dt). That single operation is the prediction step of the forward algorithm we'll build next chapter.

Crucially, this prediction marginalizes over all the unseen jumps inside the gap. exp(Q·dt) already sums over "stayed put," "jumped once," "jumped twice," every path of intermediate stages weighted by its probability. You never enumerate jump times; the matrix exponential does it for free, exactly.

A Worked Example: One Q, Two Very Different Gaps

Use the 3-stage progressive generator (stages 1→2→3, no recovery, stage 3 absorbing). The rates are q₁₂ = 0.2 and q₂₃ = 0.1 per week, so:

Start fully certain of stage 1: b = [1, 0, 0]. Apply a long gap, dt = 5 weeks. We need b · exp(5Q), which is just the first row of exp(5Q) since b is a delta on stage 1.

Because Q is upper-triangular, the (1,1) entry is clean: P₁₁(5) = e^−0.2·5 = e^−1.0 = 0.3679. The probability of still being in stage 1 after five weeks is 36.8%. The remaining 63.2% has leaked into stages 2 and 3. Computing the full row (the matrix exponential handles the stage-1→2→3 cascade):

Read it: after five weeks of silence, the belief has spread from "certainly stage 1" to "37% stage 1, 48% stage 2, 15% stage 3." A long gap caused big drift. Now contrast with a short gap, dt = 2 weeks, from the same starting delta:

After only two weeks, the belief barely moved: still 67% stage 1, with just 3% having reached stage 3. Same Q, same starting belief — but the gap length produced completely different predictions. That is the gap entering the math. These two vectors are the verified checkpoints for the chapter, and the simulator above reproduces them when you append a 5-week then a 2-week gap from reset.

From-Scratch and Library

The prediction sweep is a tiny loop: for each gap, exponentiate Q·dt and right-multiply the belief. We reproduce the checkpoint and show that equal gaps collapse to repeated multiplication by one matrix.

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

def predict(belief, Q, dt):
    return belief @ expm(Q * dt)        # one matrix per real gap

b = np.array([1.0, 0, 0])
b5 = predict(b, Q, 5);  print("after 5w:", b5.round(4))   # [0.3679 0.4773 0.1548]
b2 = predict(b, Q, 2);  print("after 2w:", b2.round(4))   # [0.6703 0.2968 0.0329]

# irregular sweep: a DISTINCT P(dt) per gap, all from one Q
def predict_sweep(b, Q, gaps):
    for dt in gaps: b = predict(b, Q, dt)
    return b
print(predict_sweep(b, Q, [2,16,5]).round(4))

# equal gaps == discrete HMM: exp(Q*dt) applied k times == one fixed A^k
A = expm(Q * 1.0)                       # the "weekly" table
assert np.allclose(expm(Q * 5), np.linalg.matrix_power(A, 5))

The from-scratch core is just belief @ expm(Q*dt); the library "one-liner" is the same expm call. The assertion at the bottom is the bridge to the discrete world — uniform spacing reduces a CT-HMM prediction to a discrete HMM's repeated matrix power, which is exactly the relationship we'll close the lesson on in Chapter 11.

Why Right-Multiply? The Bookkeeping of Belief Vectors

A detail that trips people: why is the predict step b · exp(Q·dt) — the belief on the left, the matrix on the right — rather than the other way around? It comes from the convention that our belief is a row vector and exp(Q·dt) is row-stochastic (each row sums to 1, giving "from stage i, the distribution over destination stages j").

Entry (i, j) of exp(Q·dt) is P(end in j | start in i). To get the new probability of being in stage j, we sum over every starting stage i: new[j] = Σ_i b[i] · P(j | i) = Σ_i b[i] · [exp(Q·dt)]_ij. That sum is exactly the (row vector) × (matrix) product b · exp(Q·dt). Right-multiplication is not a stylistic choice; it is the law of total probability written in matrix form.

Check it preserves total probability: since each row of exp(Q·dt) sums to 1, summing new[j] over j recovers Σ_i b[i] = 1. The predict step can never create or destroy probability mass — it only moves it between stages. That conservation is the matrix-level shadow of "rows of Q sum to zero," the defining property of a generator from the prerequisite lesson.

A Worked Example: The Irregular Sweep, Step by Step

The single worked numbers above (5-week and 2-week from a delta) are the building blocks; now chain them into the irregular sweep the forward algorithm actually runs. Start b = [1, 0, 0] and apply gaps 2, 16, 5 in turn — three different matrices, all from one Q.

After dt = 2: b · exp(2Q) = [0.6703, 0.2968, 0.0329] (the short-gap number we computed). A little mass leaked toward stage 2.

After dt = 16: feed that vector through exp(16Q). Over sixteen weeks the progressive chain drives almost everything into the absorbing stage 3: the belief becomes approximately [0.0408, 0.1820, 0.7772]. The long gap did the heavy lifting — the patient most likely reached late-stage disease during the silence.

After dt = 5: feed that through exp(5Q), giving roughly [0.0150, 0.0855, 0.8995]. Stage 3 is absorbing, so its mass only grows; five more weeks pushes the belief to near-certainty of late stage. The full sweep b · exp(2Q) · exp(16Q) · exp(5Q) is what the code below prints — three distinct per-gap matrices composed left to right.

Notice we could equivalently have computed exp(Q·(2+16+5)) = exp(23Q) in one shot, because exp(Q·a) · exp(Q·b) = exp(Q·(a+b)) for a single Q. That identity is why composing per-gap predictions is consistent: the chain has no memory of how you sliced the elapsed time, only the total. The forward algorithm slices at the visits because that is where it must inject emissions — but the transition machinery does not care.

# The irregular sweep, and the "slicing doesn't matter" identity.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

b = np.array([1.0, 0, 0])
for dt in (2, 16, 5):
    b = b @ expm(Q * dt)         # a DISTINCT matrix per gap
    print(f"after dt={dt:>2}: {b.round(4)}")
# after dt= 2: [0.6703 0.2968 0.0329]
# after dt=16: [0.0408 0.182  0.7772]
# after dt= 5: [0.015  0.0855 0.8995]

# sliced (2+16+5) vs one shot exp(23Q): identical
one_shot = np.array([1.0,0,0]) @ expm(Q * 23)
assert np.allclose(b, one_shot)        # slicing the gap is free

The Eigen View: Why It Decays Toward Equilibrium

There is a deeper reason long gaps drive belief toward a fixed distribution, and it is worth a paragraph because it explains the "ribbon bulges toward equilibrium" behavior you will see in every later simulation. Diagonalize Q = V Λ V⁻¹, where Λ holds the eigenvalues. Then exp(Q·dt) = V exp(Λ·dt) V⁻¹, and exp(Λ·dt) is just the scalar e^λ·dt applied to each eigenvalue.

A valid generator always has one eigenvalue equal to 0 (its eigenvector is the equilibrium / stationary distribution) and all others with negative real part. As dt grows, every e^λ·dt with λ < 0 decays to zero, leaving only the λ = 0 mode — the equilibrium. That is the precise mechanism: long gaps exponentially forget the starting belief and converge to the chain's stationary distribution, at rates set by the non-zero eigenvalues.

For our progressive chain the equilibrium is "all mass in the absorbing stage 3," which is why long gaps push belief toward stage 3 regardless of where it started. For a chain with recovery, the equilibrium would be a genuine spread across stages. Either way, the lesson is the same: the gap controls how far belief slides toward equilibrium, and the slide is exponential, not linear — which is exactly why a fixed-grid linear model gets it wrong.

# The eigen view: exp(Q*dt) = V exp(Lambda*dt) V^-1, decay set by eigenvalues.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

w, V = np.linalg.eig(Q)
print("eigenvalues:", np.sort(w.real).round(3))   # [-0.2 -0.1  0. ] one zero, rest negative

def expQt_eig(dt):
    return (V @ np.diag(np.exp(w * dt)) @ np.linalg.inv(V)).real
for dt in (5, 50):
    assert np.allclose(expQt_eig(dt), expm(Q * dt))   # matches expm
print("equilibrium (dt large):", (np.array([1.0,0,0]) @ expm(Q * 200)).round(3))
# equilibrium (dt large): [0. 0. 1.]  the 0-eigenvalue mode = all mass absorbed

The eigenvalues are exactly 0, −0.1, and −0.2 — the negatives of the rates on a triangular generator. As dt grows, the −0.1 and −0.2 modes die off (their e^λdt → 0), leaving only the zero mode: the stationary [0, 0, 1]. The slowest non-zero eigenvalue, −0.1, sets the time scale — about 1/0.1 = 10 weeks — over which belief "forgets" where it started. This is the quantitative content of "long gaps blur toward equilibrium."

Two Gaps That Sum but Split Differently

One more numeric to internalize the "gap is signal" idea. Compare two ways to span 10 weeks from a stage-1 delta: one 10-week gap, versus two 5-week gaps. Because exp(5Q)·exp(5Q) = exp(10Q), the final belief is identical — the chain only sees total elapsed time. But if a reading lands at the 5-week mark, the two-gap path lets the emission sharpen the belief midway, and the trajectories diverge.

This is the precise reason visit placement matters even when total observation time is fixed: each visit is a chance to inject an emission and re-sharpen. Two visits bracketing a likely transition pin it down; two visits both before it learn nothing about it. The transition machinery is indifferent to slicing, but the inference is not, because emissions can only enter at visits. Hold that distinction — it is the seed of the sampling-design lesson in the showcase.

# Slicing is free for the transition: split a 10w gap any way, same endpoint.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
b = np.array([1.0, 0, 0])

one = b @ expm(Q * 10)                 # one 10-week gap
two = b @ expm(Q * 5) @ expm(Q * 5) # two 5-week gaps
print(one.round(4)); print(two.round(4))
assert np.allclose(one, two)               # identical: chain sees only total time

# But inject an emission at the 5-week mark and the paths diverge:
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
mid = b @ expm(Q * 5)
mid = mid * B[:, 1]; mid /= mid.sum()   # a Mid reading at week 5
withobs = mid @ expm(Q * 5)
print("no obs @5:", one.round(3))
print("obs @5:  ", withobs.round(3))   # different: the reading re-sharpened belief

The assertion proves the transition cannot tell how you sliced the time. But the moment a reading enters at the 5-week mark, the two trajectories part — the observed path is sharper, because the emission pulled belief back toward the stage the reading favored. This is the cleanest possible demonstration that visits add information and gaps remove it, and that where you place a visit (not just how many you take) shapes the final belief.

The Hinge, in One Sentence

If Chapter 4 collapses to a single takeaway, it is this: replace the discrete HMM's one frozen transition matrix A with a function dt → exp(Q·dt), evaluated fresh at every real gap. That one substitution — a fixed matrix becomes a matrix-valued function of elapsed time — is the entire difference between a discrete HMM and a CT-HMM. Every algorithm in the remaining chapters is a familiar HMM routine with this one swap made.

Because it is a swap and not a rebuild, you inherit all your HMM intuition for free. Filtering still alternates predict and update; smoothing still runs a backward pass; Viterbi still walks a trellis. Only the edge weights change, from A to exp(Q·dt). Keep that framing and the rest of the lesson is less "new algorithms" and more "the algorithms you know, taught to respect the clock."

What Happens as the Gap Grows: A Degradation Sweep

Concept-plus-Realization: what does the predict step do as the input gap degrades toward the extremes? At dt → 0, exp(Q·dt) → I (the identity) — no time, no drift, belief unchanged. At dt → ∞, exp(Q·dt) → the stationary projection — belief forgets its start entirely and collapses to equilibrium. Everything interesting happens in between.

# The predict step across the full range of gaps, from a stage-1 delta.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
b = np.array([1.0, 0, 0])
for dt in (0.0, 0.5, 2, 5, 16, 100):
    print(f"dt={dt:>5}: {(b @ expm(Q*dt)).round(3)}")
# dt=  0.0: [1.    0.    0.   ]  identity: no drift
# dt=  0.5: [0.905 0.093 0.002]
# dt=    2: [0.67  0.297 0.033]
# dt=    5: [0.368 0.477 0.155]
# dt=   16: [0.041 0.182 0.777]
# dt=  100: [0.    0.    1.   ]  equilibrium: all absorbed

The sweep traces the full life of a belief under silence: identity at dt = 0, then progressive leakage from stage 1 into 2 and 3, then total collapse into the absorbing stage 3 at large dt. A gap of 5 weeks is the interesting middle — belief genuinely spread across all three stages, maximally uncertain. This single sweep is the engine behind every "ribbon bulges in long gaps" behavior you will see; the predict step is doing exactly this, once per real gap.

Verifying the Identity exp(Q(a+b)) = exp(Qa)exp(Qb)

The composition law we leaned on — that slicing a gap is free — rests on a property of the matrix exponential that holds because all the factors share one Q (matrices that commute). It is worth one explicit check, because for different generators it would fail.

# exp(Q(a+b)) = exp(Qa) exp(Qb) holds because both share Q (commute).
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

lhs = expm(Q * 7)                  # one 7-week gap
rhs = expm(Q * 3) @ expm(Q * 4)    # a 3w gap then a 4w gap
assert np.allclose(lhs, rhs)        # identical: slicing the gap is free

# CONTRAST: two DIFFERENT generators do NOT compose this way
Q2 = np.array([[-0.5,0.5,0],[0,-0.3,0.3],[0,0,0]])
print(np.allclose(expm(Q+Q2), expm(Q) @ expm(Q2)))   # False -- they don't commute

The assertion passes for one shared Q; the contrast at the bottom fails for two different generators, which is the general rule (exp(X+Y) = exp(X)exp(Y) only when X and Y commute). For a CT-HMM the gaps always share one Q, so we are safely in the commuting case — which is precisely why composing per-gap predictions, and slicing a long gap into shorter ones, always gives the consistent answer.

Step back and name what Chapter 4 delivered. We now have an exact, grid-free transition for any real gap, a predict step that pushes belief across it, a proof that slicing is free, and an eigen-level understanding of why long gaps blur toward equilibrium. That is the complete continuous-time transition machinery — the hinge the lesson named in its opening. From here on, every algorithm is a known HMM routine with this one operator swapped in. The hardest new idea is behind us; what remains is wiring it into filtering, smoothing, decoding, and learning.

One sentence to carry forward: the predict step is b · exp(Q·dt), it conserves total probability, it marginalizes over every hidden jump inside the gap, and it reduces to a fixed matrix when the gap is fixed. That single operation, called once per real gap, is the engine of every remaining chapter — the emission update of Chapter 5 simply alternates with it.

# The predict step in one function -- the lesson's one new operation.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
def predict(belief, dt):
    return belief @ expm(Q * dt)   # one matrix per real gap
b = np.array([1.0, 0, 0])
print(predict(b, 5).round(4))    # [0.3679 0.4773 0.1548]
print(predict(b, 2).round(4))    # [0.6703 0.2968 0.0329]

Chapter 5: The Forward Algorithm (Filtering)

In the discrete HMM you computed a forward message α by alternating "multiply by the fixed transition matrix A" and "multiply by the emission likelihood," renormalizing as you go (see hmm.html for the full discrete recursion and the scaling trick). Here the only change is the transition: instead of one fixed A, each step uses that gap's own exp(Q·dt). Everything else is identical. So this chapter goes deep on the continuous-time delta and defers the shared scaffolding to the discrete lesson.

Now stitch the pieces together. We have a way to push belief forward across a gap (Chapter 4: predict) and a way to sharpen it with a reading (Chapter 3: update). The forward algorithm just alternates them, visit after visit: stretch the belief over the gap, then squeeze it with the new measurement.

This is exactly the predict-update rhythm of the Kalman filter and Bayes filter — only the transition operator is exp(Q·dt) rather than a fixed linear motion model. If you've internalized "predict grows uncertainty, update shrinks it," you already understand the shape of this algorithm.

The result, the forward message α at each visit, is the filtered belief: the posterior over the hidden stage given all observations up to and including now. It uses the past but not the future — a future clear reading can't yet inform it. That's what we'll fix in the next chapter with smoothing.

One operation deserves special care: after multiplying by the emission likelihood, you must renormalize. The emission multiply scales the belief vector by likelihoods (numbers usually well below 1), so over many visits the raw numbers shrink toward zero and underflow. Dividing by the sum at each step keeps the belief a proper distribution — and the sum you divide by carries real information.

That normalizer is the probability of the observation given everything before it. The product of all the normalizers is the data likelihood; its logarithm, accumulated as a running sum, is the data log-likelihood — the exact quantity we'll maximize when we learn Q in Chapter 9. So renormalizing isn't just numerical hygiene; it hands you the objective function for free.

A Worked Example: Three Visits, Step by Step

Use the 3-stage progressive Q from Chapter 4. The observation sequence: at t=0 we see Lo, at t=5 we see Mid, at t=7 we see Hi. So the gaps are dt=5 (0→5) and dt=2 (5→7). The prior is π = [0.8, 0.15, 0.05] (probably stage 1).

Visit 0 (t=0, obs Lo). No gap yet — the prior is the predicted belief. Update with Lo, whose likelihood column is [0.7, 0.2, 0.05]. Numerator = π ⊙ [0.7, 0.2, 0.05] = [0.8×0.7, 0.15×0.2, 0.05×0.05] = [0.56, 0.03, 0.0025]. Sum = 0.5925. Normalize:

Gap dt=5: predict. Stretch α₀ across five weeks: α₀ · exp(5Q) = [0.3477, 0.4818, 0.1705]. The Lo-driven confidence in stage 1 (0.945) has bled toward stages 2 and 3 over the long gap.

Visit 1 (t=5, obs Mid). Mid's likelihood column is [0.25, 0.6, 0.3]. Numerator = [0.3477, 0.4818, 0.1705] ⊙ [0.25, 0.6, 0.3] = [0.0869, 0.2891, 0.0512]. Sum = 0.4272. Normalize:

Gap dt=2: predict. α₁ · exp(2Q) = [0.1364, 0.6145, 0.2491]. Short gap, small drift.

Visit 2 (t=7, obs Hi). Hi's likelihood column is [0.05, 0.2, 0.65]. Numerator = [0.1364, 0.6145, 0.2491] ⊙ [0.05, 0.2, 0.65] = [0.0068, 0.1229, 0.1619]. Normalize:

The belief marched from "certainly stage 1" through "probably stage 2" to "leaning stage 3" — exactly the disease progression, recovered from three noisy readings at irregular times. The data log-likelihood is the sum of the logs of the three normalizers (0.5925, 0.4272, and the t=7 sum), which works out to −2.6063. These five vectors and that log-likelihood are the verified checkpoints the code reproduces.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

def forward(obs, gaps, Q, B, pi):
    # obs: list of symbol indices; gaps[k] = time between obs[k-1] and obs[k]
    alphas = []; logL = 0.0
    a = pi * B[:, obs[0]]            # visit 0: prior x emission
    c = a.sum(); logL += np.log(c); a /= c
    alphas.append(a)
    for k in range(1, len(obs)):
        a = (a @ expm(Q * gaps[k-1]))  # PREDICT across the gap
        a = a * B[:, obs[k]]            # UPDATE with the reading
        c = a.sum(); logL += np.log(c); a /= c   # renormalize, accumulate logL
        alphas.append(a)
    return alphas, logL

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05])
alphas, logL = forward([0,1,2], [None,5,2], Q, B, pi)  # Lo,Mid,Hi
for a in alphas: print(a.round(4))
print("logL =", round(logL, 4))
# a0 = [0.9451 0.0506 0.0042]
# a1 = [0.2035 0.6768 0.1197]
# a2 = [0.0234 0.4214 0.5552]   logL = -2.6063

There is no single "library function" for the CT-HMM forward pass in NumPy/SciPy — the loop above is the standard implementation, with expm doing the only heavy lifting. Specialized packages (e.g. clinical-staging CT-HMM libraries) wrap exactly this recursion.

Why the Normalizer IS the Likelihood

The claim that "the product of the normalizers is the data likelihood" deserves a derivation, because it is the bridge to the learning chapter and beginners take it on faith too often. Define the unnormalized forward vector at step k, before dividing, as α̃_k. The normalizer is c_k = Σ_j α̃_k[j].

By construction, α̃_k[j] = P(observations o_0..k, hidden stage at k = j). Summing over j marginalizes out the current stage, leaving c_k = P(o_0..k) / P(o_0..k−1) once you account for the prior normalization — that is, each c_k is the probability of the new observation given all previous ones, P(o_k | o_0..k−1).

So the chain rule of probability — P(o_0..K) = P(o₀)·P(o₁|o₀)·… — is literally the product of the normalizers. Taking logs turns the product into the sum we accumulate. That is why logL += np.log(c) in the code is not an afterthought: it computes the exact objective the EM of Chapter 9 maximizes, for free, as a side effect of keeping the belief normalized.

The Underflow Problem, Quantified

Why renormalize at all, rather than carry unnormalized α̃ and divide once at the end? Because each step multiplies by an emission likelihood — numbers typically well below 1 — so α̃ shrinks geometrically. Over a long sequence it underflows to literal zero in floating point, and 0/0 at the final normalize gives NaN. The belief silently dies.

Concretely: if each step's normalizer is around 0.4, then after 800 visits the unnormalized product is 0.4⁸⁰⁰ ≈ 10⁻³¹⁸, right at the edge of double-precision (the smallest normal double is ≈ 10⁻³⁰⁸). A few visits more and it is exactly zero. Renormalizing each step keeps the working vector O(1) and pushes all the shrinkage into the log-likelihood accumulator, which handles huge negative numbers gracefully. This is the same scaling trick the discrete HMM uses; continuous time changes nothing about the numerics.

# Underflow demo: unnormalized forward dies; normalized survives.
import numpy as np
np.seterr(under='ignore')
emis = np.array([0.4, 0.3, 0.3])      # a typical per-step likelihood

raw = np.array([1.0, 0, 0]); logL = 0.0; norm = raw.copy()
for step in range(1000):
    raw = raw * emis                  # NEVER renormalized
    norm = norm * emis                # renormalized version
    c = norm.sum(); logL += np.log(c); norm /= c
print("raw (unnormalized) max:", raw.max())   # 0.0  -- underflowed to zero
print("normalized belief:    ", norm.round(3))  # still a valid distribution
print("logL (well-behaved):  ", round(logL, 2))  # large negative, but finite

The unnormalized vector is dead zeros; the normalized one is a clean distribution and the log-likelihood is a finite large-negative number. Run this in production without normalizing and your filter returns NaN after a few hundred visits — a bug that does not show up on toy data with three visits, which is exactly why it is worth demonstrating here.

The Log-Likelihood of Our Three Visits, Term by Term

Let us account for the −2.6063 log-likelihood we quoted, so it is not a magic number. The three normalizers from the worked example are the sums we computed before each normalize: c₀ = 0.5925 (visit 0), c₁ = 0.4272 (visit 1), and c₂ = the visit-2 numerator sum.

Visit 2's numerator was [0.0068, 0.1229, 0.1619], summing to c₂ = 0.2916. Now take logs: log 0.5925 = −0.5234, log 0.4272 = −0.8504, log 0.2916 = −1.2325. Add them: −0.5234 − 0.8504 − 1.2325 = −2.6063. The three per-visit "surprise" terms sum to the total data log-likelihood — exactly the verified checkpoint.

Read each term as "how surprised was the model by this reading, given everything before it." A less-likely reading contributes a more-negative log term. The total is the model's overall surprise at the whole sequence — and maximizing it (least surprised) over Q is precisely what learning does. The forward pass thus serves double duty: it produces the filtered beliefs and the scalar that scores how well the current (Q, B) explains the data.

The Predict-Update Rhythm You Already Know

If you have met the Kalman filter or Bayes filter, the forward algorithm is the same two-beat rhythm, and naming the correspondence makes it click. Predict: push belief forward through the dynamics, which spreads uncertainty — Kalman grows the covariance, the CT-HMM smears the categorical belief toward equilibrium via exp(Q·dt). Update: fold in a measurement, which shrinks uncertainty — Kalman applies the gain, the CT-HMM multiplies by the emission likelihood and renormalizes.

The only structural difference is what "uncertainty" is — a Gaussian covariance for Kalman, a categorical distribution for the CT-HMM — and what the transition operator is — a linear F for Kalman, exp(Q·dt) for the CT-HMM. The skeleton is identical, which is why Chapter 11 places them in the same family. If you trust "predict spreads, update sharpens," you already trust the forward algorithm; we have only changed the costumes.

From-Scratch: Watch Both Phases, Visit by Visit

The instrumented forward pass below prints the belief after each PREDICT and after each UPDATE separately, so you see the spread-then-sharpen rhythm in the actual numbers — not just the post-update α's we tabulated. This is the algorithm laid bare.

import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]

a = pi * B[:, obs[0]]; a /= a.sum()
print("v0 UPDATE:", a.round(3))                 # [0.945 0.051 0.004]
for k in range(1, len(obs)):
    pred = a @ expm(Q * gaps[k])             # PREDICT: spread over the gap
    print(f"v{k} PREDICT:", pred.round(3))
    a = pred * B[:, obs[k]]; a /= a.sum()  # UPDATE: sharpen with reading
    print(f"v{k} UPDATE: ", a.round(3))
# v1 PREDICT: [0.348 0.482 0.171]  spread by the 5-week gap
# v1 UPDATE:  [0.204 0.677 0.12 ]  sharpened by the Mid reading
# v2 PREDICT: [0.136 0.615 0.249]  spread by the 2-week gap
# v2 UPDATE:  [0.023 0.421 0.555]  sharpened by the Hi reading

Read the PREDICT lines and the UPDATE lines as alternating breaths. Each PREDICT raises uncertainty (the belief flattens, mass leaks to later stages); each UPDATE lowers it (the belief sharpens toward the stage the reading favors). The 5-week gap spreads more than the 2-week gap, exactly as exp(Q·dt) dictates. That alternation — spread, sharpen, spread, sharpen — is filtering, and it is the same two-beat you would see in any predict-update filter.

How Spacing Changes the Filtered Belief

Because the predict step's spread scales with the gap, the same readings at different spacings produce different filtered beliefs. Dense readings keep the belief sharp; sparse readings let it blur between updates. Let us run the identical Lo/Mid/Hi sequence at two spacings and compare the final filtered belief.

# Same readings, different gap spacing -> different filtered belief.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]

def final_belief(gaps):
    a = pi * B[:, obs[0]]; a /= a.sum()
    for k in range(1, len(obs)):
        a = (a @ expm(Q * gaps[k])) * B[:, obs[k]]; a /= a.sum()
    return a

print("dense  (gaps 1,1):", final_belief([None,1,1]).round(3))
print("sparse (gaps 5,2):", final_belief([None,5,2]).round(3))
print("very sparse(20,2):", final_belief([None,20,2]).round(3))
# dense  (gaps 1,1): leans stage 2 -- little drift between readings
# sparse (gaps 5,2): the verified [0.023 0.421 0.555]
# very sparse(20,2): more mass in stage 3 -- the long gap drifted further

With dense readings the belief barely drifts between updates, so the readings dominate and the trajectory stays sharp. With a 20-week first gap, the belief drifts far toward the absorbing late stage before the second reading even arrives — so the same Mid reading lands on a very different prior and the final belief shifts. This is the filtered-belief consequence of the sampling-design idea: when you observe, not just what you observe, shapes the answer. The forward algorithm handles it automatically because the gap lives inside expm(Q * gaps[k]).

To recap the chapter precisely: the forward algorithm maintains the filtered belief α by alternating PREDICT (α ← α · exp(Q·dt)) and UPDATE (α ← α ⊙ B[:,obs], renormalize), accumulating logL += log(Σ) at each step. The PREDICT is the only continuous-time-specific line; the UPDATE and the log-likelihood bookkeeping are identical to the discrete HMM. The filtered belief uses the past only — which is exactly the limitation the backward pass of Chapter 6 removes.

And the data log-likelihood the forward pass hands you for free is not a footnote: it is the objective that learning maximizes (Chapter 9) and the score that compares competing models (a higher logL means a better fit). Every forward pass thus does double duty — it tells you where the patient probably is and how well your model explains the data — which is why filtering is the workhorse routine the rest of the lesson builds on.

# Using logL to COMPARE two candidate Q's on the same data.
import numpy as np
from scipy.linalg import expm
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]
def logL(Q):
    a = pi * B[:,obs[0]]; L = np.log(a.sum()); a /= a.sum()
    for k in range(1,len(obs)):
        a = (a@expm(Q*gaps[k]))*B[:,obs[k]]; L += np.log(a.sum()); a /= a.sum()
    return L
Qa = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
Qb = np.array([[-0.5,0.5,0],[0,-0.3,0.3],[0,0,0]])
print("Qa logL:", round(logL(Qa),4), " Qb logL:", round(logL(Qb),4))
# the higher logL is the better-fitting rate matrix for this patient

This is the same scalar the learning chapter maximizes — here we use it merely to rank two candidate Q's, but it is exactly the objective EM climbs. That the forward pass yields both the filtered beliefs and this model-comparison score, in one sweep, is what makes it the most-used routine in the entire CT-HMM toolkit.

# The entire forward algorithm, distilled to its two lines.
#   PREDICT: a = a @ expm(Q * dt)          # spread across the gap
#   UPDATE:  a = a * B[:, obs]; a /= a.sum()  # sharpen with the reading
# accumulate logL += log(sum before normalizing) for the data likelihood.
# Everything else in this lesson wraps these two operations.

Two lines, alternated visit after visit, with the gap living only inside the PREDICT — that is filtering, and the next chapter adds a backward pass to make it smoothing.

Chapter 6: Backward Pass & Smoothing

Filtering (Chapter 5) only uses the past. But a crystal-clear reading next month should clean up our guess about last month. The backward pass carries information from the future back to each visit — and it unlocks something a discrete HMM cannot do at all: a belief at any continuous moment between visits. This is one of the two genuinely novel payoffs of continuous time, so we go slow.

The backward message β is the mirror image of the forward message, run right-to-left. Where α accumulates evidence from the start up to a visit, β accumulates evidence from the end back to a visit. It uses the very same machinery — exp(Q·dt) for transitions, B for emissions — just traversed in reverse.

Multiply them and you get the smoothed posterior γ = (α ⊙ β), normalized. This is the belief over the hidden stage at a visit given the entire record, past and future. Filtering answers "what do I believe given data so far?"; smoothing answers "what do I believe given everything I'll ever see?" Smoothing is always at least as informed.

The recursion for β, right to left, mirrors the forward predict-update but composes in the opposite order: β_k = exp(Q·dt) · (B[:,obs_k+1] ⊙ β_k+1). You first fold in the next visit's emission, then transport that backward across the gap. Initialize β at the last visit to all ones (no future evidence beyond the end).

Now the continuous-time superpower. Because the transition is defined for any duration, we can ask for the smoothed belief at a time t that lies between two visits. Propagate the smoothed belief forward from the earlier visit by exp(Q·(t−t_prev)), and you read off the belief in the dark stretch where no one looked. The belief ribbon bulges toward equilibrium in long gaps and pinches tight at each visit — a picture impossible in a fixed-grid HMM.

The honest teaching moment is not at the visits (where a clear reading dominates) but between them — and that's exactly where filtering and smoothing diverge most. We frame the worked example around the belief at t=5, midway through a long gap bracketed by an ambiguous early reading and a certain later one.

A Worked Example: The Future Fixes the Past (and the Middle)

Pin a 2-state model: Healthy (H), Diseased (D, absorbing), rate q_HD = 0.1/week, so Q = [[−0.1, 0.1], [0, 0]]. Emissions are crisp: B = [[0.9, 0.1], [0.1, 0.9]] (rows H/D, columns "looks-H"/"looks-D"). Prior [0.9, 0.1]. Two visits: t=0 we see "looks-H", t=10 we see "looks-D".

Filtered at t=0. Update prior [0.9, 0.1] with "looks-H" likelihood [0.9, 0.1]: numerator = [0.81, 0.01], sum 0.82, normalize → α₀ = [0.9878, 0.0122]. The early reading strongly says Healthy.

Backward to t=0. The future is the certain-D reading at t=10. First fold its emission: B[:,D] = [0.1, 0.9], times β₁₀ = [1, 1] gives [0.1, 0.9]. Transport back across dt=10 via exp(10Q): the H-row of exp(10Q) is [e^−1.0, 1−e^−1.0] = [0.3679, 0.6321], the D-row is [0, 1]. So β₀ = exp(10Q) · [0.1, 0.9]^T = [0.3679×0.1 + 0.6321×0.9, 0×0.1 + 1×0.9] = [0.6057, 0.9]. The future "prefers" D (0.9 > 0.6057).

Smoothed at t=0. γ₀ ∝ α₀ ⊙ β₀ = [0.9878×0.6057, 0.0122×0.9] = [0.5983, 0.0110]. Normalize → γ₀ = [0.982, 0.018]. The t=0 reading was so clear that even a certain-D future barely moves it — smoothing honestly reports the strong early evidence. (A common myth claims t=0 swings to ~[0.46, 0.54]; with a clear t=0 reading it does not.)

The real story is at t=5, mid-gap. Filtered (past only): propagate α₀ forward by exp(5Q), giving filtered@5 = [0.599, 0.401] — it still leans H, because the only evidence is the old "looks-H" reading decaying. Smoothed: multiply α@5 by β@5 (which carries the certain-D future), giving γ@5 = [0.408, 0.592] — it leans D. Between the visits, the future genuinely flips the belief. That gap, [0.599, 0.401] vs [0.408, 0.592], is the smoothing payoff made numeric.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

def backward(obs, gaps, Q, B):
    Nv = len(obs); N = Q.shape[0]
    beta = [None]*Nv; beta[-1] = np.ones(N)        # no future beyond the end
    for k in range(Nv-2, -1, -1):
        nxt = B[:, obs[k+1]] * beta[k+1]          # fold next emission
        beta[k] = expm(Q * gaps[k+1]) @ nxt        # transport back across gap
    return beta

def smooth(alphas, beta):
    g = [a * b for a, b in zip(alphas, beta)]
    return [gi / gi.sum() for gi in g]

def interp_belief(gamma_prev, Q, t_into_gap):
    return gamma_prev @ expm(Q * t_into_gap)         # between-visit belief

Q = np.array([[-0.1,0.1],[0,0]]); B = np.array([[0.9,0.1],[0.1,0.9]])
# with alphas from forward(): gamma0 = [0.982, 0.018]; filtered@5=[0.599,0.401]; gamma@5=[0.408,0.592]

Again there is no off-the-shelf "CT-HMM smooth" call — these two short loops are the implementation, with expm the only library dependency. The interp_belief helper is the piece with no discrete-HMM analogue: it answers "what stage at week 5?" exactly, no grid required.

What the Backward Message Actually Means

The forward message α_k[j] = P(observations o_0..k, stage at k = j) is a joint probability. The backward message is the conditional complement: β_k[j] = P(future observations o_k+1..K | stage at k = j). It answers "if I were in stage j right now, how probable is everything I will see later?" Unlike α, it is not a distribution over stages — it is a likelihood of the future, one number per possible current stage.

Why does the product α ⊙ β give the smoothed posterior? Because α_k[j] · β_k[j] = P(all observations, stage at k = j) by the conditional-independence of past and future given the present (the Markov property). Normalizing over j divides out P(all observations), leaving P(stage at k = j | all observations) — the smoothed belief. The Markov property is the load-bearing assumption: it is what lets past and future combine by simple multiplication.

This is also why β transports backward with exp(Q·dt) applied on the left (matrix times column vector) while α used it on the right (row vector times matrix). The forward message flows with the arrow of time; the backward message flows against it. Same matrix exponential, opposite multiplication side — a small but crucial bookkeeping detail the code makes explicit.

The Between-Visit Belief, Derived

The genuinely new continuous-time capability is reading belief at a time t strictly between visits, where no discrete HMM can speak. The recipe: take the smoothed belief γ at the earlier bracketing visit, and propagate it forward by the partial gap (t − t_prev) using exp(Q·(t−t_prev)). For the filtered between-visit belief you propagate the filtered α instead.

For the smoothed version, the fully-correct interpolation also conditions on the future visit — you propagate γ forward by the partial gap and reconcile with the backward message at t. In the simple bracketed-by-two-visits case the dominant effect is the forward propagation of the smoothed start, which is what our worked t=5 numbers use. The point is that the belief is defined and computable at every real time, not just at the visits — the ribbon is continuous because exp(Q·τ) is defined for every τ.

# The full filtered-vs-smoothed contrast at t=5, reproducing the numbers.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.1,0.1],[0,0]]); B = np.array([[0.9,0.1],[0.1,0.9]])
prior = np.array([0.9, 0.1])

# filtered at t=0: prior x emission(looks-H)
a0 = prior * B[:, 0]; a0 /= a0.sum()       # [0.9878 0.0122]
# filtered at t=5: just propagate a0 forward 5 weeks (no obs yet)
filt5 = a0 @ expm(Q * 5)                    # [0.599 0.401] leans H

# backward at t=5: future is certain-D reading at t=10 (dt=5 from t=5)
beta5 = expm(Q * 5) @ (B[:, 1] * np.ones(2))  # fold D-emission, transport back
g5 = filt5 * beta5; g5 /= g5.sum()           # smoothed at t=5
print("filtered@5:", filt5.round(3))   # [0.599 0.401]  leans HEALTHY
print("smoothed@5:", g5.round(3))      # [0.408 0.592]  leans DISEASED

The two beliefs at the very same instant disagree about which stage is more likely — filtered says Healthy (0.599), smoothed says Diseased (0.592). The future certain-D reading, carried back by β, flips the verdict in the unobserved middle. That flip is invisible to filtering and impossible to even ask about in a discrete HMM with no notion of "between." It is the chapter's payoff in two printed vectors.

Filtering vs Smoothing: When to Use Which

Filtering is for real-time decisions; smoothing is for analysis after the fact, when all the data is in. They are not competitors — they answer different questions, and smoothing is strictly more informed only because it has strictly more data (the future). The between-visit belief is the uniquely continuous-time offering: a principled answer to "what was happening when no one was looking," with honest uncertainty that bulges in long gaps and pinches at visits.

From-Scratch: Full Forward-Backward, End to End

Let us assemble the complete smoother — forward pass, backward pass, and their product — on the 3-stage progressive model with the Lo/Mid/Hi sequence, so you see γ differ from the filtered α at every visit, not just the toy 2-state case.

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]

def forward():
    a = pi * B[:, obs[0]]; A = [a / a.sum()]
    for k in range(1, len(obs)):
        a = (A[-1] @ expm(Q * gaps[k])) * B[:, obs[k]]
        A.append(a / a.sum())
    return A

def backward():
    Nv = len(obs); beta = [None]*Nv; beta[-1] = np.ones(3)
    for k in range(Nv-2, -1, -1):
        beta[k] = expm(Q * gaps[k+1]) @ (B[:, obs[k+1]] * beta[k+1])
    return beta

A, be = forward(), backward()
for k in range(len(obs)):
    g = A[k] * be[k]; g /= g.sum()
    print(f"visit {k}: filtered={A[k].round(3)}  smoothed={g.round(3)}")
# visit 0: filtered=[0.945 0.051 0.004]  smoothed=[0.876 0.116 0.008]
# visit 1: filtered=[0.204 0.677 0.12 ]  smoothed=[0.142 0.706 0.152]
# visit 2: filtered=[0.023 0.421 0.555]  smoothed=[0.023 0.421 0.555]

Two things to read off. At visit 0 the smoothed belief pulls toward later stages (0.945 → 0.876 on stage 1) because the future Mid and Hi readings reveal the patient was already progressing — the future corrected the past. At the last visit, filtered and smoothed are identical, because there is no future beyond the end to add information (β is all ones there). Smoothing changes the early and middle estimates most; it cannot improve on the final one.

From-Scratch: The Continuous Ribbon Between Visits

Finally, the piece with no discrete analogue: sample the smoothed belief on a fine grid inside a gap to draw the ribbon. Note the grid here is only for plotting — the model itself never used one; we are just evaluating the exact continuous belief at display resolution.

# Evaluate the continuous between-visit belief for plotting (2-state H/D).
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.1,0.1],[0,0]])
gamma0 = np.array([0.982, 0.018])     # smoothed belief at the t=0 visit

for tau in np.linspace(0, 5, 6):       # plot points inside the gap
    belief = gamma0 @ expm(Q * tau)  # exact belief tau weeks into the gap
    print(f"t={tau:.0f}: P(H)={belief[0]:.3f}")
# t=0: P(H)=0.982   t=1: P(H)=0.889   t=2: P(H)=0.804
# t=3: P(H)=0.727   t=4: P(H)=0.658   t=5: P(H)=0.595  (smoothly decaying)

The belief decays smoothly from the pinned visit value toward the equilibrium — a continuous curve we can evaluate at any resolution because exp(Q·τ) is defined for every τ. That smooth decay is the ribbon you see bulging through long gaps in the showcase. No discrete HMM can produce it; it has no concept of a belief between its fixed ticks. This is the concrete, code-level meaning of "the only change is the per-gap matrix exponential" — it quietly hands you a belief at every instant.

Smoothing Never Increases Uncertainty

A useful guarantee to internalize: at any visit, the smoothed belief is at least as certain as the filtered belief, because it conditions on strictly more data (the future). We can measure "certainty" by entropy — lower entropy means a sharper, more confident distribution — and check that smoothing's entropy never exceeds filtering's at the same visit.

# Smoothing conditions on more data, so it is never less certain.
import numpy as np
def entropy(p):
    p = p[p > 0]; return float(-(p * np.log(p)).sum())

# from the 3-stage forward-backward above:
filtered = [np.array([0.945,0.051,0.004]),
            np.array([0.204,0.677,0.12]),
            np.array([0.023,0.421,0.555])]
smoothed = [np.array([0.876,0.116,0.008]),
            np.array([0.142,0.706,0.152]),
            np.array([0.023,0.421,0.555])]
for k in range(3):
    print(f"visit {k}: H_filt={entropy(filtered[k]):.3f}  H_smooth={entropy(smoothed[k]):.3f}")
# at the last visit the two are equal (no future to add);
# earlier, smoothing is generally as sharp or sharper once the future is informative

The guarantee is exact at the final visit (filtered equals smoothed, so entropies match) and holds in expectation elsewhere: more data cannot, on average, make you less certain. This is the formal sense in which "the future fixes the past" — it does not just shift the belief, it tightens it, because the backward message contributes genuine evidence. When you see the smoothed ribbon pinch tighter than the filtered one in the showcase, this entropy inequality is what you are watching.

The Backward Message Is a Likelihood, Not a Belief

A subtlety worth a code check: unlike the forward α (which we normalize into a distribution), the backward β is a likelihood of the future and need not sum to anything in particular. Normalizing it would discard information. We only normalize after forming the product α ⊙ β. The snippet confirms β entries are not a distribution, yet the smoothed product is.

import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.1,0.1],[0,0]]); B = np.array([[0.9,0.1],[0.1,0.9]])

# backward message at t=0 for the certain-D future at t=10
beta0 = expm(Q * 10) @ (B[:, 1] * np.ones(2))
print("beta0 =", beta0.round(4), "  sum =", round(beta0.sum(), 4))
# beta0 = [0.6057 0.9   ]   sum = 1.5057   (NOT a distribution -- a likelihood)

alpha0 = np.array([0.9878, 0.0122])      # filtered belief at t=0
gamma0 = alpha0 * beta0; gamma0 /= gamma0.sum()   # NOW normalize
print("gamma0 =", gamma0.round(3), "  sum =", round(gamma0.sum(), 3))
# gamma0 = [0.982 0.018]   sum = 1.0  (a proper distribution)

The β vector sums to 1.5057 — meaningless as a probability, perfectly meaningful as a per-stage future-likelihood. Only the product with α, once normalized, is the smoothed belief [0.982, 0.018]. Confusing β for a distribution and normalizing it early is a common smoothing bug; the correct discipline is "multiply first, normalize once." With that, the forward-backward smoother is complete and correct, and the between-visit ribbon falls out for free from expm(Q·τ).

Chapter 6 delivered the two genuinely new continuous-time payoffs. First, smoothing: the backward message β carries future evidence back to each visit, and γ ∝ α ⊙ β is the all-data posterior — strictly more informed than filtering. Second, between-visit interpolation: because exp(Q·τ) is defined for every τ, we can read the belief at any instant in a gap, drawing a continuous ribbon no discrete HMM can produce. The t=5 contrast — filtered [0.599, 0.401] leaning H versus smoothed [0.408, 0.592] leaning D — is the chapter's thesis in two vectors: the future genuinely rewrites the unobserved middle.

# Filtered vs smoothed at every point in the gap -- the diverging ribbons.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.1,0.1],[0,0]])
a0 = np.array([0.9878, 0.0122])               # filtered at t=0
g0 = np.array([0.982, 0.018])                 # smoothed at t=0
for t in (0, 2, 5):
    filt = a0 @ expm(Q * t)
    print(f"t={t}: filtered P(H)={filt[0]:.3f}")    # decays 0.988 -> ... -> 0.599
# the smoothed ribbon, carrying the certain-D future, sits BELOW this,
# crossing under 0.5 by t=5 -- the future pulling the middle toward D.

The filtered ribbon decays gently from the strong t=0 reading; the smoothed ribbon, pulled by the certain-D reading at t=10, dips faster and crosses below 0.5 in the middle of the gap. Two curves over the same interval, disagreeing most where no one looked — the visual signature of smoothing's value, and the picture the showcase animates.

# Smoothing, distilled:
#   beta[k] = expm(Q*gap) @ (B[:,obs_next] * beta[k+1])  # carry future back
#   gamma   = normalize(alpha * beta)                    # all-data posterior
#   between-visit: gamma_prev @ expm(Q * tau)            # belief at any tau
# beta is a LIKELIHOOD (don't normalize it); normalize only the product.

Chapter 7: Decoding the Most Likely Trajectory

In the discrete HMM, Viterbi finds the single most likely sequence of hidden states by carrying a best-path score forward through a trellis, using the fixed transition matrix A on every edge and recording back-pointers (full derivation in hmm.html). Here the only change is the edge weight: each gap's transition is exp(Q·dt) instead of one fixed A. We defer the shared trellis scaffolding to the discrete lesson and focus on the continuous-time edge and the path-vs-marginal trap.

Smoothing (Chapter 6) gives a belief at each visit. But sometimes you want one committed story: "stage 1, then stage 1, then stage 2, then stage 3." That single best sequence is what decoding produces.

Continuous-time Viterbi is discrete Viterbi with the fixed A swapped for per-gap exp(Q·dt). We work in log-space (sums instead of products) for numerical stability, exactly as in the discrete case. Each trellis edge from stage i to stage j across a gap carries weight log P(dt)_ij from exp(Q·dt), plus the destination's emission log-likelihood.

Here is the trap that survives the jump to continuous time. The most likely PATH is not the chain of most-likely STATES. If you take smoothing's highest-probability stage at each visit independently and string them together, you can produce a path the chain assigns near-zero joint probability — for instance, one that demands a backward jump that a progressive (no-recovery) Q forbids outright.

Why do they differ? Smoothing maximizes each visit's marginal probability separately. Viterbi maximizes the joint probability of the whole path. Viterbi will happily reject a stage that's locally most likely if reaching it requires an improbable transition, choosing instead a globally coherent story. They answer different questions, and both are correct answers to their own question.

The recursion: initialize δ₀ = logπ + log B[:,obs₀]. For each gap, δ_k[j] = max_i(δ_k−1[i] + log P(dt)_ij) + log B[j, obs_k], storing the argmax i as a back-pointer. At the end, take the best final stage and follow back-pointers to recover the path.

A Worked Example: One Trellis Cell, in Full

Use the 3-stage progressive Q, visits at gaps dt = [5, 2], observations ['Lo', 'Mid', 'Hi'], prior π = [0.8, 0.15, 0.05]. Work in logs.

Initialize. δ₀ = logπ + log B[:,Lo]. With B[:,Lo] = [0.7, 0.2, 0.05]: δ₀ = [log 0.8 + log 0.7, log 0.15 + log 0.2, log 0.05 + log 0.05] = [−0.580, −3.507, −5.991].

Transition across dt=5. From exp(5Q), the stage-1 row gives P₁₁=0.3679, P₁₂=0.4773, P₁₃=0.1548, so log P = [−1.000, −0.740, −1.866]. (Stages 2 and 3 can't reach stage 1 or 2 backward in this progressive Q, so their entries are −∞ / very negative.)

Compute δ₁[stage 2] in full. We need max over starting stage i of (δ₀[i] + log P(5)_i,2), then add log B[2, Mid] = log 0.6 = −0.511. The candidates into stage 2: from stage 1, −0.580 + log 0.4773 = −0.580 + (−0.740) = −1.320; from stage 2, −3.507 + log P₂₂(5). The stage-1 origin wins (−1.320 is the larger/less-negative). So δ₁[2] = −1.320 + (−0.511) = −1.831, with back-pointer = stage 1.

Carry to δ₂. Repeat with log exp(2Q) and log B[:,Hi]. Backtracing from the best final stage yields the path [1, 2, 3] — the patient progressed stage 1 → 2 → 3, the only story a no-recovery chain allows that also fits Lo→Mid→Hi. The greedy per-visit argmax happens to agree here, but tweak the emissions or rates and it can demand an impossible stage-3→stage-2 reversal that Viterbi refuses.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

def viterbi(obs, gaps, Q, B, pi):
    Nv = len(obs); N = Q.shape[0]
    eps = 1e-300
    delta = np.log(pi + eps) + np.log(B[:, obs[0]] + eps)
    bptr = []
    for k in range(1, Nv):
        logP = np.log(expm(Q * gaps[k]) + eps)   # per-gap edge weights
        scores = delta[:, None] + logP          # [from i, to j]
        bp = scores.argmax(axis=0)               # best predecessor per j
        delta = scores.max(axis=0) + np.log(B[:, obs[k]] + eps)
        bptr.append(bp)
    path = [int(delta.argmax())]                # backtrace
    for bp in reversed(bptr): path.append(int(bp[path[-1]]))
    return path[::-1]

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
print(viterbi([0,1,2], [None,5,2], Q, B, np.array([0.8,0.15,0.05])))
# -> [0, 1, 2]   i.e. stage 1 -> stage 2 -> stage 3

The library shortcut: log-space and logsumexp-style stability come from scipy.special, but the decode itself is the loop above — the continuous-time change is the single line logP = np.log(expm(Q * gaps[k])).

Why Log-Space, and Why max Not sum

Two design choices in Viterbi deserve unpacking. First, why logs? A path's probability is a product of many transition and emission factors, each below 1; the product underflows for long paths, exactly as in the forward algorithm. Taking logs turns products into sums, which never underflow — you accumulate large negative numbers instead of vanishing tiny ones. Every comparison Viterbi makes is between sums of log-probabilities.

Second, why max where the forward algorithm used sum? Because the questions differ. The forward algorithm asks "what is the total probability of being in stage j, summed over all ways to get there" — so it sums over predecessors. Viterbi asks "what is the probability of the single best way to be in stage j" — so it maxes over predecessors and remembers which one won. Sum marginalizes; max optimizes. One number, one operator, all the difference between a posterior and a best path.

The back-pointer is the other half of Viterbi: at each step you record which predecessor i achieved the max into each j. Forward never needs this — it sums, so there is no "winner" to remember. After the forward sweep, Viterbi starts at the best final stage and walks the back-pointers in reverse to reconstruct the single most probable path. That reconstruction is why Viterbi returns a path while the forward algorithm returns a belief.

A Worked Example: When Path and Marginal Disagree

Let us build a case where the greedy "pick each visit's most-likely stage" path differs from Viterbi's joint-best path — the classic trap, alive in continuous time. Use a progressive Q (no recovery: stage 3 cannot return to 2). Suppose smoothing gives per-visit marginals: visit 0 most-likely stage 1, visit 1 most-likely stage 3, visit 2 most-likely stage 2.

The greedy path reads off [1, 3, 2]. But the chain forbids 3→2 (no recovery) — that transition has probability zero, so the greedy "path" has zero joint probability. It is not a path the model can produce at all. Viterbi refuses it: it will instead choose, say, [1, 2, 2] or [1, 2, 3] — whichever globally-coherent story has the highest joint probability — even though stage 3 was locally most-likely at visit 1.

This is the precise sense in which "the most likely path is not the chain of most likely states." The marginal at each visit ignores what the neighbors demand; the joint path respects every transition. With a progressive Q the trap is stark because backward jumps are outright forbidden, but it arises with any Q whenever a locally-attractive stage requires an improbable transition to reach. Viterbi is the tool when you need one coherent story; smoothing is the tool when you need honest per-visit uncertainty.

# Greedy-marginal path can be IMPOSSIBLE; Viterbi returns a legal best path.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

P = expm(Q * 5)
print("P(3 -> 2 over 5w):", P[2, 1])   # 0.0  -- a backward jump is impossible
# so a greedy path that demands stage 3 then stage 2 has joint prob 0:
greedy_logp = np.log(P[2, 1] + 1e-300)
print("greedy [.,3,2] log-prob:", round(greedy_logp, 1))  # ~ -690, i.e. ruled out
# Viterbi (from earlier) never proposes 3->2; it returns a legal monotone path.

The transition probability for the forbidden 3→2 jump is exactly zero, so its log is −∞ (here clamped to about −690 by the epsilon). Viterbi's max will never pick a predecessor through a −∞ edge, so it structurally cannot return an impossible path — while the greedy marginal readout cheerfully does. That guarantee — Viterbi paths are always legal trajectories — is exactly why you decode with Viterbi rather than stringing together argmaxes.

The Trellis, Drawn in Numbers

Let us complete the worked trellis cell we started, carrying all three stages forward so you see the full δ table, not just one entry. Recall δ₀ = [−0.580, −3.507, −5.991] from logπ + log B[:,Lo]. The first transition is across dt = 5 weeks.

From exp(5Q), the stage-1 row is [0.3679, 0.4773, 0.1548], so log P(5) from stage 1 is [−1.000, −0.740, −1.866]. Stage 2 can only reach stages 2 and 3 (progressive), and stage 3 only itself. Adding δ₀[i] to each log-edge and maxing over the origin i, then adding the destination emission log B[:,Mid] = [log 0.25, log 0.6, log 0.3] = [−1.386, −0.511, −1.204], gives δ₁.

For destination stage 2: candidates are from stage 1 (−0.580 + (−0.740) = −1.320) and from stage 2 (−3.507 + log P₂₂(5)). The stage-1 origin wins at −1.320; add the Mid emission (−0.511) to get δ₁[2] = −1.831, back-pointer stage 1. Repeating for stages 1 and 3 fills the column, and the same procedure across dt = 2 with log B[:,Hi] fills δ₂. Backtracing the winners yields [1, 2, 3].

# Full Viterbi vs greedy-marginal on the worked example.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]; eps = 1e-300

# Viterbi
delta = np.log(pi+eps) + np.log(B[:,obs[0]]+eps); bptr = []
for k in range(1, len(obs)):
    logP = np.log(expm(Q*gaps[k])+eps); sc = delta[:,None] + logP
    bptr.append(sc.argmax(0)); delta = sc.max(0) + np.log(B[:,obs[k]]+eps)
path = [int(delta.argmax())]
for bp in reversed(bptr): path.append(int(bp[path[-1]]))
print("Viterbi path:", path[::-1])      # [0, 1, 2] = stage 1,2,3
print("best log-prob:", round(float(delta.max()), 3))

The full table confirms the single best path is stage 1 → 2 → 3 with the log-probability printed, and that the back-pointers traced cleanly through legal transitions. On this gentle example the greedy per-visit argmax happens to agree, but the moment you tune the emissions so a locally-attractive stage demands a backward jump, Viterbi and greedy split — and only Viterbi stays legal. That is the whole reason the algorithm carries back-pointers instead of just taking column maxima.

Three Different Questions, Three Different Tools

It is easy to blur filtering, smoothing, and decoding together; the table pins down exactly what each returns and when to reach for it. They are not refinements of one another — they answer genuinely different questions about the hidden chain.

If you want a per-visit belief with uncertainty, smooth. If you want a single committed storyline for reporting or simulation, decode. If you must act before the future arrives, filter. Choosing the wrong one produces subtly wrong answers that still look plausible — for instance, reporting the greedy chain of smoothed argmaxes as "the patient's trajectory" can assert an impossible disease course. Match the tool to the question, and the continuous-time machinery (per-gap exp(Q·dt)) is shared across all three.

The deeper unity: all three are the same trellis walked differently. Filtering sweeps left-to-right summing; smoothing adds a right-to-left sweep and multiplies; Viterbi sweeps left-to-right maxing with back-pointers. One data structure, three traversals, three answers — and the only continuous-time change in any of them is that the trellis edges carry exp(Q·dt) instead of a fixed A.

A Tuned Example Where Greedy Goes Wrong

Let us construct the case where greedy and Viterbi genuinely disagree, so the trap is not just a warning but a concrete computation. Keep the progressive Q, but imagine an emission that makes stage 3 locally attractive at an early visit — say a misleadingly low reading at visit 1. The per-visit smoothed marginal might then favor stage 3 at visit 1, while visit 2's reading favors stage 2.

The greedy chain reads [stage 1, stage 3, stage 2]. But the progressive Q forbids 3→2: that edge has probability zero. So the greedy "trajectory" is impossible — it asserts a recovery the model rules out. Viterbi, maximizing the joint path probability, refuses the −∞ edge and returns instead [stage 1, stage 2, stage 2] or [stage 1, stage 2, stage 3] — whichever legal path scores highest. The locally-attractive stage 3 at visit 1 is rejected because reaching it and then continuing is jointly improbable.

# Greedy marginal can assert a forbidden recovery; Viterbi cannot.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])

# suppose smoothing's per-visit argmax came out [0, 2, 1] (1->3->2)
greedy = [0, 2, 1]
# joint log-prob of that path needs the 3->2 edge over the second gap:
P2 = expm(Q * 2)
print("P(stage3 -> stage2 over 2w):", P2[2, 1])      # 0.0
print("=> greedy path 1->3->2 has joint probability 0 (impossible)")
# Viterbi structurally avoids this; its returned path is always legal.

The forbidden edge has probability exactly zero, so the greedy path's joint probability is zero — it is not a trajectory the chain can produce. This is the trap in its starkest form: a progressive disease cannot un-progress, yet stringing together per-visit argmaxes can claim it did. Viterbi's joint optimization is immune by construction. Whenever you need a single coherent story — for a patient narrative, a simulation, or a report — decode it; never assemble it from independent marginals.

From-Scratch: Viterbi and Greedy Side by Side

Let us run both decoders on the same inputs and print them together, so the disagreement (when it occurs) is right there in the output. The greedy decoder takes the argmax of the smoothed marginal at each visit; Viterbi maximizes the joint path. On benign inputs they match; tune the emissions and they split.

import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]; eps = 1e-300

def greedy_smoothed():
    a = pi * B[:,obs[0]]; A = [a/a.sum()]
    for k in range(1,len(obs)):
        a = (A[-1]@expm(Q*gaps[k]))*B[:,obs[k]]; A.append(a/a.sum())
    be = [None]*len(obs); be[-1]=np.ones(3)
    for k in range(len(obs)-2,-1,-1):
        be[k] = expm(Q*gaps[k+1])@(B[:,obs[k+1]]*be[k+1])
    return [int((A[k]*be[k]).argmax()) for k in range(len(obs))]

print("greedy (per-visit argmax):", greedy_smoothed())   # [0, 1, 2] here
# Viterbi (computed earlier): [0, 1, 2] -- they agree on this benign input,
# but greedy can return an illegal path when a forbidden edge is locally tempting.

On this gentle Lo/Mid/Hi sequence both decoders return stage 1→2→3 — agreement is common when the data is clean. The value of keeping them separate is for the adversarial cases: the greedy decoder has no mechanism to reject a forbidden transition, while Viterbi's joint max structurally cannot traverse a zero-probability edge. When in doubt, decode with Viterbi; it is the only one that guarantees a legal trajectory.

To recap the chapter: decoding returns the single most likely whole-path trajectory, computed by a log-space Viterbi sweep whose only continuous-time change is per-gap exp(Q·dt) edge weights and back-pointers to reconstruct the path. The key trap — the most-likely path is not the chain of most-likely states — persists in continuous time and is starkest under a progressive Q that forbids recovery. Use smoothing for per-visit uncertainty, Viterbi for one coherent story, and never confuse the two. With filtering, smoothing, and decoding all in hand, the showcase puts them on one living timeline.

# The Viterbi edge weight: per-gap log-transition + destination emission.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])
eps = 1e-300; dt = 5; obs_next = 1   # Mid at the next visit
logP = np.log(expm(Q * dt) + eps)         # per-gap edge weights [from i, to j]
logB = np.log(B[:, obs_next] + eps)         # destination emission log-likelihood
print("edge 1->2 weight:", round(float(logP[0,1] + logB[1]), 3))   # -1.251
print("edge 1->3 weight:", round(float(logP[0,2] + logB[2]), 3))   # more negative
# Viterbi picks the predecessor maximizing delta[i] + this edge weight.

The 1→2 edge over five weeks scores higher than 1→3 because reaching stage 3 in five weeks requires passing through stage 2 — a less probable double transition. Viterbi sums these log-edge weights along candidate paths and keeps the maximum at each node; the per-gap exp(Q·dt) is the only thing that changes versus a discrete trellis, and it changes only the edge weights, never the algorithm.

# Decoding, distilled:
#   delta_0 = log(pi) + log B[:,obs0]
#   delta_k[j] = max_i( delta_{k-1}[i] + log[expm(Q*dt)]_ij ) + log B[j,obs_k]
#               store argmax i as back-pointer
#   path: start at argmax(delta_last), follow back-pointers in reverse
# maximizes JOINT path prob -> always a legal trajectory (no forbidden jumps).
# contrast: greedy = [argmax(smoothed_k) for each k] CAN be illegal.
# the ONLY continuous-time change vs discrete Viterbi: expm(Q*dt) edges.
# work in log-space (sums, not products) for numerical stability.

With filtering, smoothing, and decoding all built — each the same trellis walked summing, multiplying, or maxing — the inference half of the lesson is complete. The showcase next runs all three live on one irregular timeline, and then the final chapters learn Q from data and place the model in its family.

One closing caution worth repeating: never report the chain of per-visit smoothed argmaxes as "the trajectory." That greedy assembly can assert an impossible course; only Viterbi's joint optimization guarantees a legal one. When the deliverable is a single storyline, decode it.

Chapter 8: Showcase — The Living CT-HMM

This is the payoff. Everything you've built — the hidden chain, emissions, per-gap exp(Q·dt), forward filtering, backward smoothing — runs together on one irregular timeline. You are now the modeler: dial the progression rates, the measurement noise, and where the clinic visits fall, and watch truth and belief co-evolve.

The display has three stacked layers over a shared, to-scale time axis. Top: the true hidden continuous-time stage path (a color-coded step function). Middle: the noisy observations as dots at the irregular visit times, colored by emitted symbol. Bottom: the smoothed posterior belief ribbon, computed by full forward-backward, continuous across time — including between visits.

Because the gaps are drawn to scale, a 16-week dropout is visibly eight times wider than a 2-week interval, and the belief ribbon's smear scales with it. Starve the model with one long gap and watch the ribbon balloon toward the stationary distribution; then place a visit and watch it collapse to a sharp belief. The live readout shows the current data log-likelihood.

Try the experiments the sim invites. Push the emission noise up and the ribbon never gets crisp even at visits — the readings stop being informative. Speed up q₁₂ and the true path races through stage 1; the belief has to chase it. Cluster your visits early and leave a long tail unobserved, and the ribbon's honesty becomes visible: it admits how little it knows in the dark.

The default scenario anchors the lesson. A 3-stage progressive Q with q₁₂=0.2, q₂₃=0.1; visits at weeks [0, 2, 4, 20, 22] — a tight early cluster, a 16-week dropout, then two more. Walk the dropout: the smoothed belief at week 4 sits around [0.1, 0.6, 0.3]; over the 16-week gap, exp(Q·16) pushes it toward roughly [0.02, 0.25, 0.73] as mass piles into the absorbing late stage; the week-20 "Hi" reading then confirms close to [0.0, 0.1, 0.9].

The deepest lesson here isn't an equation — it's about sampling design. The simulation makes vivid that a visit placed right after a likely transition is worth far more than one during a stable stretch, and that during a long unobserved gap the belief is genuinely uncertain no matter how good your Q is. The information simply isn't there to recover; the model's job is to be honest about that.

No quiz here — the simulation is the test. If you can predict, before you drag a slider, which way the ribbon will move and why, you understand the CT-HMM. Push q₂₃ down and the late stages should fill more slowly; crank the noise and even the at-visit belief should refuse to commit. Make those predictions, then check them.

Chapter 9: Learning Q from Scattered Visits (CT-EM)

So far Q was handed to us. Now we learn it from data alone — only the irregular readings, no stage labels. This is EM again (E-step: infer the hidden states; M-step: re-estimate the parameters), and it shares its skeleton with discrete Baum-Welch and with the EM for Gaussian mixtures. But it is genuinely harder, and the reason is the whole point of this chapter.

In discrete Baum-Welch, the M-step is easy: you count expected transitions i→j and divide by expected visits to i. The sufficient statistics are simple counts. In continuous time that breaks, because we never see the jumps (Chapter 2). A gap that starts in stage i and ends in stage j might have gone i→j directly, or i→k→j, or jumped back and forth — all invisible.

What a rate actually needs is two continuous-time quantities, per gap: the expected total time spent in each stage and the expected number of i→j jumps. A rate is jumps-per-unit-time, so its denominator is time, not a step count. Both quantities are integrals over all the unobserved jump times within the gap — and that sounds intractable.

The beautiful fact, the Van Loan trick: both integrals fall out of one larger matrix exponential. Build a 2N×2N block matrix with Q on both diagonal blocks and an indicator matrix E in the top-right block. Exponentiate it over the gap. The top-right block of the result holds exactly the expected statistic you want — expected time-in-stage if E marks a diagonal, expected jump-mass if E marks an off-diagonal transition.

The M-step then mirrors the discrete one but with the right denominators: q_ij^new = (Σ expected i→j jumps) / (Σ expected time in i), summed over all gaps and patients, with the diagonal set to −(row sum) to keep Q a valid generator. Iterate E-step / M-step and the data log-likelihood climbs monotonically.

A Worked Example: One Gap, the 4×4 Built by Hand

Two-state H/D, q_HD unknown, current guess q_HD = 0.08. One gap of dt = 10 weeks, with the smoothed posterior starting on H and ending on D. We want the expected time-in-H and the expected number of H→D jumps over this gap.

The generator at the current guess is Q = [[−0.08, 0.08], [0, 0]]. The H→D indicator is E = [[0, 1], [0, 0]] (a single 1 in the H→D slot). Stack them into the 4×4:

Now exponentiate A·10. The top-left 2×2 block is just exp(Q·10) (the ordinary transition table). The top-right 2×2 block is the Van Loan integral. Its (H, D) entry comes out to 6.8834 — this is the expected H→D "jump mass" accumulated over the gap (before weighting by the start/end posteriors).

We can sanity-check the diagonal version: the expected time spent in H over a 10-week gap, starting in H, has the closed form (1 − e^−q·dt)/q = (1 − e^−0.8)/0.08 = (1 − 0.4493)/0.08 = 0.5507/0.08 = 6.8834 weeks. The same number falls out of the block-matrix corner with the diagonal indicator — an independent confirmation that the Van Loan machinery is correct.

The M-step. Over many gaps and patients, accumulate the expected H→D jumps in the numerator and the expected time-in-H in the denominator: q_HD^new = (Σ expected H→D jumps) / (Σ expected time in H). Starting from the guess 0.08 and iterating, this ratio climbs toward the true generating rate of 0.10 across, say, 50 patients — and the data log-likelihood increases at every iteration. The simulator above traces exactly that climb.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

def vanloan_corner(Q, E, dt):
    N = Q.shape[0]
    A = np.zeros((2*N, 2*N))
    A[:N, :N] = Q; A[N:, N:] = Q; A[:N, N:] = E   # block matrix [[Q,E],[0,Q]]
    M = expm(A * dt)
    return M[:N, N:]                                   # the top-right corner

q = 0.08
Q = np.array([[-q, q], [0, 0]])
E_jump = np.array([[0.0, 1.0], [0, 0]])         # H->D indicator
E_time = np.array([[1.0, 0], [0, 0]])             # time-in-H indicator
print("H->D jump mass:", vanloan_corner(Q, E_jump, 10)[0,1])  # 6.8834
print("time in H:    ", vanloan_corner(Q, E_time, 10)[0,0])  # 6.8834
print("closed form:  ", (1 - np.exp(-0.8)) / 0.08)            # 6.8834

def m_step(exp_jumps, exp_time):              # off-diagonal rates, then fix diagonal
    Qn = exp_jumps / exp_time[:, None]
    np.fill_diagonal(Qn, 0);
    for i in range(Qn.shape[0]): Qn[i,i] = -Qn[i].sum()
    return Qn

The "library one-liner" is again expm — applied to the stacked block matrix. SciPy has no dedicated CT-EM, but scipy.linalg.expm on the 2N×2N block is the entire trick; everything else is bookkeeping (weight each gap's corner by the smoothed start/end posteriors, sum, then form the M-step ratio).

The E-Step: Weighting the Corner by the Posterior

The Van Loan corner gives the expected statistic conditioned on a specific start and end stage. But in a CT-HMM we do not know the start and end stages of a gap — they are hidden. So the E-step weights each (start i, end j) corner by the smoothed probability that the gap actually started in i and ended in j, summed over all i, j. Those weights come straight from the forward-backward of Chapter 6.

Concretely, for one gap the expected i→j jumps = Σ_{start a, end b} P(start=a, end=b | data) · [Van Loan corner for the i→j indicator, restricted to start a, end b]. The forward message at the gap's start and the backward message at its end supply P(start=a, end=b | data). We accumulate these weighted corners over every gap and every patient — that sum is the E-step's output.

This is why learning needs inference: every EM iteration runs a full forward-backward pass (the E-step) to compute the weights, then a single ratio (the M-step) to update Q. The inference you built in Chapters 5 and 6 is not a separate topic — it is the engine inside the E-step. Learning is inference run in a loop, with a parameter update between passes.

The M-Step Ratio, Worked

Suppose after the E-step we have accumulated, across 50 patients, a total expected number of H→D jumps of 38.5 and a total expected time-in-H of 372.0 weeks. The M-step rate update is simply jumps over time: q_HD^new = 38.5 / 372.0 = 0.1035. Starting from a guess of 0.08, this single ratio already moved the estimate toward the true 0.10.

Run the loop again with q = 0.1035: the E-step recomputes the expected jumps and times under the new rate, the ratio updates again, and the estimate creeps closer to 0.10, the data log-likelihood rising every step. EM is guaranteed to never decrease the likelihood — each M-step is a maximization of a lower bound — which is why the simulator's likelihood curve climbs monotonically and never dips.

# A minimal CT-EM loop for the 2-state H/D rate, sketched.
import numpy as np
from scipy.linalg import expm

def vanloan_corner(Q, E, dt):
    N = Q.shape[0]; A = np.zeros((2*N, 2*N))
    A[:N,:N] = Q; A[N:,N:] = Q; A[:N,N:] = E
    return expm(A * dt)[:N, N:]

def em_iter(q, gaps):
    Q = np.array([[-q, q],[0,0]])
    Ej = np.array([[0.0,1.0],[0,0]])    # H->D jump indicator
    Et = np.array([[1.0,0],[0,0]])      # time-in-H indicator
    jumps = sum(vanloan_corner(Q, Ej, dt)[0,1] for dt in gaps)
    time_ = sum(vanloan_corner(Q, Et, dt)[0,0] for dt in gaps)
    return jumps / time_                       # M-step ratio

gaps = [10] * 50                              # 50 patients, one 10w gap each (sketch)
q = 0.08
for it in range(5):
    q = em_iter(q, gaps); print(f"iter {it}: q = {q:.4f}")
# the estimate climbs from 0.08 toward the data-implied rate each iteration

This sketch omits the posterior weighting (it assumes every gap genuinely starts H and ends accounted-for, valid when D is absorbing and the data says so), but it shows the skeleton: build Q from the current guess, Van-Loan each gap for expected jumps and time, take the ratio, repeat. The real implementation folds in the forward-backward weights, but the heart — expm of a block matrix, then a ratio — is exactly this.

Cross-Checking Van Loan Against Closed Forms

Trusting a block-matrix trick requires validation. For the 2-state absorbing chain, both expected statistics have closed forms we can check the corner against. The expected time in H over a gap dt starting in H is (1 − e^−q·dt)/q; the expected number of H→D jumps over the same gap is just the probability the jump happened, 1 − e^−q·dt, since there can be at most one. The table confirms the corner matches both.

# Van Loan corner vs the two closed forms, across gaps.
import numpy as np
from scipy.linalg import expm
q = 0.08; Q = np.array([[-q, q],[0,0]])
Et = np.array([[1.0,0],[0,0]])      # time-in-H indicator
Ej = np.array([[0.0,1.0],[0,0]])      # H->D jump indicator

def corner(E, dt):
    A = np.zeros((4,4)); A[:2,:2] = Q; A[2:,2:] = Q; A[:2,2:] = E
    return expm(A * dt)[:2, 2:]
for dt in (5, 10, 20):
    t_vl = corner(Et, dt)[0,0]; t_cf = (1 - np.exp(-q*dt)) / q
    j_vl = corner(Ej, dt)[0,1]; j_cf = 1 - np.exp(-q*dt)
    print(f"dt={dt}: time {t_vl:.4f}/{t_cf:.4f}  jumps {j_vl:.4f}/{j_cf:.4f}")
# every Van-Loan value matches its closed form to machine precision

Both statistics, at every gap, match their closed forms exactly — the block-matrix trick is not an approximation but an exact computation of the integrals. For a 2-state chain we could skip Van Loan and use the closed forms, but they only exist for trivially simple Q. The whole point of Van Loan is that the same expm-of-a-block-matrix recipe works for any number of stages and any transition structure, where no closed form exists. The 2-state check is just our confidence that the general machine is correct.

Why the Likelihood Climbs Monotonically

EM's signature guarantee is that the data log-likelihood never decreases across iterations. This is not luck — it is structural. Each M-step maximizes a lower bound on the log-likelihood that touches the true log-likelihood at the current parameters, so improving the bound cannot worsen the true objective. In practice this gives a clean diagnostic: if your implementation ever shows a dip in logL, you have a bug, full stop.

# The forward pass gives logL; track it across EM iterations.
import numpy as np
from scipy.linalg import expm

def data_logL(q, obs_seqs, gap_seqs, B, pi):
    Q = np.array([[-q, q, 0],[0,-q/2,q/2],[0,0,0]])
    total = 0.0
    for obs, gaps in zip(obs_seqs, gap_seqs):
        a = pi * B[:, obs[0]]; c = a.sum(); total += np.log(c); a /= c
        for k in range(1, len(obs)):
            a = (a @ expm(Q * gaps[k])) * B[:, obs[k]]
            c = a.sum(); total += np.log(c); a /= c
    return total

# As EM nudges q toward the true rate, data_logL should rise each step.
# A dip would signal a broken E-step or M-step -- a reliable bug detector.

The function reuses the exact forward pass from Chapter 5 — learning literally calls inference. Plotting data_logL at the end of each EM iteration gives the rising curve the simulator draws, and the monotonicity is your contract with the algorithm. When that curve climbs and flattens, EM has converged to a local optimum of the likelihood; the q it returns is the maximum-likelihood rate given the scattered data.

To recap the chapter's logic: learning Q is hard because the jumps are hidden, so we cannot count transitions as discrete Baum-Welch does. CT-EM instead computes the expected continuous-time sufficient statistics — expected jumps and expected time-in-stage — via the Van Loan block-matrix exponential, weights them by the smoothed posteriors, and forms the rate as jumps-over-time in the M-step. Iterate, and the likelihood climbs to a local maximum. The whole machine is one expm of a 2N×2N matrix wrapped in a forward-backward weighting and a ratio — inference, run in a loop, with a parameter update between passes.

The single sentence to remember: a rate is jumps-per-unit-time, so the M-step denominator is expected time-in-stage, not a step count — and that is precisely what discrete Baum-Welch lacks and what Van Loan supplies. Endpoint-counting (treat each gap as one i→j transition) ignores both intermediate jumps and elapsed time, giving a biased Q. Recovering the true continuous-time sufficient statistics is the entire reason the chapter exists, and the matrix exponential is, once again, the tool that makes it tractable.

Notice the beautiful recursion in the architecture: learning contains inference. The E-step runs the full forward-backward of Chapters 5 and 6 to weight each gap's Van Loan corner by the smoothed start/end posteriors; the M-step is then a single ratio. So everything you built for inference is reused, unchanged, inside the learning loop — which is why a CT-HMM library is mostly a forward-backward routine plus an expm-of-a-block-matrix, wrapped in an iteration. Master inference and learning is a short step away.

# CT-EM, distilled to its skeleton:
#   E-step: forward-backward gives smoothed posteriors;
#           Van Loan corner exp([[Q,E],[0,Q]]*dt)[:N,N:] gives expected stats
#   M-step: q_ij = sum(expected i->j jumps) / sum(expected time in i)
#           then set q_ii = -sum_j!=i q_ij  (keep Q a valid generator)
# key: denominator is TIME, not a step count -- that is the CT twist.
# guarantee: logL never decreases; a dip means an implementation bug.
# the E-step LITERALLY calls the Chapter 5/6 forward-backward routines.
# repeat until the data log-likelihood stops rising.

That skeleton — forward-backward for weights, Van Loan for expected stats, a ratio for the update — is the whole of continuous-time learning, and it reuses every inference routine you already built. With learning in hand, only the edge cases (Chapter 10) and the family map (Chapter 11) remain, and you can now fit a CT-HMM end to end.

Chapter 10: Edge Cases — Absorbing, Identifiability, Stiff Q

The clean theory meets messy reality. Three failure modes show up constantly in real CT-HMMs, and knowing them is the difference between a model that ships and one that silently lies. We'll take them one at a time, each with a concrete quantity that visibly breaks.

Absorbing states. Disease models have a "death" or "failure" stage you can't leave — its generator row is all zeros (no outgoing rate). Once the belief mass enters it, exp(Q·dt) keeps it there forever: the column for that stage tends to 1, and observations after absorption add no transition information (there's nothing left to transition to). The stationary distribution is degenerate — all mass on the absorbing stage.

Non-identifiability. Two genuinely different rate matrices can produce nearly identical observation likelihoods under a sparse visit schedule. If you only see a patient at two far-apart times, you constrain a single transition probability over one long gap — one number — and many Q's reproduce that one number. High likelihood does not mean you recovered the true Q.

Stiffness. When rates span orders of magnitude (a fast acute phase at rate 5, a slow chronic phase at rate 0.01), the matrix exponential becomes a numerical hazard for naive methods. A crude discretization like (I + Q·dt/k)^k with small k can produce negative or oscillating "probabilities," while a proper expm stays stable and valid.

The thread through all three: continuous-time models don't fail loudly. An absorbing state quietly stops learning. A non-identifiable Q quietly fits the data while being wrong. A stiff Q quietly returns garbage probabilities. The remedy is always the same posture — check structure and numerics before trusting a fitted model.

The simulator below lets you trigger each mode and watch its signature quantity break: absorbed mass piling up, two likelihood curves nearly coinciding then separating, and divergent transition entries from a bad numerical method.

A Worked Example: Why Sparse Visits Lose Identifiability

Two-state H/D. Model A has q_HD = 0.10. With only two visits 20 weeks apart, the data constrains exactly one quantity: the chance of having moved H→D over that gap. For model A that's P(start H, end D) = 1 − e^−0.10·20 = 1 − e^−2.0 = 1 − 0.1353 = 0.8647.

Now build model B with a different rate structure — say a slightly different q_HD plus a fast reversible H↔intermediate blip — tuned so that its H→D probability over 20 weeks is also ≈ 0.8647. The two far-apart visits see the same single number, 0.8647, so they assign A and B the same data likelihood. The rates are non-identifiable from this schedule.

Add one visit at week 5. Now model A predicts P(still H at week 5) = e^−0.10·5 = e^−0.5 = 0.6065, while model B — with its different early dynamics — predicts something else. The intermediate observation breaks the tie: the two models now disagree on a measurable quantity, and the data can finally tell them apart. Sampling design restores identifiability.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

# (A) absorbing: row of zeros never lets mass leave
Qa = np.array([[-0.1, 0.1], [0, 0]])
print(expm(Qa * 100)[1])           # [0. 1.] — D is a trap forever

# (B) identifiability: same single-gap prob, different rate
print(1 - np.exp(-0.10 * 20))      # 0.8647 over the long gap
print(np.exp(-0.10 * 5))          # 0.6065 at the added mid-visit — separates models

# (C) stiffness: naive discretization vs expm
Qs = np.array([[-5.0, 5.0], [0.01, -0.01]])
dt, k = 1.0, 2
naive = np.linalg.matrix_power(np.eye(2) + Qs * dt / k, k)
print("naive:", naive[0])           # can go negative — invalid probabilities
print("expm: ", expm(Qs * dt)[0])     # stays a valid distribution

The library move is, as always, to use scipy.linalg.expm rather than rolling your own discretization — SciPy's scaling-and-squaring + Pade implementation is robust to stiffness, which the naive (I + Q·dt/k)k is not.

Absorbing States, in Detail

An absorbing stage has an all-zero generator row: no outgoing rate, so once entered it is never left. In a disease model this is "death" or "end-stage failure." The diagnostic signature is that the corresponding column of exp(Q·dt) tends to 1 as dt grows — all mass eventually piles in — and the stationary distribution is the degenerate "all mass on the absorbing stage."

The practical consequence for inference: observations after absorption carry no transition information, because there is nowhere left to transition. They can still inform the emission (confirming you are in the absorbing stage), but they cannot help estimate any rate, because no rate governs a stage with no outgoing transitions. For learning, this means absorbing-stage data constrains only the rates into the absorbing stage, never any rate out of it — there are none.

# Absorbing stage: its column of exp(Q*dt) -> 1; nothing escapes.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])
for dt in (5, 20, 100):
    print(dt, "P(*->stage3):", expm(Q * dt)[:, 2].round(3))
# 5   P(*->stage3): [0.155 0.394 1.   ]
# 20  P(*->stage3): [0.733 0.865 1.   ]
# 100 P(*->stage3): [0.998 0.999 1.   ]  everything ends absorbed

Non-Identifiability, in Detail

Identifiability asks: does the data uniquely determine Q? With sparse irregular visits, often no. Two distant visits constrain exactly one number per stage-pair — the transition probability over that one gap — and many Q's reproduce any single number. High likelihood then certifies "fits the data," not "recovered the truth," and those are very different claims.

The remedy is information, supplied three ways: more or better-timed observations (an intermediate visit constrains a second number, as our worked example showed), priors or regularization on Q (encode prior knowledge of plausible rates), or structural constraints (force progressive-only with no recovery, collapsing the parameter space). Each adds a constraint the sparse data lacked. The honest practice is to report the uncertainty of a fitted Q, not just its point value — a wide posterior on a rate is the model telling you the data could not pin it down.

# Two rates, same single-gap probability -> indistinguishable from one gap.
import numpy as np
# model A: rate 0.10; model B: rate chosen to match over a 20w gap
pA = 1 - np.exp(-0.10 * 20)            # 0.8647
# any qB with 1-exp(-qB*20)=0.8647 gives qB=0.10 too here, but with a
# different model STRUCTURE (e.g. a hidden fast blip) the SAME 0.8647
# can be reproduced by a different effective rate. Add a mid-visit:
mid_A = np.exp(-0.10 * 5)              # 0.6065 -- a SECOND constraint
print("one gap sees:", round(pA, 4))   # 0.8647
print("mid-visit adds:", round(mid_A, 4))  # 0.6065 -- now models must agree on TWO numbers

Stiffness, in Detail

A generator is "stiff" when its rates span orders of magnitude — a fast acute transition at rate 5 alongside a slow chronic one at 0.01. The fast mode forces tiny steps for any explicit method, while the slow mode requires a long horizon, and the ratio is the stiffness. Naive discretization like (I + Q·dt/k)^k produces negative or oscillating "probabilities" unless k is enormous; expm's scaling-and-squaring with a Pade approximant stays valid regardless.

The thread through all three failure modes: continuous-time models fail quietly. An absorbing stage silently stops contributing learning signal; a non-identifiable Q silently fits while being wrong; a stiff Q silently returns invalid probabilities. None throws an error. The defensive posture is to check structure (which stages are absorbing?), check identifiability (is the data dense enough?), and check numerics (use expm, validate that rows sum to 1 and entries are non-negative) before trusting any fitted model.

Stiffness, Quantified

Let us measure how badly the naive method fails on a stiff Q as the slice count varies, so "use expm" is a recommendation backed by numbers rather than a slogan. With rates 5 and 0.01 — a 500× spread — the crude (I + Q·dt/k)^k needs a large k just to stay non-negative, while expm is correct immediately.

# Stiff Q: naive discretization needs many slices to even be valid.
import numpy as np
from scipy.linalg import expm
Qs = np.array([[-5.0, 5.0], [0.01, -0.01]])
exact = expm(Qs * 1.0)
for k in (1, 2, 8, 64):
    naive = np.linalg.matrix_power(np.eye(2) + Qs / k, k)
    ok = (naive >= -1e-9).all() and np.allclose(naive.sum(1), 1)
    print(f"k={k:>2} valid={ok}  err={np.abs(naive-exact).max():.4f}")
print("expm valid=True err=0.0000  (always)")
# k= 1 valid=False err=...   negative entries -- nonsense
# k= 2 valid=False ...
# k=64 valid=True  err small  -- finally usable, at 64x the cost

A Pre-Flight Checklist for a Fitted CT-HMM

Because the failures are silent, a short checklist before trusting any fitted model saves real grief. Run through it every time, the same way a pilot checks instruments before takeoff.

Each check targets one of the failure modes: structure catches accidental absorption, identifiability catches sparse-data ambiguity, numerics catches stiffness, and the monotone-likelihood check catches EM implementation bugs. None is expensive, and skipping them is how a model that "fits beautifully" quietly ships a wrong Q into a clinical decision. The honest CT-HMM practitioner treats these checks as mandatory, not optional — the math is exact, but only if you respect its preconditions.

The Three Remedies for Non-Identifiability, in Code

When the data cannot pin Q, three remedies add the missing constraint. The code below sketches each: denser visits (a second observation per patient), a structural prior (forbid recovery so fewer parameters exist), and regularization (penalize rates far from a prior belief). Each shrinks the set of Q's that explain the data.

import numpy as np

# (1) DENSER VISITS: a mid-visit adds a second constraint.
q = 0.10
long_gap = 1 - np.exp(-q * 20)   # 0.8647  -- one number, many Q fit
mid_gap  = np.exp(-q * 5)        # 0.6065  -- a SECOND number to match
print("constraints:", round(long_gap, 4), round(mid_gap, 4))

# (2) STRUCTURAL PRIOR: progressive-only Q has fewer free rates.
def progressive_Q(rates):     # only forward transitions allowed
    n = len(rates) + 1; Q = np.zeros((n, n))
    for i, r in enumerate(rates): Q[i, i+1] = r; Q[i, i] = -r
    return Q                       # no backward rates to estimate -> fewer params
print(progressive_Q([0.2, 0.1]))

# (3) REGULARIZATION: pull the estimate toward a prior rate q0.
def reg_mstep(jumps, time_, q0, strength):
    return (jumps + strength * q0) / (time_ + strength)   # MAP-style shrinkage
print("regularized q:", round(reg_mstep(38.5, 372.0, 0.1, 50), 4))

Denser visits add empirical constraints; the structural prior removes parameters that the data was never going to identify (you cannot learn a recovery rate from a no-recovery disease); regularization biases the estimate toward prior knowledge when data is thin, trading a little bias for a lot of variance reduction. In practice you combine all three. The unifying principle: identifiability is about whether the data, the structure, and the prior together pin down Q — not the data alone.

And always report uncertainty. A rate learned from two patients with one gap each is barely constrained; the same rate from five hundred patients with dense visits is tight. A point estimate hides that difference; a posterior or bootstrap interval exposes it. The most dangerous CT-HMM is one that reports a confident-looking Q learned from data that could not possibly have constrained it — the checklist and these remedies exist to keep you honest about exactly that.

From-Scratch: A Bootstrap Confidence Interval on a Rate

The honest way to report a learned rate is with an interval, and the bootstrap gives one cheaply: resample patients with replacement, refit q on each resample, and look at the spread. A wide interval is the model admitting the data did not pin the rate down — exactly the identifiability warning made quantitative.

# Bootstrap a confidence interval on a learned rate (sketch).
import numpy as np
rng = np.random.default_rng(0)

# pretend per-patient (jumps, time) sufficient stats from the E-step:
patients = [(0.8, 7.0), (1.0, 9.0), (0.6, 6.5), (0.9, 8.0),
            (0.7, 7.5), (1.1, 9.5), (0.5, 6.0), (0.95, 8.5)]
data = np.array(patients)

def fit_q(rows): return rows[:,0].sum() / rows[:,1].sum()   # M-step ratio

print("point estimate q =", round(fit_q(data), 4))
boot = [fit_q(data[rng.integers(0, len(data), len(data))]) for _ in range(2000)]
lo, hi = np.percentile(boot, [2.5, 97.5])
print(f"95% CI = [{lo:.4f}, {hi:.4f}]")   # width reveals how well data pins q

The point estimate is a single number; the bootstrap interval around it is the honesty. With eight patients the interval is appreciable; with eight hundred it would tighten. Reporting the interval — not just the point — is what separates a defensible CT-HMM fit from one that overclaims. Combined with the structural and numerical checks above, this is the full discipline for trusting (or distrusting) a fitted continuous-time model.

To recap the chapter: three failure modes lurk in real CT-HMMs, and all fail quietly. Absorbing states stop contributing learning signal (their column of exp(Q·dt) tends to 1). Non-identifiability lets many Q's fit sparse data equally well (high likelihood is necessary, not sufficient). Stiffness makes naive discretization return invalid probabilities (always use expm). The remedies — denser or better-timed visits, structural priors, regularization, and reported uncertainty — all add the constraint the data alone lacked. The discipline is: check structure, identifiability, and numerics before trusting any fitted model, because none of these failures throws an error on its own.

# A one-call validity check to run on any exp(Q*dt) before trusting it.
import numpy as np
from scipy.linalg import expm
def valid_transition(P, tol=1e-9):
    rows_ok = np.allclose(P.sum(axis=1), 1, atol=tol)   # each row a distribution
    nonneg  = (P >= -tol).all()                          # no negative probabilities
    return bool(rows_ok and nonneg)

Qgood = np.array([[-0.2,0.2],[0,0]])
Qstiff = np.array([[-5.0,5.0],[0.01,-0.01]])
print("expm good: ", valid_transition(expm(Qgood * 10)))     # True
naive = np.linalg.matrix_power(np.eye(2) + Qstiff, 1)
print("naive stiff:", valid_transition(naive))               # False -- caught!

The check is two lines and catches the most dangerous quiet failure — a transition matrix with negative entries or rows that do not sum to 1, which signals either a broken Q or a numerically unstable method. Run it on every exp(Q·dt) in a pipeline and stiffness bugs surface immediately instead of silently corrupting beliefs. It costs nothing and is the cheapest insurance a CT-HMM implementation can buy.

# Edge-case defenses, distilled:
#   absorbing:        expect column of exp(Q*dt) -> 1; no learning signal out
#   non-identifiable: many Q fit sparse data; add visits / priors / structure
#   stiff:            naive (I+Q*dt/k)^k goes negative; ALWAYS use expm
# validate: rows of exp(Q*dt) sum to 1, no negatives; logL climbs in EM;
#           report a bootstrap interval on any learned rate.
# all three fail SILENTLY -- no exception, just wrong/biased answers.
# remedies stack: more/better-timed visits + priors + structure.
# verification is mandatory, not optional, for a CT-HMM you trust.

The unifying habit: continuous-time models fail silently, so verification is not optional. Run the structural, identifiability, and numerical checks every time, and report intervals rather than point estimates — that discipline is what turns a model that merely fits into one you can actually trust.

With the failure modes understood and defended against, you have not just the algorithms but the judgment to deploy them responsibly. The final chapter steps back to place the CT-HMM in the broader family of latent-state filters, where every piece you built reappears as one costume of a single underlying idea.

The recurring lesson across all three edge cases is humility about what data can tell you. A model is only as trustworthy as the data that constrained it and the checks you ran on it; a confident-looking Q from thin, sparse, or numerically mishandled data is the most dangerous output a CT-HMM can produce. Treat verification as part of the model, not an afterthought.

None of these failures is exotic; all three appear in routine clinical and reliability datasets. Internalize the checklist now and they become reflexes rather than surprises.

Chapter 11: Connections — Where CT-HMM Lives

Step back and place the CT-HMM on the map. It is not an exotic one-off — it is the continuous-time member of the latent-state filtering family, sitting in a clean grid alongside models you've likely already met. Seeing the grid is how you stop memorizing algorithms and start recognizing one idea wearing different clothes.

The family organizes along two axes: is the hidden state discrete or continuous, and is time discrete or continuous? Four cells. Discrete state, discrete time is the discrete HMM. Discrete state, continuous time is the CT-HMM (this lesson). Continuous state, discrete time is the Kalman filter. Continuous state, continuous time is Kalman-Bucy / SDE filtering.

Every cell uses the same predict-update rhythm; only the transition operator changes. The discrete HMM uses a fixed matrix A. The CT-HMM uses exp(Q·dt), recomputed per gap. The Kalman filter uses a linear matrix F. They are siblings, not strangers — the Bayes filter is the abstract parent of all of them.

The CT-HMM is the discrete HMM with the clock set free. And it is the discrete-state cousin of the Kalman filter: where Kalman tracks a continuous position with Gaussian uncertainty, the CT-HMM tracks a discrete stage with a categorical belief. They even share the learning engine — the EM machinery of Chapter 9 is the same EM that fits a Gaussian mixture.

In the wild, this is the workhorse of disease-progression modeling (the kidney example is real clinical practice), equipment-degradation and reliability, and customer-lifecycle analysis — anywhere a hidden condition decays on its own schedule and you only get to peek at irregular moments. The continuous-time framing isn't academic; it's what makes irregular real-world data tractable.

The single most important relationship to internalize: the discrete HMM is a special case, recovered when observations happen to be evenly spaced. The worked example below proves it numerically — sample a CT-HMM on a uniform grid and you get exactly a discrete HMM. Free the grid and you get the general CT-HMM. This is the closing arc of the whole lesson.

A Worked Example: The Discrete HMM Falls Out

Take a 2-state CT-HMM with Q = [[−0.2, 0.2], [0, 0]] (Healthy leaks to absorbing Diseased at rate 0.2). Discretize on a uniform 1-week grid: the per-step transition matrix is A = exp(Q·1) = [[0.8187, 0.1813], [0, 1]]. That is now a fixed discrete-HMM transition matrix.

Run the discrete HMM with this A for 5 steps. The (Healthy→Healthy) entry of A⁵ is 0.8187⁵ = 0.3679.

Now ask the CT-HMM directly for the 5-week transition: exp(Q·5), whose (Healthy→Healthy) entry is e^−0.2·5 = e^−1.0 = 0.3679. Identical. The discrete HMM running 5 weekly steps and the CT-HMM evaluating one 5-week gap give the exact same answer — because A⁵ = (exp(Q·1))⁵ = exp(Q·5).

That equality is the whole relationship. The discrete HMM is precisely a CT-HMM sampled on a uniform grid; the CT-HMM is what you get when you free the grid to arbitrary real-valued gaps. There is no approximation involved — uniform sampling makes the per-gap exp(Q·dt) collapse into one repeated matrix, and you're back in discrete-HMM land.

From-Scratch and Library

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2, 0.2], [0, 0]])
A = expm(Q * 1.0)                       # discretize on a weekly grid
print(A)                                  # [[0.8187 0.1813] [0. 1.]]

lhs = np.linalg.matrix_power(A, 5)[0, 0]   # discrete HMM, 5 weekly steps
rhs = expm(Q * 5)[0, 0]                # CT-HMM, one 5-week gap
print(lhs, rhs)                          # 0.367879  0.367879
assert np.isclose(lhs, rhs)              # discrete HMM == uniformly-sampled CT-HMM

# Related lessons (the latent-state filtering family):
#   ctmc.html         — the hidden engine: Q, holding times, exp(Qt)
#   hmm.html          — the discrete-time cousin (fixed A)
#   kalman-filter.html— continuous state, discrete time
#   bayes-filter.html — the abstract predict-update parent
#   em-gmm.html       — the shared EM learning engine

The library move is the one-line np.isclose(matrix_power(A,5)[0,0], expm(Q*5)[0,0]) — a numerical proof that the two worldviews coincide on a uniform grid. That single assertion is the bridge between everything you knew before this lesson and everything you learned in it.

The Four Cells, One at a Time

The latent-state filtering family is a 2×2 grid: state (discrete or continuous) crossed with time (discrete or continuous). Walking the cells makes the CT-HMM's place precise. Discrete state, discrete time is the discrete HMM: a fixed transition matrix A, applied once per tick. The casino-dice and part-of-speech taggers live here.

Discrete state, continuous time is the CT-HMM of this lesson: the transition is exp(Q·dt), recomputed per real gap. Disease staging, equipment degradation, and customer lifecycles live here — anywhere a categorical hidden condition drifts on its own clock and is seen at irregular moments.

Continuous state, discrete time is the Kalman filter: the hidden state is a real-valued vector with Gaussian uncertainty, transitioned by a linear matrix F at fixed steps. Robot tracking and sensor fusion live here. Continuous state, continuous time is Kalman-Bucy / SDE filtering: a stochastic differential equation drives the state, observed continuously or at irregular times.

Every cell runs the identical predict-update rhythm; only the transition operator and the representation of uncertainty change. That is the unifying insight: you are not memorizing four algorithms but recognizing one idea — recursive Bayesian filtering — in four costumes. The Bayes filter is the abstract parent that all four specialize.

The Shared Learning Engine

The kinship runs deeper than inference: these models share a learning engine too. The CT-EM of Chapter 9 is the same expectation-maximization that fits a Gaussian mixture and the same Baum-Welch that fits a discrete HMM. In every case the E-step infers the hidden variables under the current parameters, and the M-step re-estimates the parameters from those inferred quantities.

What changes per model is only the sufficient statistics. A Gaussian mixture needs expected responsibilities; a discrete HMM needs expected transition counts; a CT-HMM needs expected transition jump-mass and time-in-stage (the continuous-time analogue of counts, recovered by Van Loan). The skeleton — infer hidden, re-estimate, repeat, likelihood climbs — is identical. Learn EM once and you have learned it for the whole family.

Where CT-HMMs Earn Their Keep

This is not an academic curiosity. CT-HMMs are the workhorse of disease-progression modeling — the kidney-staging example is real clinical practice, with eGFR readings at irregular visits driving stage estimates and treatment decisions. They model equipment degradation (a machine drifts through health states, inspected sporadically) and customer lifecycles (a user moves through engagement states, observed only when they act).

The common thread is irregular observation of a hidden condition that decays on its own schedule. Wherever your data arrives at uneven times and the underlying state keeps moving between observations, the CT-HMM is the honest model — not a discrete HMM with phantom steps, not a last-reading heuristic, but a model that treats the gap as the signal it is. The continuous-time framing is precisely what makes messy real-world timing tractable.

A Second Limiting Case: Uniform Sampling Recovers the Discrete HMM

The worked example above proved A⁵ = exp(Q·5) for one Q. Let us make the "discrete HMM is a uniformly-sampled CT-HMM" claim fully explicit with a second angle that the build spec singles out. Take the absorbing Q = [[−0.2, 0.2], [0, 0]] and form A = exp(Q·1) = [[0.8187, 0.1813], [0, 1]] — a fixed weekly transition matrix.

Now A⁵[0,0] = 0.8187⁵ = 0.3679, which equals exp(Q·5)[0,0] = e^−1.0 = 0.3679. Running the discrete HMM with this A for five weekly steps gives bit-for-bit the same answer as the CT-HMM evaluating one five-week gap. The discrete machine is not an approximation of the continuous one here — on a uniform grid it is exactly the continuous one, because A^k = (exp(Q·1))^k = exp(Q·k).

# Discrete HMM on a uniform grid == CT-HMM, exactly.
import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2, 0.2], [0, 0]])
A = expm(Q * 1.0)                    # weekly discrete-HMM matrix
print("A =", A.round(4))                  # [[0.8187 0.1813] [0. 1.]]
for weeks in (5, 10):
    disc = np.linalg.matrix_power(A, weeks)[0,0]
    cont = expm(Q * weeks)[0,0]
    print(f"{weeks}w: discrete={disc:.6f}  continuous={cont:.6f}")
# 5w:  discrete=0.367879  continuous=0.367879
# 10w: discrete=0.135335  continuous=0.135335
assert np.allclose(np.linalg.matrix_power(A, 5), expm(Q * 5))

The discrete and continuous numbers match to six digits at every horizon. So the discrete HMM is the special case "all gaps equal to the grid spacing," and the CT-HMM is the general case "gaps are arbitrary real numbers." Freeing the grid is the entire generalization — no new probability theory, just permission for dt to vary. That is the closing arc: everything you knew about discrete HMMs is still true, on a uniform grid; the CT-HMM simply lifts the uniform-grid restriction.

What Changes and What Stays, Across the Family

Reading the table column by column, the predict-update skeleton and the EM learning loop are constant; only the state representation and the transition operator vary. That invariance is the payoff of seeing the family as one: master the skeleton once, and each model is a small substitution away. The CT-HMM's substitution — a fixed matrix becomes a matrix-valued function of the gap — is the one this whole lesson built.

The Whole Pipeline, in One Program

To close, here is everything the lesson built, assembled into one runnable pipeline: generate irregular data from a hidden Q, filter it, smooth it, decode it, and (sketched) learn Q back. Each call is a chapter; seeing them composed is the point of the connections chapter.

import numpy as np
from scipy.linalg import expm

Q = np.array([[-0.2,0.2,0],[0,-0.1,0.1],[0,0,0]])   # Ch 1
B = np.array([[0.7,0.25,0.05],[0.2,0.6,0.2],[0.05,0.3,0.65]])    # Ch 3
pi = np.array([0.8,0.15,0.05]); obs = [0,1,2]; gaps = [None,5,2]

def forward():                                # Ch 5
    a = pi * B[:,obs[0]]; A = [a/a.sum()]; logL = np.log(a.sum())
    for k in range(1, len(obs)):
        a = (A[-1] @ expm(Q*gaps[k])) * B[:,obs[k]]
        logL += np.log(a.sum()); A.append(a/a.sum())
    return A, logL

def backward():                               # Ch 6
    be = [None]*len(obs); be[-1] = np.ones(3)
    for k in range(len(obs)-2, -1, -1):
        be[k] = expm(Q*gaps[k+1]) @ (B[:,obs[k+1]] * be[k+1])
    return be

A, logL = forward(); be = backward()
gamma = [(a*b)/(a*b).sum() for a,b in zip(A, be)]   # Ch 6 smoothing
print("filtered last:", A[-1].round(3))
print("smoothed v0: ", gamma[0].round(3))
print("data logL:   ", round(logL, 4))            # -2.6063, the verified checkpoint

One file, every algorithm, one shared expm. The forward pass produces the filtered beliefs and the −2.6063 log-likelihood; the backward pass and product produce the smoothed beliefs; Viterbi (Chapter 7) and CT-EM (Chapter 9) bolt on with the same per-gap matrix exponential. This is the concrete sense in which the CT-HMM is "the HMM with the clock set free": every routine you know, with A swapped for exp(Q·dt), composed into a working pipeline.

Where to Go Next

If this lesson landed, three directions extend it naturally. Continuous emissions (Gaussian per stage, sketched in Chapter 3) handle real-valued sensors like eGFR directly. Hidden semi-Markov models relax the exponential-holding-time assumption when sojourn times are not memoryless. And full Bayesian treatment of Q — a posterior over rates rather than a point estimate — addresses the identifiability concerns of Chapter 10 head-on.

But the core you now own — one rate matrix Q, the matrix exponential exp(Q·dt) evaluated per real gap, and the predict-update-decode-learn routines built on it — is the load-bearing structure under all of them. You can read a clinical-staging CT-HMM paper, implement the forward-backward from scratch, and reason about when a fitted Q is trustworthy. That is the whole skill, assembled from one object and one operation.

One Last Numeric: The Family in Three Lines

To leave you with the unity in code, here are the three discrete-state-or-discrete-time family members' transition steps side by side — the same predict, three operators. Run it and see that the discrete HMM and the uniformly-sampled CT-HMM produce identical steps, while the CT-HMM generalizes to any gap.

import numpy as np
from scipy.linalg import expm
Q = np.array([[-0.2,0.2],[0,0]])
b = np.array([1.0, 0])

A = expm(Q * 1.0)                         # discrete HMM weekly matrix
print("discrete HMM, 1 step:  ", (b @ A).round(4))
print("CT-HMM, 1-week gap:    ", (b @ expm(Q * 1.0)).round(4))  # identical
print("CT-HMM, 3.7-week gap:  ", (b @ expm(Q * 3.7)).round(4))  # a gap A^k cannot do
# discrete HMM, 1 step:   [0.8187 0.1813]
# CT-HMM, 1-week gap:     [0.8187 0.1813]   same as discrete on a uniform grid
# CT-HMM, 3.7-week gap:   [0.4771 0.5229]   a fractional gap no integer power gives

The first two lines are bit-identical — on a uniform grid the CT-HMM is the discrete HMM. The third line is the generalization: a 3.7-week gap, which no integer power A^k can express, but exp(Q·3.7) handles exactly. That fractional gap is the entire reason the CT-HMM exists, and it is the one capability that lifts the discrete HMM into the continuous-time world. Everything else — emissions, filtering, smoothing, decoding, learning — you carried over unchanged.

So the map closes where it opened. Chapter 0 argued that real data arrives at irregular times and a grid is the wrong tool; Chapter 11 proves the grid is merely the special case of equal gaps, and the CT-HMM is the honest general model. One rate matrix, one matrix exponential per gap, and the whole family of latent-state filters opens up. That is a durable mental model — carry it forward to the Kalman filter, the HMM, and beyond.

# The complete CT-HMM toolkit, named -- every routine the lesson built.
#   predict(b, dt)      = b @ expm(Q*dt)            # Ch 4: the gap enters
#   update(b, obs)      = normalize(b * B[:,obs])   # Ch 3/5: a reading sharpens
#   forward(obs, gaps)  = alternate predict+update  # Ch 5: filtering + logL
#   backward(obs, gaps) = right-to-left beta        # Ch 6: future evidence
#   smooth              = normalize(alpha * beta)   # Ch 6: all-data posterior
#   interp(g, tau)      = g @ expm(Q*tau)           # Ch 6: between-visit belief
#   viterbi(obs, gaps)  = log-space max + backptr   # Ch 7: best path
#   vanloan_corner      = expm([[Q,E],[0,Q]]*dt)    # Ch 9: expected stats
# Every line above shares ONE object (Q) and ONE operation (expm).

That single block is the lesson, indexed. Eight routines, each a chapter, all resting on the rate matrix Q and the matrix exponential. If you can read this list and recall what each line does and why, you have internalized the CT-HMM — not as a collection of algorithms, but as one idea (a hidden chain on its own clock) expressed through one operation (exp(Q·dt)) at every irregular gap. "What I cannot create, I do not understand" — and now you can create all of it.

# One predict step, three family members -- same shape, different operator.
import numpy as np
from scipy.linalg import expm

def hmm_predict(belief, A):          # discrete HMM: fixed matrix
    return belief @ A
def cthmm_predict(belief, Q, dt):    # CT-HMM: per-gap exponential
    return belief @ expm(Q * dt)
def kalman_predict(x, F):           # Kalman: linear map (mean part)
    return F @ x

Q = np.array([[-0.2,0.2],[0,0]]); A = expm(Q)   # uniform-grid A from Q
b = np.array([1.0, 0])
print("HMM   :", hmm_predict(b, A).round(4))           # [0.8187 0.1813]
print("CT-HMM:", cthmm_predict(b, Q, 1.0).round(4))    # [0.8187 0.1813] (same!)
print("CT-HMM:", cthmm_predict(b, Q, 3.7).round(4))    # [0.4771 0.5229] (any gap)

Three functions, one shape: a belief goes in, a propagated belief comes out, and the only difference is the operator — a fixed A, a per-gap exp(Q·dt), or a linear F. The HMM and the uniform-grid CT-HMM print identical results because, on equal spacing, they are the same model. Swap the operator and you move between cells of the family; keep the predict-update skeleton and you reuse everything you know. That is the unifying lesson to carry forward.

And the third line — a 3.7-week gap — is the whole reason the CT-HMM earns its own name: no integer power of a weekly A can express a fractional gap, but exp(Q·3.7) does it exactly. That single capability, evaluated once per irregular gap, is the entire generalization from discrete to continuous time. Everything else in the family — the predict-update rhythm, the EM learning loop, the belief representation — you carried over unchanged. Master the skeleton, swap the operator, and you have mastered the family.

This is why seeing the family as one unit is worth more than memorizing four algorithms. The discrete HMM, the CT-HMM, the Kalman filter, and Kalman-Bucy filtering are not four things to learn but one thing — recursive Bayesian filtering — wearing four costumes that differ only in the transition operator and the belief representation. Learn the predict-update-decode-learn skeleton once, and each model becomes a one-line substitution. The CT-HMM's substitution, the per-gap matrix exponential, is the one this lesson built from scratch.

So the lesson ends where it began, with the gap. Chapter 0 insisted that the time between observations is signal, not noise to discretize away. Eleven chapters later, that single commitment has produced a complete inference-and-learning toolkit — one rate matrix, one matrix exponential per gap, and every routine in the latent-state family within reach. Carry the gap forward as a first-class quantity, and irregular real-world data stops being a problem and becomes just another input the math handles exactly.

If you want to test your grip, try the Feynman move: close this page and re-derive the predict step, explain why exp(Q·dt) marginalizes over hidden jumps, and sketch why discrete Baum-Welch's counts fail in continuous time. If you can do those three from memory, you own the CT-HMM — not as memorized formulas, but as a small set of ideas you could rebuild from scratch. That is the bar this lesson set out to clear, and it is the bar worth holding yourself to before moving on.

And if you cannot quite do all three yet, the Teach Mode on this lesson is built for exactly that: explaining each idea aloud, at a whiteboard, is the fastest way to find the seams in your understanding and close them. The model is small enough to hold entirely in your head — one matrix, one operation — which is precisely what makes it worth mastering completely.

That compactness is the gift of the continuous-time framing: a single rate matrix Q and a single operation, the matrix exponential, generate the entire toolkit. Hold those two and you hold the lesson.

# The family, distilled: same predict-update, swap the operator.
#   discrete HMM : transition = fixed A          | belief = categorical
#   CT-HMM       : transition = expm(Q*dt) / gap | belief = categorical
#   Kalman       : transition = linear F         | belief = Gaussian
#   Kalman-Bucy  : transition = SDE              | belief = Gaussian
# learning: discrete=Baum-Welch counts; CT-HMM=Van Loan; all are EM.
# discrete HMM == CT-HMM on a uniform grid (A = expm(Q*Delta)).
# abstract parent of all four: the Bayes filter (predict-update).
# CT-HMM applications: disease staging, equipment wear, customer lifecycle.
# the through-line: treat the GAP as signal, never noise to discretize.
# master the predict-update skeleton once; each model swaps one operator.
# next steps: continuous (Gaussian) emissions; semi-Markov holds; Bayesian Q.