Probabilistic Reasoning & Bayesian Networks

Before computing anything, we must agree on what probability means. There are two main interpretations:

Frequentist: Probability is the long-run frequency of an event in repeated trials. P(heads) = 0.5 means: if you flip a fair coin infinitely many times, half will be heads. Clean, but useless for one-off events like "what is the probability it will rain tomorrow?"

Bayesian (Subjective): Probability is a degree of belief — a number in [0,1] representing how confident a rational agent is in a proposition. P(rain tomorrow) = 0.7 means you would accept a bet that pays $1 if it rains, at a cost of up to $0.70. Cox's theorem proves this is the only coherent way to assign degrees of belief: any other system either violates transitivity or allows Dutch books (bets you lose regardless of outcome).

The axioms (Kolmogorov):

1. P(A) ≥ 0 for all events A
2. P(Ω) = 1 (the sample space has probability 1)
3. P(A ∪ B) = P(A) + P(B) if A ∩ B = ∅

From these three axioms, everything else follows:

• P(¬A) = 1 − P(A)
• P(A ∪ B) = P(A) + P(B) − P(A ∩ B) (inclusion-exclusion)
• P(A, B) = P(A | B) · P(B) ← the product rule, crucial everywhere

Joint distributions: For n binary random variables, the full joint distribution is a table of 2^n numbers summing to 1. For n=30, that is over a billion entries — intractable to store or reason over directly. This exponential curse is precisely why Bayesian networks were invented: to factor the joint into local conditionals.

The marginalization identity: P(A) = Σ_b P(A, B=b). You sum out the variable you don't care about. This is called the sum rule or law of total probability.

from itertools import product

# Full joint over 3 binary variables: Cloudy(C), Sprinkler(S), Rain(R)
# Table indexed by (c, s, r)
joint = {
    (0,0,0): 0.192, (0,0,1): 0.048,
    (0,1,0): 0.144, (0,1,1): 0.216,
    (1,0,0): 0.288, (1,0,1): 0.072,
    (1,1,0): 0.040, (1,1,1): 0.000,
}
assert abs(sum(joint.values()) - 1.0) < 1e-9

# Marginalise out C and R to get P(S)
P_S1 = sum(joint[(c,1,r)] for c in (0,1) for r in (0,1))
print(f'P(Sprinkler=True) = {P_S1:.3f}')  # 0.400

# Conditional: P(R=1 | S=1)
P_S1_and_R1 = sum(joint[(c,1,1)] for c in (0,1))
P_R1_given_S1 = P_S1_and_R1 / P_S1
print(f'P(Rain=T | Sprinkler=T) = {P_R1_given_S1:.3f}')  # 0.540

Bayes' theorem is the single most important formula in probabilistic reasoning:

$$P(H \mid E) = \frac{P(E \mid H) \cdot P(H)}{P(E)}$$

In English: posterior = likelihood × prior / evidence.

Each term has a name and a role:

• Prior P(H): your belief in hypothesis H before seeing evidence E
• Likelihood P(E | H): how probable the observed evidence is if H is true — this is the model's voice
• Posterior P(H | E): your updated belief after seeing E — this is the output
• Evidence P(E): a normalizing constant, equal to Σ_h P(E|H=h)·P(H=h) — ensures the posterior sums to 1

Worked example — medical diagnosis: Disease D has prevalence 0.001. A test has sensitivity P(T⁺|D) = 0.99 and specificity P(T⁻|¬D) = 0.99 (false positive rate = 0.01).

P(D | T⁺) = P(T⁺|D)·P(D) / P(T⁺)

P(T⁺) = P(T⁺|D)·P(D) + P(T⁺|¬D)·P(¬D) = 0.99×0.001 + 0.01×0.999 = 0.01098

P(D | T⁺) = (0.99 × 0.001) / 0.01098 ≈ 0.0902

Only 9% chance of disease even with a positive test! This counterintuitive result — called the base-rate fallacy — is why medical screening of low-prevalence conditions needs confirmatory testing. The rare prior dominates.

Sequential (online) Bayesian updating: Each posterior becomes the next prior. After seeing evidence E1, E2, ..., En (assuming they are conditionally independent given H):

P(H | E1..En) ∝ P(H) · Π P(Ei | H)

This is how a spam filter accumulates evidence word-by-word.

The normalizing constant trick: Often you don't need P(E) explicitly. Compute the unnormalized posteriors for each hypothesis and then divide by their sum:

P(Hi | E) = α · P(E | Hi) · P(Hi) where α = 1 / Σ_j P(E|Hj)·P(Hj)

def bayes_update(prior: dict, likelihoods: dict, evidence: str) -> dict:
    """
    prior: {hypothesis: probability}
    likelihoods: {hypothesis: {evidence_value: probability}}
    evidence: the observed evidence value
    Returns updated posterior.
    """
    unnorm = {h: prior[h] * likelihoods[h][evidence] for h in prior}
    total = sum(unnorm.values())
    return {h: v / total for h, v in unnorm.items()}

# Sequential spam filter
prior = {'spam': 0.3, 'ham': 0.7}
# P(word='free' | class)
word_likelihoods = {
    'spam': {'free': 0.60, 'meeting': 0.05},
    'ham':  {'free': 0.05, 'meeting': 0.30},
}
after_free = bayes_update(prior, word_likelihoods, 'free')
print(after_free)  # {'spam': ~0.837, 'ham': ~0.163}
after_both = bayes_update(after_free, word_likelihoods, 'meeting')
print(after_both)  # spam probability drops after 'meeting'

Independence: A and B are independent (A ⊥ B) iff P(A, B) = P(A)·P(B), equivalently P(A|B) = P(A). Knowing B tells you nothing about A.

Conditional independence: A and B are conditionally independent given C (A ⊥ B | C) iff:

P(A, B | C) = P(A | C) · P(B | C)

This is much more nuanced than marginal independence. Two variables can be marginally dependent but conditionally independent, or vice versa.

The toothache-catch-cavity example (from AIMA): Toothache (T) and Probe-Catch (C) are both caused by Cavity (V). Once you know whether there's a cavity, T and C are conditionally independent — the cavity explains the correlation between them. But marginally, T and C are correlated (both more likely when cavity is present).

P(T, C | V) = P(T | V) · P(C | V) ✓ conditional independence P(T, C) ≠ P(T) · P(C) ✗ NOT marginally independent

The chain rule (product rule generalized): P(X1, X2, ..., Xn) = P(X1) · P(X2|X1) · P(X3|X1,X2) · ... · P(Xn|X1,...,Xn-1)

This is exact and always valid. The trick is that conditional independence lets us drop variables from the conditioning sets:

P(Sprinkler | Cloudy, Rain) = P(Sprinkler | Cloudy) ← if Rain doesn't cause Sprinkler

So instead of conditioning on all predecessors in the chain rule, we condition only on the Markov blanket. This is the entire theoretical basis of Bayesian networks: we encode which variables are relevant parents for each node, replacing a full joint with a product of local CPTs.

Naïve Bayes is the simplest application: assume all features Xi are conditionally independent given the class Y:

P(Y | X1,...,Xn) ∝ P(Y) · Π P(Xi | Y)

This reduces an exponential table to n linear tables. Despite the obviously wrong independence assumption, Naïve Bayes classifiers work remarkably well on text classification — the decision boundary often correct even when probabilities are miscalibrated.

import math
from collections import defaultdict

class NaiveBayes:
    def fit(self, X, y):
        n = len(y)
        self.classes = set(y)
        self.log_prior = {}
        self.log_lik = {}  # log P(feature | class)
        for c in self.classes:
            idx = [i for i,yi in enumerate(y) if yi == c]
            self.log_prior[c] = math.log(len(idx) / n)
            # Laplace smoothing: add 1 to each count
            counts = defaultdict(int)
            for i in idx:
                for feat, val in X[i].items():
                    counts[(feat, val)] += 1
            vocab = {feat for row in X for feat in row}
            self.log_lik[c] = {
                (f, v): math.log((counts[(f,v)] + 1) / (len(idx) + 2))
                for f in vocab for v in (True, False)
            }
    def predict(self, x):
        scores = {}
        for c in self.classes:
            s = self.log_prior[c]
            for feat, val in x.items():
                s += self.log_lik[c].get((feat, val), -10)
            scores[c] = s
        return max(scores, key=scores.get)

A Bayesian network (Bayes net, belief network, DAG model) is a pair (G, Θ) where:

• G is a directed acyclic graph (DAG) with one node per random variable
• Θ is a set of conditional probability tables (CPTs): one table per node, giving P(Xi | Parents(Xi))

The key semantic claim: each node Xi is conditionally independent of its non-descendants given its parents.

Given this, the joint factorizes as:

P(X1, ..., Xn) = Π_{i=1}^{n} P(Xi | Parents(Xi))

This is the Bayes net factorization. It is not just a computational trick — it is a modeling claim about the generative process of the world.

How to construct a Bayes net:

1. List the relevant variables
2. Order them (roughly causal: causes before effects)
3. For each variable, select its minimal set of direct causes as parents
4. Fill in the CPT: P(Xi | each combination of parent values)

Example — the classic Burglary network (AIMA Ch. 13): Variables: Burglary (B), Earthquake (E), Alarm (A), JohnCalls (J), MaryCalls (M) Edges: B→A, E→A, A→J, A→M

Joint: P(B,E,A,J,M) = P(B)·P(E)·P(A|B,E)·P(J|A)·P(M|A)

Storage savings: Full joint = 2^5 = 32 numbers. Bayes net CPTs: 1+1+4+2+2 = 10 numbers (minus constraints).

For n variables each with k parents out of n-1: CPT size is 2^k per node vs 2^n for the full joint. With k=3, a 30-variable network needs 30×8 = 240 numbers instead of 10^9.

Building a CPT: For Alarm given Burglary and Earthquake:

These come from domain experts, data, or both. The network is a model — its validity depends on whether the independence assumptions are warranted.

# Represent the Burglary Bayes net
# Using pgmpy (pip install pgmpy)
from pgmpy.models import BayesianNetwork
from pgmpy.factors.discrete import TabularCPD
from pgmpy.inference import VariableElimination

model = BayesianNetwork([('B','A'),('E','A'),('A','J'),('A','M')])

cpd_B = TabularCPD('B', 2, [[0.999],[0.001]])
cpd_E = TabularCPD('E', 2, [[0.998],[0.002]])
cpd_A = TabularCPD('A', 2,
    [[0.999, 0.06, 0.71, 0.05],
     [0.001, 0.94, 0.29, 0.95]],
    evidence=['B','E'], evidence_card=[2,2])
cpd_J = TabularCPD('J', 2, [[0.95,0.10],[0.05,0.90]],
    evidence=['A'], evidence_card=[2])
cpd_M = TabularCPD('M', 2, [[0.99,0.30],[0.01,0.70]],
    evidence=['A'], evidence_card=[2])

model.add_cpds(cpd_B, cpd_E, cpd_A, cpd_J, cpd_M)
assert model.check_model()

infer = VariableElimination(model)
# P(Burglary | JohnCalls=True, MaryCalls=True)
q = infer.query(['B'], evidence={'J':1, 'M':1})
print(q)  # P(B=True | J,M) ≈ 0.284

The power of Bayes nets is that you can read off conditional independence relationships directly from the graph structure, without touching any numbers. This is done via d-separation (directional separation).

Consider a path between nodes X and Y through a sequence of nodes. The path is blocked given evidence set E if any node Z along the path satisfies one of these patterns:

Three fundamental connection patterns (for a node Z on a path X → Z → Y or X ← Z ← Y or X ← Z → Y or X → Z ← Y):

1. Chain (serial): X → Z → Y. Blocked if Z is in E. (Z is the "relay": once you know Z, X tells you nothing more about Y.)

2. Fork (common cause): X ← Z → Y. Blocked if Z is in E. (Z is the common cause: conditioning on Z removes the X-Y correlation.)

3. Collider (v-structure): X → Z ← Y. Blocked if Z and all its descendants are NOT in E. Unblocked when Z (or a descendant) IS in E — this is called Berkson's paradox or explaining away.

Colliders are the counterintuitive case: two independent causes of Z become dependent once you observe Z. Example: if Alarm = True, then observing Burglary = False makes Earthquake more probable — the earthquake "explains away" the alarm.

D-separation definition: X and Y are d-separated by E if every path between X and Y is blocked. If d-separated, then X ⊥ Y | E in every distribution consistent with the graph.

Markov blanket: The minimal set that renders a node Xi conditionally independent of all others: MB(Xi) = Parents(Xi) ∪ Children(Xi) ∪ Co-parents(Xi). For inference, you only need the Markov blanket.

Why this matters: D-separation lets you answer questions like "does observing symptoms A and B make disease C and D independent?" purely by graph inspection — no arithmetic needed. This is what makes Bayes nets an expressive modeling language in addition to a computational tool.

# Demonstrating explaining away (collider effect) numerically
import numpy as np

# Burglary and Earthquake are independent priors
p_B = 0.001
p_E = 0.002

# P(Alarm=T | B, E)
def p_alarm(b, e):
    table = {(1,1):0.95, (1,0):0.94, (0,1):0.29, (0,0):0.001}
    return table[(b,e)]

# P(B=1 | Alarm=True) vs P(B=1 | Alarm=True, Earthquake=False)
# Marginal: ignoring E
def p_B_given_alarm():
    num = sum(p_alarm(1,e)*p_B*(p_E if e else 1-p_E) for e in [0,1])
    den = sum(p_alarm(b,e)*(p_B if b else 1-p_B)*(p_E if e else 1-p_E)
              for b in [0,1] for e in [0,1])
    return num/den

def p_B_given_alarm_and_no_E():
    # Earthquake ruled out — burglar probability goes UP
    num = p_alarm(1,0)*p_B*(1-p_E)
    den = sum(p_alarm(b,0)*(p_B if b else 1-p_B)*(1-p_E) for b in [0,1])
    return num/den

print(f"P(B|Alarm)         = {p_B_given_alarm():.4f}")   # ~0.376
print(f"P(B|Alarm, ¬Eq)    = {p_B_given_alarm_and_no_E():.4f}")  # ~0.657
# Earthquake=False "explains away" less, so Burglary becomes more likely

Inference in a Bayes net means answering queries of the form: given observed evidence E=e, compute the posterior P(Q | E=e) for query variable(s) Q.

The naive approach — sum over the full joint — is exponential. Variable elimination (VE) exploits the factored structure to push summations inside products, drastically reducing cost.

The key operation: Summing out (marginalizing) a variable Z from a product of factors:

Σ_z f1(Z, ...) · f2(Z, ...) · f3(...) = [Σ_z f1(Z,...) · f2(Z,...)] · f3(...)

You only include factors that contain Z inside the sum. The result is a new factor over the remaining variables.

Variable Elimination Algorithm:

Given: Query Q, evidence E=e, hidden variables H = {H1,...,Hk}
1. Initialize factors = all CPTs, with evidence variables clamped to observed values
2. For each hidden variable Hi (in some elimination order):
   a. Collect all factors containing Hi
   b. Multiply them together into a product factor
   c. Sum out Hi from the product → new factor (no longer contains Hi)
   d. Replace collected factors with the new factor
3. Multiply all remaining factors
4. Normalize to get P(Q | e)

Complexity: The cost is exponential in the treewidth of the graph (the maximum factor size encountered during elimination). For tree-structured networks, VE is linear. For graphs with many loops, treewidth can be large.

Elimination ordering matters: finding the optimal ordering is NP-hard, but greedy heuristics (min-degree, min-fill) work well in practice.

Junction tree algorithm: A more systematic exact inference algorithm that preprocesses the graph into a junction tree (clique tree), enabling a single message-passing pass to answer all queries. This is what underlies efficient implementations in production Bayes net toolkits.

The key insight from VE: you never materialize the full joint. Instead you stream through factors, summing out variables as soon as all factors containing them have been collected.

import numpy as np
from functools import reduce

# Factor represented as (variables_tuple, values_array)
# Simple VE for P(B | J=True, M=True) in Burglary net
# (Manual, pedagogical version)

# P(B): shape (2,)  [P(B=F), P(B=T)]
f_B = np.array([0.999, 0.001])

# P(E): shape (2,)
f_E = np.array([0.998, 0.002])

# P(A|B,E): shape (2,2,2)  dims: [A, B, E]
f_A = np.array([
    [[0.999, 0.71], [0.06, 0.05]],   # A=F
    [[0.001, 0.29], [0.94, 0.95]],   # A=T
])  # [a, b, e]

# P(J|A): shape (2,2) dims [J,A], clamp J=True → take row 1
f_J_given_A = np.array([0.05, 0.90])  # P(J=T|A=F), P(J=T|A=T)

# P(M|A): shape (2,) clamped to M=True
f_M_given_A = np.array([0.01, 0.70])  # P(M=T|A=F), P(M=T|A=T)

# Step 1: Combine J and M evidence with A
# f_AJM[a] = P(J=T|a) * P(M=T|a)
f_AJM = f_J_given_A * f_M_given_A  # shape (2,)

# Step 2: Sum out E  → factor over (A, B)
# f_ABE[a,b] = sum_e P(A=a|B=b,E=e)*P(E=e)
f_AB = np.einsum('abe,e->ab', f_A, f_E)  # shape (A, B)

# Step 3: Combine with J,M evidence on A, sum out A
# f_B[b] = sum_a f_AB[a,b] * f_AJM[a]
f_B_unnorm = np.einsum('ab,a->b', f_AB, f_AJM) * f_B
f_B_posterior = f_B_unnorm / f_B_unnorm.sum()
print(f'P(B=T | J=T, M=T) = {f_B_posterior[1]:.4f}')  # ~0.2842

For large, densely connected networks, exact inference is intractable (NP-hard in general). Monte Carlo sampling provides an approximate alternative: estimate P(Q | e) by drawing samples and counting.

Prior sampling (simplest): Sample all variables in topological order, each from P(Xi | parents). To estimate P(Q=q | E=e): draw N samples, keep only those consistent with evidence e, count fraction where Q=q. Problem: most samples may be inconsistent with evidence — rejection sampling wastes computation when P(E=e) is small.

Likelihood weighting (smarter): Fix evidence variables to their observed values (don't sample them). For each hidden variable, sample normally. But because we skipped sampling the evidence, we need to weight each sample by the probability of the evidence given the sampled parents:

w = Π_{i ∈ evidence} P(Ei=ei | Parents(Ei))

Then P(Q=q | e) ≈ Σ_{samples s: Q(s)=q} w(s) / Σ_s w(s)

This is a form of importance sampling. It's consistent (converges to the true answer) and more efficient than rejection sampling, but weights can become very small (weight degeneracy).

Markov Chain Monte Carlo (MCMC) / Gibbs sampling (most powerful): Construct a Markov chain over assignments to all non-evidence variables. At each step, pick a variable Xi and resample it from P(Xi | Markov blanket(Xi)). After a burn-in period, the chain mixes to the stationary distribution = the posterior.

Gibbs sampling only requires computing P(Xi | MB(Xi)), which involves only the CPT of Xi and the CPTs of its children — local, cheap computations. The chain is guaranteed to converge (under mild conditions) but the burn-in and mixing time must be managed carefully.

Particle filters extend sampling to dynamic/temporal settings (see HMMs below): maintain a population ("particles") of hypothesized states, propagate them forward in time, and reweight by evidence.

import random

def likelihood_weighting(net, query_var, evidence, n_samples=100000):
    """
    net: dict of {var: {'parents': [...], 'cpd': callable(parents_dict)->float}}
    evidence: dict of {var: value}
    Returns: P(query_var=True | evidence)
    """
    topo_order = topological_sort(net)  # assume implemented
    weighted_true = 0.0
    total_weight = 0.0
    for _ in range(n_samples):
        sample = {}
        w = 1.0
        for var in topo_order:
            pa = {p: sample[p] for p in net[var]['parents']}
            p_true = net[var]['cpd'](pa)  # P(var=True | parents)
            if var in evidence:
                val = evidence[var]
                w *= p_true if val else (1 - p_true)
                sample[var] = val
            else:
                sample[var] = random.random() < p_true
        total_weight += w
        if sample[query_var]:
            weighted_true += w
    return weighted_true / total_weight

# With enough samples this converges to the exact answer.
# The variance scales as 1/N, so error is O(1/sqrt(N)).

A Hidden Markov Model extends Bayesian reasoning to sequences. The setting: a system evolves through hidden states X_1, X_2, ..., X_T that we cannot observe directly. At each time step, the system emits an observable O_t that depends on the current hidden state.

The HMM is defined by three components:

1. Initial distribution π: P(X_1 = s) for each state s
2. Transition model A: P(X_{t+1} = s' | X_t = s) — a |S|×|S| matrix
3. Emission/observation model B: P(O_t = o | X_t = s) — a |S|×|O| matrix (for discrete) or a Gaussian for continuous

The Markov assumption: X_{t+1} ⊥ X_{1:t-1} | X_t. The future depends on the past only through the current state.

The observation independence assumption: O_t ⊥ {X_{1:t-1}, O_{1:t-1}} | X_t. The current observation depends only on the current state, not history.

These two assumptions make the joint: P(X_{1:T}, O_{1:T}) = P(X_1) · Π_{t=2}^T P(X_t | X_{t-1}) · Π_{t=1}^T P(O_t | X_t)

This is simply a Bayesian network unrolled through time — a dynamic Bayesian network (DBN). Each time slice has the same structure, with edges between adjacent time slices (the transition model).

Example — robot localization: Hidden states = positions on a grid. Observations = sensor readings (wall distances). The robot moves stochastically (transition model = noisy motion). Sensor readings have noise (observation model). The HMM lets us track where the robot probably is over time.

Other applications: speech recognition (states = phonemes, observations = acoustic feature vectors), part-of-speech tagging (states = POS tags, observations = words), protein sequence modeling, financial regime detection, CpG island prediction in genomics.

import numpy as np

class HMM:
    def __init__(self, pi, A, B):
        """
        pi: (n_states,) initial distribution
        A:  (n_states, n_states) transition matrix A[i,j] = P(X_t+1=j | X_t=i)
        B:  (n_states, n_obs)    emission matrix B[i,k] = P(O_t=k | X_t=i)
        """
        self.pi = np.array(pi)
        self.A  = np.array(A)
        self.B  = np.array(B)
        self.n_states = len(pi)

    def sample(self, T):
        """Sample a hidden state sequence and observations."""
        states = [np.random.choice(self.n_states, p=self.pi)]
        obs    = [np.random.choice(self.B.shape[1], p=self.B[states[0]])]
        for _ in range(T - 1):
            s = np.random.choice(self.n_states, p=self.A[states[-1]])
            o = np.random.choice(self.B.shape[1], p=self.B[s])
            states.append(s)
            obs.append(o)
        return states, obs

# Example: 2 hidden states (sunny/rainy), 2 observations (walk/home)
pi = [0.6, 0.4]
A  = [[0.7, 0.3], [0.4, 0.6]]
B  = [[0.5, 0.5], [0.1, 0.9]]  # P(walk|sunny)=0.5, P(walk|rainy)=0.1
hmm = HMM(pi, A, B)
states, obs = hmm.sample(T=10)
print('States:', states)
print('Obs:   ', obs)

Given an HMM and a sequence of observations O_{1:T}, there are three canonical inference problems:

Problem 1: Evaluation — P(O_{1:T}): What is the probability of this observation sequence under the model? Solved by the Forward algorithm.

Problem 2: Filtering/Smoothing — P(X_t | O_{1:t}): What state are we probably in now? Solved by forward (filtering) or forward-backward (smoothing).

Problem 3: Decoding — argmax_{X_{1:T}} P(X_{1:T} | O_{1:T}): What is the most likely hidden state sequence? Solved by the Viterbi algorithm.

Forward algorithm (dynamic programming): Define α_t(s) = P(O_{1:t}, X_t=s).

Base: α_1(s) = π_s · B_{s,O_1} Recursion: α_{t+1}(s') = [Σ_s α_t(s) · A_{s,s'}] · B_{s',O_{t+1}} Result: P(O_{1:T}) = Σ_s α_T(s)

Time: O(T · |S|^2). Without this DP, naive enumeration is O(|S|^T).

Viterbi algorithm (replace sum with max): Define δ_t(s) = max_{X_{1:t-1}} P(X_{1:t-1}, X_t=s, O_{1:t}).

Base: δ_1(s) = π_s · B_{s,O_1} Recursion: δ_{t+1}(s') = [max_s δ_t(s) · A_{s,s'}] · B_{s',O_{t+1}} ψ_{t+1}(s') = argmax_s δ_t(s) · A_{s,s'} (backpointer) Trace back ψ to recover the most likely path.

Viterbi vs. Forward: swap Σ for max, store backpointers. The structure is identical — this is the general max-product / min-sum message passing, also used in LDPC codes, CRFs, and neural CTC decoding.

Backward algorithm and smoothing: Define β_t(s) = P(O_{t+1:T} | X_t=s) — the probability of future observations given current state.

Base: β_T(s) = 1 Recursion: β_t(s) = Σ_{s'} A_{s,s'} · B_{s',O_{t+1}} · β_{t+1}(s')

Smoothed posterior: P(X_t=s | O_{1:T}) = α_t(s) · β_t(s) / P(O_{1:T})

This is the forward-backward algorithm (Baum-Welch E-step), the core of HMM training.

import numpy as np

def viterbi(hmm, obs):
    """
    hmm: HMM object with pi, A, B
    obs: list of observation indices
    Returns: best_path (list of state indices)
    """
    T = len(obs)
    S = hmm.n_states
    # Work in log space to avoid underflow
    log_A = np.log(hmm.A + 1e-300)
    log_B = np.log(hmm.B + 1e-300)
    log_pi = np.log(hmm.pi + 1e-300)

    delta = np.zeros((T, S))   # max log-prob
    psi   = np.zeros((T, S), dtype=int)  # backpointers

    delta[0] = log_pi + log_B[:, obs[0]]
    for t in range(1, T):
        for s in range(S):
            trans = delta[t-1] + log_A[:, s]
            psi[t, s] = np.argmax(trans)
            delta[t, s] = trans[psi[t, s]] + log_B[s, obs[t]]

    # Trace back
    path = [np.argmax(delta[T-1])]
    for t in range(T-1, 0, -1):
        path.append(psi[t, path[-1]])
    path.reverse()
    return path

# Forward algorithm (log-sum-exp for numerical stability)
def forward(hmm, obs):
    from scipy.special import logsumexp
    T, S = len(obs), hmm.n_states
    log_alpha = np.log(hmm.pi) + np.log(hmm.B[:, obs[0]])
    for t in range(1, T):
        log_alpha = np.array([
            logsumexp(log_alpha + np.log(hmm.A[:, s])) + np.log(hmm.B[s, obs[t]])
            for s in range(S)
        ])
    return logsumexp(log_alpha)  # log P(observations)

When hidden states are continuous rather than discrete, the sum over states in the forward algorithm becomes an integral — generally intractable. But if:

1. State transitions are linear Gaussian: X_{t+1} = F·X_t + w, w ~ N(0, Q)
2. Observations are linear Gaussian: O_t = H·X_t + v, v ~ N(0, R)
3. Initial state is Gaussian: X_0 ~ N(μ_0, Σ_0)

Then the posterior P(X_t | O_{1:t}) is exactly Gaussian for all t, and the forward algorithm reduces to the Kalman filter — a closed-form recursive update.

Kalman filter equations (one time step):

Predict step (propagate prior forward):

• μ_{t|t-1} = F · μ_{t-1|t-1}
• Σ_{t|t-1} = F · Σ_{t-1|t-1} · F^T + Q

Update step (incorporate observation O_t):

• Innovation (residual): y_t = O_t − H · μ_{t|t-1}
• Innovation covariance: S_t = H · Σ_{t|t-1} · H^T + R
• Kalman gain: K_t = Σ_{t|t-1} · H^T · S_t^{-1}
• μ_{t|t} = μ_{t|t-1} + K_t · y_t
• Σ_{t|t} = (I − K_t · H) · Σ_{t|t-1}

The Kalman gain K_t acts as an optimal blend: when sensor noise R is small, K is large (trust the sensor); when process noise Q is small, K is small (trust the prediction).

Applications: GPS receivers, IMU fusion in drones and phones, target tracking in radar systems, econometrics, climate modeling. The Kalman filter is one of the most widely deployed algorithms in engineering — it runs inside virtually every navigation system.

Nonlinear extensions: When dynamics or observations are nonlinear, use the Extended Kalman Filter (EKF, linearize via Jacobians) or the Unscented Kalman Filter (UKF, propagate sigma points). For highly non-Gaussian posteriors, use particle filters (sequential Monte Carlo).

import numpy as np

def kalman_filter(F, H, Q, R, mu0, Sigma0, observations):
    """
    F: state transition matrix
    H: observation matrix  
    Q: process noise covariance
    R: observation noise covariance
    Returns: list of (mu, Sigma) filtered posteriors
    """
    mu, Sigma = mu0.copy(), Sigma0.copy()
    results = []
    for z in observations:
        # --- Predict ---
        mu_pred    = F @ mu
        Sigma_pred = F @ Sigma @ F.T + Q
        # --- Update ---
        y = z - H @ mu_pred              # innovation
        S = H @ Sigma_pred @ H.T + R     # innovation covariance
        K = Sigma_pred @ H.T @ np.linalg.inv(S)  # Kalman gain
        mu    = mu_pred    + K @ y
        Sigma = (np.eye(len(mu)) - K @ H) @ Sigma_pred
        results.append((mu.copy(), Sigma.copy()))
    return results

# 1D constant-velocity tracking example
# State = [position, velocity], Observation = position with noise
dt = 1.0
F = np.array([[1, dt],[0, 1]])   # constant velocity model
H = np.array([[1, 0]])           # observe position only
Q = np.eye(2) * 0.1              # small process noise
R = np.array([[5.0]])            # sensor noise std=sqrt(5)
mu0 = np.array([0.0, 1.0])       # initial: pos=0, vel=1
S0  = np.eye(2) * 1.0

np.random.seed(42)
true_pos = [i * 1.0 for i in range(20)]
obs = [p + np.random.randn()*np.sqrt(5) for p in true_pos]
filtered = kalman_filter(F, H, Q, R, mu0, S0, obs)
for t, (mu, _) in enumerate(filtered[:5]):
    print(f't={t}: est_pos={mu[0]:.2f}, est_vel={mu[1]:.2f}, true_pos={true_pos[t]}')

So far we've assumed the CPTs are given. In practice, they must be learned from data.

Complete data (all variables observed): Given N i.i.d. training samples, the MLE for each CPT entry is simply:

P̂(Xi = xi | Pa(Xi) = pa) = Count(Xi=xi, Pa=pa) / Count(Pa=pa)

With Laplace (add-1) smoothing to avoid zero probabilities:

P̂(Xi=xi | Pa=pa) = (Count(Xi=xi, Pa=pa) + α) / (Count(Pa=pa) + α·|dom(Xi)|)

Because the log-likelihood factorizes over nodes (given the structure), each CPT can be learned independently — a huge computational simplification.

Incomplete data (latent variables): When some variables are hidden (like the hidden states in HMMs), MLE requires integrating out the hidden variables — leading to a non-convex optimization. The EM algorithm (Expectation-Maximization) handles this:

E-step: Given current parameters θ^{(old)}, compute the expected complete-data log-likelihood — equivalently, compute the posterior P(hidden | observed, θ^{(old)}) for each training example. For HMMs, this is exactly the forward-backward algorithm, computing P(X_t | O_{1:T}).

M-step: Update θ to maximize the expected log-likelihood computed in the E-step. For HMMs, this gives closed-form updates for π, A, and B (the Baum-Welch algorithm).

EM guarantees non-decreasing likelihood at each step and converges to a local maximum (not necessarily global). Initialization matters — run multiple restarts.

Structure learning: Choosing the DAG structure (which edges to include) is harder than parameter learning. Methods include:

• Score-based: search over structures, scoring each by BIC = log-likelihood − (k/2)·log(N) (penalizes complexity)
• Constraint-based: run conditional independence tests on data (PC algorithm, FCI)
• Hybrid: PC to prune, then score-based to refine

import numpy as np
from collections import Counter

def baum_welch(hmm, sequences, n_iter=20):
    """
    Train HMM parameters via EM on a list of observation sequences.
    Modifies hmm.pi, hmm.A, hmm.B in place.
    """
    S = hmm.n_states
    O = hmm.B.shape[1]
    for iteration in range(n_iter):
        # Accumulators
        gamma_sum  = np.zeros(S)
        xi_sum     = np.zeros((S, S))
        emit_sum   = np.zeros((S, O))
        gamma_init = np.zeros(S)
        for obs in sequences:
            T = len(obs)
            # Forward
            alpha = np.zeros((T, S))
            alpha[0] = hmm.pi * hmm.B[:, obs[0]]
            for t in range(1, T):
                alpha[t] = (alpha[t-1] @ hmm.A) * hmm.B[:, obs[t]]
            # Backward
            beta = np.ones((T, S))
            for t in range(T-2, -1, -1):
                beta[t] = hmm.A @ (hmm.B[:, obs[t+1]] * beta[t+1])
            # Gamma: P(X_t=s | obs)
            gamma = alpha * beta
            gamma /= gamma.sum(axis=1, keepdims=True) + 1e-300
            # Xi: P(X_t=s, X_{t+1}=s' | obs)
            for t in range(T-1):
                xi = alpha[t:t+1].T * hmm.A * hmm.B[:, obs[t+1]] * beta[t+1]
                xi_sum += xi / (xi.sum() + 1e-300)
            gamma_sum  += gamma.sum(axis=0)
            gamma_init += gamma[0]
            for t in range(T):
                emit_sum[:, obs[t]] += gamma[t]
        # M-step
        hmm.pi = gamma_init / gamma_init.sum()
        hmm.A  = xi_sum / xi_sum.sum(axis=1, keepdims=True)
        hmm.B  = emit_sum / emit_sum.sum(axis=1, keepdims=True)
    return hmm