Representation, inference, and learning in probabilistic graphical models

Why this subject matters now

Ten years ago graphical models were the dominant formalism for structured prediction, and a practitioner was expected to run junction trees and train CRFs by hand. Today most of that surface area has been absorbed into neural networks, and it is tempting to conclude the subject is historical. The opposite is true. The ideas did not disappear, they became the load-bearing abstractions underneath the systems that replaced them. The evidence lower bound that trains every VAE and underlies the modern derivation of diffusion models is a graphical-model identity. The reparameterization trick is a variational-inference technique. Probabilistic programming languages such as Stan, Pyro, and NumPyro are graphical models with the graph made implicit by a program, and their inference engines are the algorithms on this page. Graph neural networks compute a learned generalization of belief-propagation message passing. Causal inference, which every serious A/B-testing, medicine, and economics pipeline now needs, is stated entirely in the language of directed graphical models. What a practitioner is expected to know today is not how to hand triangulate a 50-node network, but why the ELBO is a bound, why the partition function makes undirected learning hard, when a sampler has actually converged, and why conditioning on a collider manufactures spurious correlation. Those are the questions this page answers with derivations rather than slogans.

Why factorization beats the curse of dimensionality

Start with the raw object. A joint distribution over \(n\) discrete variables, each with \(K\) states, is a table of \(K^n\) nonnegative numbers summing to one, hence \(K^n - 1\) free parameters. The exponent is the whole problem. For 30 binary variables the table has \(2^{30} = 1{,}073{,}741{,}824\) entries. Stored in float64 that is 8.6 GB for a toy model with thirty yes/no facts. For 100 binary variables it is \(2^{100} \approx 1.27 \times 10^{30}\) entries, more numbers than could be stored on all hardware ever manufactured. The statistical problem is worse than the storage problem. Estimating \(K^n - 1\) parameters requires on the order of \(K^n\) observations, because most cells of the table will never be observed even once. A dataset of a billion rows covers less than one entry per parameter at \(n = 30\).

Conditional independence is the structural remedy. If knowing \(Z\) renders \(X\) uninformative about \(Y\), written \(X \perp Y \mid Z\), then \(p(x, y \mid z) = p(x \mid z)\, p(y \mid z)\), and the joint over the three no longer needs a full three-way table. A graphical model is nothing more than a bookkeeping system for such statements. The graph records which conditional independences the modeler asserts, and the factorization those assertions license. The payoff compounds. If every variable depends directly on at most \(d\) others, the representation costs \(O(n K^{d+1})\) numbers instead of \(K^n\), linear in \(n\) and exponential only in the local neighborhood size. The same locality is what makes inference and learning tractable, because algorithms can work factor by factor instead of touching the exponential joint. Everything on this page, from variable elimination to the ELBO to do-calculus, is a way of exploiting that locality without ever materializing the full table.

Bayesian networks

The chain-rule factorization

The chain rule of probability holds for any distribution and any variable ordering.

So far nothing is saved. The last conditional still has \(n-1\) conditioning variables. A Bayesian network is the assertion that each variable depends directly on only a subset of its predecessors, its parents \(\mathrm{pa}(X_i)\).

with the parent sets forming a directed acyclic graph. The DAG plus one conditional probability distribution (CPD) per node is the entire model. The running example on this page is the four-node sprinkler network, used because it is small enough to verify every computation by hand and rich enough to exhibit all three kinds of dependence structure.

        Cloudy (C)
        /        \
       v          v
 Sprinkler (S)   Rain (R)        p(C, S, R, W) = p(C) p(S|C) p(R|C) p(W|S,R)
       \          /
        v        v
        WetGrass (W)

The CPDs used throughout are \(p(C{=}1) = 0.5\), then \(p(S{=}1 \mid C{=}0) = 0.5\) and \(p(S{=}1 \mid C{=}1) = 0.1\) (cloudy days discourage running the sprinkler), then \(p(R{=}1 \mid C{=}0) = 0.2\) and \(p(R{=}1 \mid C{=}1) = 0.8\), and finally \(p(W{=}1 \mid S, R)\) equal to 0, 0.9, 0.9, 0.99 for \((S,R) = (0,0), (0,1), (1,0), (1,1)\). The full joint would need \(2^4 - 1 = 15\) free parameters. The network needs \(1 + 2 + 2 + 4 = 9\). The saving is modest at four nodes and decisive at forty.

d-separation, reading independence off the graph

The factorization implies conditional independences beyond the obvious ones, and d-separation is the sound and complete graphical test for them (Pearl, 1988, with completeness in the sense that any independence not implied by d-separation fails for some parameterization of the graph). Every path between two variables is built from three primitive connections, and each behaves differently under conditioning.

chain:     A ──> B ──> C      blocked when B observed
fork:      A <── B ──> C      blocked when B observed
collider:  A ──> B <── C      blocked when B (and all its
                              descendants) UNobserved;
                              opened by observing B

Chain. \(A \to B \to C\), for example Cloudy \(\to\) Rain \(\to\) WetGrass. Marginally \(A\) and \(C\) are dependent, because information flows through \(B\). Conditioning on \(B\) blocks the path, and the factorization proves it. By Bayes' rule on the first two factors, \(p(a, c \mid b) = p(a)p(b \mid a)p(c \mid b) / p(b) = p(a \mid b)\, p(c \mid b)\). Once the intermediate cause is known, the origin adds nothing.

Fork. \(A \leftarrow B \to C\), a common cause, for example Sprinkler \(\leftarrow\) Cloudy \(\to\) Rain. The two effects are marginally dependent (both track the weather) but independent given the cause, since \(p(a, c \mid b) = p(b)p(a \mid b)p(c \mid b)/p(b) = p(a \mid b)\,p(c \mid b)\) directly from the factorization. Forks are the graphical face of confounding. Ignore \(B\) and the \(A\)-\(C\) correlation looks causal when it is not.

Collider. \(A \to B \leftarrow C\), a common effect, for example Sprinkler \(\to\) WetGrass \(\leftarrow\) Rain. This one inverts the other two. Marginally \(A\) and \(C\) are independent, again by the factorization, where \(p(a, c) = \sum_b p(a)p(c)p(b \mid a, c) = p(a)p(c)\) since the CPD sums to one. But conditioning on \(B\), or on any descendant of \(B\), couples them. This is explaining away, and it deserves a numeric demonstration because it is the least intuitive behavior in the subject. In the sprinkler network the prior probability the sprinkler ran is \(p(S{=}1) = 0.5 \cdot 0.5 + 0.5 \cdot 0.1 = 0.3\). Seeing wet grass raises it to \(p(S{=}1 \mid W{=}1) = 0.4298\) (computed in Problem 2 below). Now additionally learn that it rained, and \(p(S{=}1 \mid W{=}1, R{=}1) = 0.1945\). The rain explains the wet grass, so the evidence for the sprinkler collapses to below its prior. The two causes have become negatively dependent even though they are marginally independent. The intuition that conditioning simply adds information does not survive this example, which is why colliders cause so many real statistical accidents. Conditioning on a downstream outcome (hospital admission, dataset inclusion, being hired) silently correlates its independent causes.

A path is blocked by a conditioning set \(Z\) if it contains a chain or fork whose middle node is in \(Z\), or a collider whose middle node and all of that node's descendants are outside \(Z\). Sets \(X\) and \(Y\) are d-separated by \(Z\) when every path between them is blocked, and then \(X \perp Y \mid Z\) in every distribution that factorizes over the graph. The converse direction, that d-connection implies dependence, holds only for generic ("faithful") parameterizations. A fine-tuned CPD can cancel a dependence the graph permits.

The Markov blanket

The Markov blanket of a node \(X\) is the smallest set that shields it from the rest of the network, namely its parents, its children, and its children's other parents (the co-parents). Conditioned on the blanket, \(X\) is independent of everything else. The three components are each forced by one connection type. Parents block incoming chains and forks, children must be included because they are direct effects, and co-parents must be included because conditioning on a child opens a collider path to them. The blanket is what makes Gibbs sampling local. The conditional \(p(x_i \mid x_{-i})\) needed at every Gibbs step depends only on blanket variables, so

a product of the few CPDs that mention \(x_i\). In the sprinkler network the blanket of Sprinkler is \(\{C, W, R\}\), parent, child, and co-parent through the collider at \(W\). The Gibbs sampler in the implementation section is built from exactly these products.

I-maps, perfect maps, and equivalence classes

The relationship between graphs and distributions has an exact vocabulary. A DAG \(G\) is an I-map of a distribution \(p\) if every independence that d-separation reads off \(G\) actually holds in \(p\). The factorization theorem says \(p\) factorizes over \(G\) if and only if \(G\) is an I-map of \(p\). A complete DAG is trivially an I-map of everything (it asserts nothing), so the useful notion is a minimal I-map, from which no edge can be removed without asserting a false independence. A perfect map goes further, requiring that the graph's independences coincide exactly with the distribution's. Not every distribution has a perfect directed map. The standard counterexample is the four-cycle of dependences \(X \perp Y \mid \{W, Z\}\) and \(W \perp Z \mid \{X, Y\}\), which an undirected square captures and no DAG can.

Distinct DAGs can encode identical independence statements. \(A \to B \to C\), \(A \leftarrow B \leftarrow C\), and \(A \leftarrow B \to C\) all assert exactly \(A \perp C \mid B\), while the collider \(A \to B \leftarrow C\) asserts the different statement \(A \perp C\). The general result, due to Verma and Pearl (1990), is that two DAGs are Markov equivalent, meaning they imply the same d-separations, if and only if they share the same skeleton (edges ignoring direction) and the same v-structures (colliders whose parents are not directly connected). Equivalence classes matter for learning. Observational data can at best identify the class, not the member, which is why structure-learning algorithms return a partially directed graph (a CPDAG) and why causal direction needs interventions or extra assumptions. It also means "the arrows are causal" is an additional assumption laid on top of a Bayesian network, not something the probability model itself asserts. The causality section returns to this.

Markov random fields

The Gibbs distribution and the partition function

Some dependence structures have no natural direction. Neighboring pixels in an image constrain each other symmetrically. Atoms in a lattice interact without one causing the other. A Markov random field replaces the DAG with an undirected graph and the CPDs with nonnegative potential functions \(\psi_c\) over the cliques \(c\) of the graph.

The potentials are not probabilities. They are unnormalized compatibility scores, and the partition function \(Z\) restores normalization by summing the product over every joint configuration. That innocuous-looking sum is the central computational object of the undirected world. It couples all the potentials globally, since a change to one potential anywhere changes \(Z\). Computing it exactly is #P-hard in general, and its intractability is the reason undirected learning is hard (the gradient of the log-likelihood contains an expectation under the model, derived in the learning section) and the reason variational inference is organized around bounding \(\log Z\). Bayesian networks dodge the problem by construction. Each CPD is locally normalized, so the product of CPDs sums to one automatically and \(Z = 1\).

Conditional independence in an MRF is plain graph separation, \(X \perp Y \mid Z\) whenever every path from \(X\) to \(Y\) passes through \(Z\). The Markov blanket is just the neighbors. The two directions of the representation are tied together by the Hammersley-Clifford theorem. A strictly positive distribution satisfies the Markov properties of an undirected graph \(G\) if and only if it factorizes as a Gibbs distribution over the cliques of \(G\). The positivity condition is essential (distributions with hard zeros can satisfy the Markov properties without factorizing), and the proof, via the Möbius inversion argument in Besag's 1974 exposition, is out of scope here. Koller and Friedman give it in full. The practical content is that "local Markov structure" and "product of clique potentials" are the same modeling assumption.

The standard concrete MRF, used in the implementation section, is the Ising model, with binary spins \(x_i \in \{-1, +1\}\) on a grid,

with local fields \(h_i\) biasing individual spins and a coupling \(J > 0\) rewarding neighbor agreement. A \(3 \times 3\) grid has only \(2^9 = 512\) configurations, so \(Z\) and all marginals can be enumerated exactly, which makes it a good test bed for approximate inference, because every approximation can be scored against ground truth.

Factor graphs

Both formalisms lose information about how a distribution factorizes. An MRF clique of size three might carry one three-way potential or three pairwise ones. The undirected graph cannot say which. A factor graph makes the factorization explicit. It is a bipartite graph with a variable node per variable, a factor node per factor, and an edge whenever a factor mentions a variable. Here is the sprinkler network as a factor graph.

  (C)───[p(C)]
  (C)───[p(S|C)]───(S)
  (C)───[p(R|C)]───(R)
  (S)───[p(W|S,R)]───(R)
             │
            (W)

  variables in (parens), factor nodes in [brackets]; each factor
  connects to exactly the variables in its scope, so the same (C)
  node has edges to three factors

Factor graphs are the native data structure of message-passing algorithms. The sum-product updates below are stated most cleanly as messages flowing between variable and factor nodes, and libraries that implement belief propagation (and error-correcting decoders, where factor graphs are the standard notation) operate on this representation. Both a Bayesian network (one factor per CPD) and an MRF (one factor per clique potential) convert to a factor graph mechanically, which is why inference code only needs to be written once.

Conditional random fields

A conditional random field (Lafferty, McCallum, and Pereira, 2001) is an MRF over labels \(y\) whose potentials are allowed to look at an input \(x\) that is never modeled.

The derivation of why this form is the right one is short. A generative sequence model such as an HMM spends capacity modeling \(p(x)\), which is wasted when \(x\) is always observed, and its per-state emission factorization forbids features that span the whole input (is this word capitalized and is the previous tag a verb and does a gazetteer contain the next two words). Modeling \(p(y \mid x)\) directly removes both constraints. Any feature of the entire input may appear in any potential, because the potentials are only ever evaluated at the observed \(x\). The price is that the partition function \(Z(x)\) must be recomputed for every input. For a linear-chain CRF that is cheap. The label structure is a chain, so \(Z(x)\) and the marginals \(p(y_t, y_{t+1} \mid x)\) needed for the gradient come from exactly the forward-backward recursions derived in the temporal section, with transition scores \(\exp(\sum_k \lambda_k f_k)\) in place of transition probabilities. Training maximizes the conditional log-likelihood, whose gradient is the difference between empirical and expected feature counts, \(\partial \ell / \partial \lambda_k = \sum_t f_k(\text{data}) - \E_{p(y \mid x)}[\sum_t f_k]\), a convex problem. The same layer survives today on top of neural encoders. A BiLSTM or transformer produces the per-position scores and a CRF layer supplies the label-consistency structure and the global normalization.

From directed to undirected via moralization and triangulation

Exact inference algorithms are stated on undirected structures, so a Bayesian network must first be converted. Two steps. Moralization comes first. For every node, connect ("marry") all pairs of its parents, then drop edge directions. The added edges are forced by the factorization. The CPD \(p(w \mid s, r)\) is a single factor touching all three variables, so the undirected graph must contain a clique over \(\{W, S, R\}\), and leaving \(S\) and \(R\) unconnected would assert an independence the factor violates. Moralizing the sprinkler network adds one edge, \(S - R\). Information is lost in the direction of independences. The moral graph can no longer express the marginal independence of \(S\) and \(R\), which is exactly the collider statement undirected graphs cannot represent. Triangulation follows. Add further edges until every cycle of length four or more has a chord. Triangulated (chordal) graphs are precisely the graphs whose cliques can be arranged into a junction tree, which the exact-inference section uses. Finding the triangulation that minimizes the largest clique is NP-hard. Practical systems use greedy heuristics (min-fill, min-degree) and accept a possibly suboptimal clique size.

Exact inference

Variable elimination, worked in full

Every inference query is a sum over the joint. The query \(p(R \mid W{=}1)\) on the sprinkler network is, naively,

a sum of \(2 \times 2 = 4\) products for each value of \(R\). Variable elimination is the observation that the distributive law lets the sums be pushed inside the product, so each sum touches only the factors that mention its variable. Eliminate \(C\) first. The factors mentioning \(C\) are \(p(c)\), \(p(s \mid c)\), \(p(r \mid c)\). Multiply them and sum out \(c\) to produce a new factor

Its four entries were computed by hand in Problem 2 and are worth writing out, because seeing the intermediate factor as a concrete table is the entire point of the algorithm.

\(\tau_1\) is a distribution over \((S, R)\) here only by the accident that \(C\) was a root, which makes the entries sum to 1. In general intermediate factors are unnormalized and have no probabilistic interpretation on their own. Next eliminate \(S\). The remaining factors mentioning \(S\) are \(\tau_1(s, r)\) and the evidence-reduced \(p(W{=}1 \mid s, r)\).

No summations remain. \(\tau_2\) is the unnormalized posterior. Its normalizer \(0.189 + 0.4581 = 0.6471\) is exactly \(p(W{=}1)\), and

The implementation section reproduces this computation in PyTorch and JAX, printing the same intermediate tables, and confirms it against brute-force enumeration of the joint to an error of \(2.2 \times 10^{-16}\), one unit in the last place of float64.

Induced width, or why elimination order is everything

Each elimination step multiplies the factors touching a variable and sums it out, creating a factor over that variable's current neighbors. Graphically, eliminating a node connects all its neighbors (fill-in edges) and removes it. The cost of the step is exponential in the size of the created factor. For an elimination ordering \(\prec\), the largest factor created has \(w(\prec) + 1\) variables, where \(w(\prec)\) is the induced width of the ordering. Total time and memory are \(O(n \cdot K^{w+1})\). The minimum over orderings is the treewidth of the (moralized) graph. Everything about exact inference is contained in that quantity. Chains and trees have treewidth 1 and admit linear-time exact inference. An \(m \times m\) grid has treewidth \(m\), so exact inference on a \(100 \times 100\) Ising model would need factors with \(2^{101}\) entries and is out of the question. A fully connected graph has treewidth \(n - 1\), no better than the raw joint. Finding an optimal ordering is NP-hard, but the exponential dependence means even heuristic orderings decide feasibility. Induced width 20 on binary variables is a million-entry factor and fine, while induced width 60 is impossible on any machine.

Junction trees, all marginals at once

Variable elimination answers one query, and re-running it per query wastes shared work. The junction tree algorithm amortizes, computing marginals for every clique simultaneously for twice the cost of one elimination. The construction moralizes, triangulates, takes the maximal cliques of the triangulated graph as nodes, and connects them into a tree maximizing separator sizes (a maximum spanning tree over shared-variable counts). The defining property to verify is the running intersection property. For any variable, the set of cliques containing it forms a connected subtree. RIP is what makes local message passing globally consistent. If a variable appeared in two cliques separated by cliques not containing it, no sequence of local messages could reconcile the two copies. For the sprinkler network the triangulated moral graph has maximal cliques \(\{C, S, R\}\) and \(\{S, R, W\}\) with separator \(\{S, R\}\).

  ┌─────────┐   {S,R}   ┌─────────┐
  │ C, S, R │───────────│ S, R, W │
  └─────────┘ separator └─────────┘
  holds p(C), p(S|C),   holds p(W|S,R)
  p(R|C)

Each CPD is assigned to one clique containing its scope, giving initial clique potentials. Message passing then runs a two-pass schedule. Pick any clique as root, pass messages from the leaves inward (collect), then from the root back outward (distribute). The message from clique \(i\) to neighbor \(j\) marginalizes the clique's current potential down to the separator.

and after both passes each clique's potential times its incoming messages equals the unnormalized joint marginal over that clique's variables. The tree is then calibrated, meaning neighboring cliques agree on their separator marginals. In the sprinkler tree the single message from left to right is exactly the \(\tau_1(s, r)\) table computed above, and the right clique's calibrated potential \(\tau_1(s,r)\, p(W{=}1 \mid s,r)\) contains every downstream answer. Two passes over the tree, and \(p(C \mid W{=}1)\), \(p(S \mid W{=}1)\), \(p(R \mid W{=}1)\) are all available by summing small tables. The complexity is again exponential in the largest clique, which is the treewidth plus one when the triangulation is optimal. The junction tree does not evade the treewidth barrier, it just shares the cost across queries.

Belief propagation on trees, derived and proved exact

When the graph is already a tree, cliques are edges and the junction tree machinery collapses to the sum-product algorithm on the nodes themselves. Write the model with node and edge potentials, \(p(x) \propto \prod_i \psi_i(x_i) \prod_{(i,j)} \psi_{ij}(x_i, x_j)\). Define the message from \(i\) to a neighbor \(j\).

Claim. On a tree, after messages have flowed once in each direction along every edge (the same two-pass leaf-root-leaf schedule), the belief \(b_i(x_i) \propto \psi_i(x_i) \prod_{k \in \mathrm{nb}(i)} m_{k \to i}(x_i)\) equals the exact marginal \(p(x_i)\).

Proof. Root the tree at \(i\). Removing \(i\) splits the tree into disjoint subtrees \(T_k\), one hanging from each neighbor \(k\), and every factor of the model belongs to exactly one subtree, to node \(i\) itself, or to an edge \((k, i)\). The marginal is \(p(x_i) \propto \psi_i(x_i) \prod_k \big[ \sum_{x_{T_k}} \psi_{ki}(x_k, x_i) \prod_{f \in T_k} \psi_f \big]\), where the inner sums factor across subtrees precisely because they share no variables. So it suffices to show by induction on subtree height that \(m_{k \to i}(x_i)\) equals the bracketed sum over its subtree. In the base case, a leaf \(k\) has no neighbors besides \(i\), and the message definition gives \(m_{k \to i}(x_i) = \sum_{x_k} \psi_k(x_k) \psi_{ki}(x_k, x_i)\), which is the subtree sum. For the inductive step, with internal \(k\), substitute the inductive hypothesis for each incoming message \(m_{c \to k}\) from \(k\)'s children. The message definition then interchanges the sum over \(x_k\) with the product of complete subtree sums, which is valid because the children's subtrees are disjoint, reconstructing the full sum over \(T_k\). Since the two-pass schedule computes messages in exactly the order the induction requires (children before parents inward, parents before children outward), every message is available and every belief is exact. \(\square\)

The proof used the tree property twice, and both uses fail on loopy graphs. Subtrees are no longer disjoint, so the product of partial sums double-counts shared factors. Running the same updates anyway is loopy belief propagation, the subject of the next section. On a chain, the sum-product messages specialize to the forward-backward algorithm and the max-product variant (sum replaced by max) to Viterbi. Deriving them once here means the HMM section gets them nearly for free.

Approximate inference with variational methods

Loopy belief propagation and the Bethe free energy

Run the sum-product updates on a graph with cycles, iterating until the messages stop changing, and on many graphs the beliefs converge to good approximations of the true marginals. The algorithm was used this way (notably in turbo decoding) before anyone knew what it was computing. The answer came from statistical physics. Define, for a pairwise MRF, the Bethe free energy over candidate beliefs \(b_i\), \(b_{ij}\).

where \(d_i\) is the degree of node \(i\). The entropy term counts each edge's joint entropy once and corrects for the \(d_i\) times each node is counted by its edges. On a tree this correction is exact and \(F_{\text{Bethe}}\) equals the true Gibbs free energy, whose minimum over all valid distributions is \(-\ln Z\). On a loopy graph the counting argument breaks, and \(F_{\text{Bethe}}\) is an uncontrolled approximation, not a bound in either direction. The theorem of Yedidia, Freeman, and Weiss (2005) says the fixed points of loopy BP are exactly the stationary points of \(F_{\text{Bethe}}\) subject to the local consistency constraints \(\sum_{x_j} b_{ij}(x_i, x_j) = b_i(x_i)\) and normalization. The proof forms the Lagrangian of the constrained problem and identifies the exponentiated Lagrange multipliers with BP messages. Stationarity conditions reproduce the message updates. This explains both the successes and the failures. LBP is minimizing a sensible physics-derived objective, so its answers are often accurate on sparse, weakly coupled graphs. But nothing forces the iteration to converge (messages can oscillate), the beliefs are only locally consistent (they may correspond to no joint distribution, a point the marginal-polytope view of Wainwright and Jordan makes precise), and on strongly coupled loops the Bethe entropy miscounts badly. Generalized BP climbs a hierarchy of better entropy approximations (Kikuchi cluster expansions) at higher message cost.

Mean-field variational inference, derived from KL

Variational inference turns inference into optimization. Pick a tractable family \(\mathcal{Q}\), then find \(q^* = \argmin_{q \in \mathcal{Q}} \KL(q \,\|\, p)\), where \(p\) is the intractable target (a posterior, or a Gibbs distribution known up to \(Z\)). Write the unnormalized target \(\tilde p(x) = Z \, p(x)\). Then

so minimizing the KL is the same as maximizing a lower bound \(\mathcal{F}(q)\) on the log partition function, with the gap equal to the KL itself. Note which direction of KL this is. \(\KL(q \| p)\) penalizes \(q\) for putting mass where \(p\) has little, so the optimum is mode-seeking and typically underestimates the target's spread. The mean-field family is the fully factorized one, \(q(x) = \prod_i q_i(x_i)\). To derive the coordinate update, isolate one coordinate \(j\) and write the ELBO as a functional of \(q_j\) with the others fixed.

A KL is minimized, at zero, when its arguments coincide, so the optimal coordinate is

the coordinate-ascent update, in which the log of the optimal factor is the expected log joint under all the other factors. Iterating over coordinates increases the ELBO monotonically (each step is an exact minimization), so the procedure converges to a local optimum. For the Ising model the update has a closed form. Everything in \(\ln \tilde p\) involving \(x_i\) is \(x_i \big( h_i + J \sum_{j \in \mathrm{nb}(i)} x_j \big)\). Taking the expectation over the neighbors replaces each \(x_j\) by its mean \(m_j = \E_{q_j}[x_j]\), so \(q_i(x_i) \propto \exp\big( x_i (h_i + J \sum_j m_j) \big)\) and, for spins in \(\{-1, +1\}\),

These are the classical self-consistency equations of mean-field physics, here derived purely from KL minimization. The implementation section runs them on a \(3 \times 3\) grid against exact enumeration, and the honest result is worth previewing. At weak coupling \(J = 0.15\) the marginals are accurate to about 0.013 and the ELBO sits 0.12 nats below the true \(\ln Z\). At moderate coupling \(J = 0.4\) mean field breaks the near-symmetry, commits to the all-up mode, and its marginals are off by up to 0.26 with a 0.70-nat gap. The factorized family cannot represent the correlations that carry the probability mass, and the mode-seeking KL makes it confidently wrong rather than diffusely wrong.

Structured variational families repair this by keeping tractable dependence. Factor the joint over tractable blocks (a chain, a tree, one full time-series per block) rather than single variables. The same derivation goes through, with the coordinate update producing a distribution over a block that is itself handled by exact inference inside the block. The forward-backward algorithm frequently reappears as the inner loop. The art, going back to Saul and Jordan's structured mean-field work in the 1990s, is choosing which dependencies to keep. Each retained edge buys fidelity and costs inference.

Expectation propagation

EP (Minka, 2001) approximates in the other KL direction, locally. Write the target as a product of sites, \(p(x) \propto \prod_a f_a(x)\), and approximate each site by an exponential-family term \(\tilde f_a\), so \(q(x) \propto \prod_a \tilde f_a(x)\). The refinement loop removes one site from the approximation (the cavity \(q^{\setminus a} = q / \tilde f_a\)), splices in the true site to form the tilted distribution \(\propto f_a(x)\, q^{\setminus a}(x)\), and chooses the new \(q\) to minimize \(\KL(\text{tilted} \,\|\, q)\), which for an exponential family means matching moments. The new \(\tilde f_a\) is whatever ratio makes the update hold. Because the KL runs from target to approximation, EP is inclusive. It prefers covering the tilted distribution's mass to locking onto a mode, and it tends to overestimate rather than underestimate variance, the mirror image of mean field's failure. EP is the classical method of choice for Gaussian-process classification and is the fixed-point view of loopy BP, since BP is EP with the approximating family chosen as fully factorized discrete distributions. It converges on no general guarantee, and damping is standard practice.

The ELBO and the reparameterization trick

The same bound, written for a latent-variable model \(p_\theta(x, z) = p_\theta(x \mid z)\, p(z)\) with an amortized approximate posterior \(q_\phi(z \mid x)\) produced by a network, is the bridge from this subject to modern deep generative models. Start from the marginal likelihood and insert \(q_\phi\), giving

an identity checked by expanding the KL and cancelling, since the expectation of \(\ln p_\theta(x,z) - \ln q_\phi\) plus \(\E_q[\ln q_\phi - \ln p_\theta(z \mid x)]\) equals \(\E_q[\ln p_\theta(x, z) - \ln p_\theta(z \mid x)] = \ln p_\theta(x)\). Since the KL is nonnegative, the first term is a lower bound, tight exactly when \(q_\phi\) equals the true posterior. Rearranged into its practical form,

reconstruction quality minus a regularizer keeping the posterior near the prior. Maximizing over \(\theta\) fits the model, maximizing over \(\phi\) tightens the bound, and both can be done jointly by stochastic gradient ascent, provided the gradient with respect to \(\phi\) can pass through the sampling of \(z\). The naive estimator (the score-function or REINFORCE estimator, \(\E_q[f(z) \nabla_\phi \ln q_\phi(z)]\)) is unbiased but high-variance. The reparameterization trick (Kingma and Welling, 2014, and concurrently Rezende, Mohamed, and Wierstra, 2014) removes the randomness from the parameter path. For a Gaussian posterior \(q_\phi = \N(\mu_\phi(x), \diag \sigma^2_\phi(x))\), write \(z = \mu_\phi(x) + \sigma_\phi(x) \odot \varepsilon\) with \(\varepsilon \sim \N(0, I)\), so

an expectation of an ordinary chain-rule gradient, estimated with one or a few samples at low variance. The Gaussian KL term is available in closed form, \(\KL(\N(\mu, \sigma^2) \| \N(0,1)) = \tfrac{1}{2}(\mu^2 + \sigma^2 - \ln \sigma^2 - 1)\) per dimension, obtained by integrating the log-density ratio against the first Gaussian. That is the entire VAE, mean-field amortized variational inference with the coordinate updates replaced by a learned inference network and SGD. Diffusion models are trained on the same bound applied to a fixed Markov chain of latents, which is why the graphical-model derivation is worth owning rather than renting.

Approximate inference by sampling

Rejection and importance sampling, with variance

Monte Carlo replaces sums with samples, \(\E_p[f] \approx \frac{1}{n} \sum_i f(x_i)\), unbiased with variance \(\Var_p[f]/n\) when \(x_i \sim p\). The problem is drawing from \(p\). Rejection sampling uses an envelope, a proposal \(q\) easy to sample and a constant \(M\) with \(\tilde p(x) \le M q(x)\) everywhere. Draw \(x \sim q\) and accept with probability \(\tilde p(x) / (M q(x))\). Accepted draws are exact samples from \(p\) (the acceptance thins \(q\) precisely by the density ratio), and the acceptance rate is \(\int \tilde p / M = Z / M\). The method dies in high dimension because \(M\) grows exponentially. To cover a standard Gaussian in \(d = 100\) dimensions with a proposal Gaussian only 10% wider per axis, the density ratio at the origin forces \(M = 1.1^{100} \approx 13{,}780\), an acceptance rate of about \(7 \times 10^{-5}\), and the mismatch in real problems is far worse than 10%.

Importance sampling keeps every draw and reweights.

unbiased, with variance \(\frac{1}{n} \big( \E_q[f^2 w^2] - \E_p[f]^2 \big)\). The variance is finite only when \(q\) has heavier tails than \(p \cdot |f|\). A proposal slightly lighter-tailed than the target gives an estimator that is technically unbiased and practically useless, dominated by rare enormous weights. Minimizing the variance over \(q\) (by Cauchy-Schwarz, or calculus of variations with the normalization constraint) gives the optimal proposal \(q^*(x) \propto |f(x)|\, p(x)\), which is unattainable (it requires the answer) but sets the design principle. Put proposal mass where the integrand is large, not merely where \(p\) is large. When only \(\tilde p\) is known the weights are self-normalized, \(\hat\mu = \sum_i f(x_i) \bar w_i\) with \(\bar w_i = w_i / \sum_j w_j\), which introduces \(O(1/n)\) bias but is consistent. The standard health metric is the effective sample size \(n_{\text{eff}} = 1 / \sum_i \bar w_i^2\), which equals \(n\) for uniform weights and collapses toward 1 as one weight dominates. It is the honest count of how many samples the estimate is really built on, and it reappears below both in MCMC diagnostics and in particle filtering.

Gibbs sampling, derived

Markov chain Monte Carlo abandons independent samples. Construct a Markov chain with transition kernel \(T\) whose stationary distribution is \(p\), run it, and average over the trajectory. Gibbs sampling is the special case where the kernel cycles through variables, resampling each from its full conditional \(p(x_i \mid x_{-i})\), which in a graphical model depends only on the Markov blanket, as derived earlier. Stationarity is direct. If \(x \sim p\) and coordinate \(i\) is replaced by a draw from \(p(\cdot \mid x_{-i})\), the joint law of the result is \(p(x_{-i})\, p(x_i' \mid x_{-i}) = p\) again. Each coordinate update also satisfies detailed balance with respect to \(p\) (it is a Metropolis-Hastings step whose acceptance probability works out to exactly 1, checked in Problem 5's framework), and if every conditional is strictly positive the chain is irreducible and aperiodic, so \(p\) is its unique limiting distribution. The failure mode is correlation. Gibbs moves one coordinate at a time, so strongly coupled variables produce a chain that diffuses slowly along the coupled directions. The implementation section runs Gibbs twice, on the sprinkler posterior (198,000 kept samples estimate \(p(R{=}1 \mid W{=}1)\) as 0.7082 against the exact 0.70793) and on the \(3 \times 3\) Ising model at the same \(J = 0.4\) where mean field failed, where 49,000 Gibbs sweeps get every marginal right to within 0.003. Sampling buys asymptotic correctness at the price of compute and correlation, the exact trade mean field refuses.

Metropolis-Hastings and the detailed-balance proof

The general recipe starts from state \(x\), proposes \(x' \sim q(x' \mid x)\), and accepts with probability

otherwise stay at \(x\). Only the ratio of unnormalized densities appears, which is why MCMC works when \(Z\) is unknown, because the partition function cancels. Claim. The resulting kernel \(T\) satisfies detailed balance, \(p(x)\, T(x \to x') = p(x')\, T(x' \to x)\) for all \(x \ne x'\). Proof. For \(x \ne x'\) the kernel is \(T(x \to x') = q(x' \mid x)\, \alpha(x, x')\). Consider the case \(\tilde p(x')\, q(x \mid x') \le \tilde p(x)\, q(x' \mid x)\) (the other case is symmetric). Then \(\alpha(x, x') = \frac{\tilde p(x') q(x \mid x')}{\tilde p(x) q(x' \mid x)}\) and \(\alpha(x', x) = 1\), so

using \(\tilde p = Z p\) so that \(p(x)/\tilde p(x) = 1/Z\) and \(\tilde p(x')/Z = p(x')\). Detailed balance implies stationarity by summing both sides over \(x\), giving \(\sum_x p(x) T(x \to x') = p(x') \sum_x T(x' \to x) = p(x')\). \(\square\) With irreducibility and aperiodicity (satisfied by any proposal that can reach everywhere and the nonzero rejection probability), the ergodic theorem gives \(\frac{1}{n}\sum_t f(x_t) \to \E_p[f]\) almost surely. Detailed balance is sufficient, not necessary. Well-designed non-reversible chains can mix faster, but reversibility is what makes correctness a two-line argument.

Hamiltonian Monte Carlo and NUTS

Random-walk proposals explore a \(d\)-dimensional target in steps of size \(O(1/\sqrt d)\) if the acceptance rate is to stay reasonable, so mixing time degrades badly with dimension. HMC (introduced as hybrid Monte Carlo for lattice field theory in 1987, with Neal's 2011 exposition the standard reference) fixes this with physics. Augment the state \(x\) with a momentum \(r \sim \N(0, M)\) and define the Hamiltonian

so the joint density \(\propto e^{-H}\) has \(p\) as its \(x\)-marginal. The physical picture is a frictionless puck sliding on the surface \(U(x)\). Give it a random shove (resample \(r\)) and let it glide. It exchanges potential for kinetic energy, traveling far across the landscape along a single energy contour instead of diffusing. Hamilton's equations \(\dot x = M^{-1} r\), \(\dot r = -\nabla U(x)\) conserve \(H\) exactly. Simulated for time \(L\varepsilon\) they yield a distant proposal that would be accepted with probability 1 under exact integration. Discretization uses the leapfrog integrator.

chosen because it is volume-preserving (each substep is a shear) and time-reversible (negate \(r\), rerun, recover the start). Those two properties are exactly what make the MH correction valid with acceptance probability \(\min(1, e^{-\Delta H})\), where \(\Delta H\) is the leapfrog's energy error, typically \(O(\varepsilon^2)\) per unit time, so acceptance stays high even for long trajectories. The gradient \(\nabla U\) is supplied by automatic differentiation, which is why HMC became universally practical exactly when autodiff did.

Two tuning knobs remain, step size \(\varepsilon\) and trajectory length \(L\). Too few steps and the sampler is a slow random walk. Too many and the trajectory makes a U-turn, wastefully retracing toward its start. The No-U-Turn Sampler (Hoffman and Gelman, 2014) removes \(L\). It doubles the trajectory outward in both time directions, building a balanced binary tree of leapfrog states, and stops when any subtree makes a U-turn, detected by the criterion \((x^{+} - x^{-}) \cdot r < 0\) at the trajectory's ends. Sampling a state from the tree with a scheme that preserves detailed balance yields a sampler with no hand-tuned trajectory length, and dual averaging adapts \(\varepsilon\) toward a target acceptance rate (0.8 is the customary default) during warmup. NUTS is the default sampler in Stan, PyMC, NumPyro, and BlackJAX. For continuous, differentiable, moderate-dimensional posteriors it is the strongest general-purpose tool available.

Convergence diagnostics, R-hat and effective sample size

An MCMC estimate without a convergence diagnostic is a rumor. The Gelman-Rubin statistic runs \(m\) independent chains from dispersed starts, each of length \(n\), and compares between-chain to within-chain variance for each scalar quantity. With chain means \(\bar\theta_j\), grand mean \(\bar\theta\), and within-chain variances \(s_j^2\),

If the chains have all reached the same stationary distribution, \(B\) estimates the same variance as \(W\) and \(\hat R \to 1\). Chains stuck in different regions inflate \(B\) and push \(\hat R\) above 1. The modern refinement (Vehtari et al., 2021) applies \(\hat R\) to rank-normalized and folded quantities and splits each chain in half to detect within-chain drift, with 1.01 as the accepted threshold. The companion quantity is the effective sample size. Correlated draws carry less information than independent ones, and for a chain with lag-\(t\) autocorrelations \(\rho_t\) the variance of the mean is inflated by \(1 + 2\sum_{t \ge 1} \rho_t\), so

estimated in practice by summing empirical autocorrelations until they become noise. Reporting a posterior mean with \(n_{\text{eff}} = 40\) from a million draws is not pedantry. It is the difference between an estimate and an anecdote.

Sequential Monte Carlo and particle filtering

State-space models need posteriors over a growing latent trajectory, updated online as observations arrive. Sequential Monte Carlo represents \(p(x_{1:t} \mid y_{1:t})\) by \(N\) weighted particles and updates them recursively. Propagate each particle through a proposal \(x_t^{(i)} \sim q(x_t \mid x_{t-1}^{(i)}, y_t)\) and update its weight by the importance-sampling identity applied to the one-step-extended target.

The bootstrap filter (Gordon, Salmond, and Smith, 1993) takes \(q = p(x_t \mid x_{t-1})\), reducing the weight update to the likelihood \(p(y_t \mid x_t^{(i)})\). Left alone, the weights degenerate. Their variance grows with \(t\) until one particle carries everything, the same pathology importance sampling always develops in high (here, growing) dimension. The fix is resampling. When \(n_{\text{eff}} = 1/\sum_i (\bar w_t^{(i)})^2\) drops below a threshold (commonly \(N/2\)), draw \(N\) particles from the current weighted set with replacement and reset weights to uniform, killing improbable trajectories and cloning probable ones. The cost is path degeneracy (ancestral histories coalesce), which particle smoothing and more sophisticated kernels address. SMC is the standard tool for nonlinear, non-Gaussian state spaces, such as robot localization, tracking, econometric volatility models, and, in its particle-MCMC form, parameter inference wrapped around a filter. A useful byproduct is that the running product of average unnormalized weights is an unbiased estimate of the marginal likelihood \(p(y_{1:t})\).

Temporal models

Hidden Markov models and forward-backward, with scaling

An HMM is the simplest dynamic Bayesian network, a latent chain \(z_1 \to z_2 \to \cdots \to z_T\) with an observation hanging off each state, \(p(z_{1:T}, x_{1:T}) = \pi_{z_1} \prod_{t=2}^T A_{z_{t-1} z_t} \prod_{t=1}^T B_{z_t}(x_t)\), where \(\pi\) is the initial distribution, \(A\) the transition matrix, and \(B\) the emission model. The chain has treewidth 1, so everything is exactly computable in \(O(T K^2)\). The algorithms below are belief propagation on a chain, specialized. The running example is small enough to verify by hand, with two states (Rainy, Sunny), three observations (walk, shop, clean), and \(\pi = (0.6, 0.4)\),

The forward pass computes \(\alpha_t(k) = p(x_{1:t}, z_t = k)\). Marginalizing the previous state and using the two conditional independences of the chain (\(x_t\) depends only on \(z_t\), and \(z_t\) only on \(z_{t-1}\)) gives

Run literally, this underflows. Each step multiplies by numbers below one, so \(\alpha_T\) is of order \(10^{-T}\) and float64 fails near \(T \approx 700\) even for this toy. The standard fix (Rabiner's 1989 tutorial is the canonical treatment) normalizes each step. Compute the unnormalized update, let \(c_t\) be its sum, and store \(\hat\alpha_t = \alpha_t^{\text{unnorm}}/c_t\). Then \(\hat\alpha_t(k) = p(z_t = k \mid x_{1:t})\) (the filtered posterior, a bonus), each \(c_t = p(x_t \mid x_{1:t-1})\), and the likelihood is recovered in log space.

For the observation sequence (walk, shop, clean), \(c_1 = 0.6 \cdot 0.1 + 0.4 \cdot 0.6 = 0.30\), giving \(\hat\alpha_1 = (0.2, 0.8)\). At step 2 the predictive is \(\hat\alpha_1 A = (0.2 \cdot 0.7 + 0.8 \cdot 0.4,\ 0.2 \cdot 0.3 + 0.8 \cdot 0.6) = (0.46, 0.54)\). Multiplying by the shop column of \(B\) gives \((0.184, 0.162)\), so \(c_2 = 0.346\) and \(\hat\alpha_2 = (0.5318, 0.4682)\). At step 3 the predictive is \((0.5595, 0.4405)\), and the clean column gives \((0.2798, 0.0440)\), so \(c_3 = 0.3238\). The likelihood is \(0.30 \times 0.346 \times 0.3238 = 0.033612\), which the brute-force sum over all \(2^3 = 8\) state paths confirms exactly.

The backward pass computes \(\beta_t(k) = p(x_{t+1:T} \mid z_t = k)\) by the mirror-image recursion \(\beta_T = \mathbf{1}\), \(\beta_t(j) = \sum_k A_{jk} B_k(x_{t+1}) \beta_{t+1}(k)\), scaled by the same \(c_{t+1}\) so the products stay \(O(1)\). Combining,

with the scale factors cancelling by construction. The identity itself is just \(p(z_t, x_{1:T}) = p(x_{1:t}, z_t)\, p(x_{t+1:T} \mid z_t)\), the chain's conditional independence once more. For the three-step example the smoothed posteriors are \(\gamma_1 = (0.2317, 0.7683)\), \(\gamma_2 = (0.6241, 0.3759)\), \(\gamma_3 = (0.8640, 0.1360)\). The walk observation says Sunny, the shop and clean observations progressively say Rainy, and smoothing lets the later evidence revise the earlier belief. These numbers are verified against brute-force enumeration to \(2.2 \times 10^{-16}\) and against hmmlearn to \(2.2 \times 10^{-14}\) in the implementation section.

Viterbi and a worked trellis

Forward-backward answers the question of each state's marginal. The most probable single path is a different question with a different answer (the sequence of marginal argmaxes can even be an impossible path with probability zero). Replace the sum in the forward recursion by a max, in log space for stability. Then \(\delta_t(k) = \max_{z_{1:t-1}} \log p(z_{1:t-1}, z_t{=}k, x_{1:t})\), with the recursion \(\delta_t(k) = \log B_k(x_t) + \max_j \big[ \delta_{t-1}(j) + \log A_{jk} \big]\) and a backpointer \(\mathrm{bp}_t(k) = \argmax_j [\cdot]\). The optimal path is read off backward from \(\argmax_k \delta_T(k)\). Correctness is the same subtree induction as belief propagation with \((\max, +)\) replacing \((+, \times)\), since both are semiring instances of the same message-passing algorithm.

Baum-Welch is EM

Learning HMM parameters from observations alone is a latent-variable problem, and Baum-Welch is exactly the EM algorithm (derived in the learning section) applied to the chain. The E-step needs the posterior over latents, which forward-backward supplies in the form of the singleton marginals \(\gamma_t(k)\) and the pairwise marginals

The M-step maximizes the expected complete-data log-likelihood, which (as in the fully observed case below) separates into independent multinomial problems, giving expected-count updates.

Each iteration provably does not decrease the likelihood, but the surface is riddled with local optima and the labels are unidentifiable up to permutation, so multiple restarts are standard practice, and in modern use the E-step recursions survive inside neural models (the CRF layer's gradient is the same \(\gamma, \xi\) computation) more often than the full generative fit.

Kalman filtering, derived

Replace the discrete chain with linear-Gaussian dynamics and the sums become integrals that stay Gaussian.

The tool is the Gaussian conditioning lemma. Let \((a, b)\) be jointly Gaussian with means \((\mu_a, \mu_b)\) and covariance blocks \(\Sigma_{aa}, \Sigma_{ab}, \Sigma_{bb}\). Writing the joint exponent with the precision matrix \(\Lambda = \Sigma^{-1}\) and collecting terms in \(a\), the conditional exponent is \(-\tfrac12 a\T \Lambda_{aa} a + a\T\big(\Lambda_{aa}\mu_a - \Lambda_{ab}(b - \mu_b)\big) + \text{const}\), which is a Gaussian in \(a\) with covariance \(\Lambda_{aa}^{-1}\) and mean \(\mu_a - \Lambda_{aa}^{-1}\Lambda_{ab}(b - \mu_b)\). The block matrix inversion (Schur complement) identities \(\Lambda_{aa}^{-1} = \Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1} \Sigma_{ba}\) and \(-\Lambda_{aa}^{-1}\Lambda_{ab} = \Sigma_{ab}\Sigma_{bb}^{-1}\) convert this to the usable form

Now filter. Suppose \(z_{t-1} \mid x_{1:t-1} \sim \N(\mu_{t-1}, P_{t-1})\). Predict by applying the linear dynamics, so \(z_t \mid x_{1:t-1} \sim \N(\bar\mu_t, \bar P_t)\) with \(\bar\mu_t = F\mu_{t-1}\) and \(\bar P_t = F P_{t-1} F\T + Q\) (linear maps of Gaussians are Gaussian, and covariances add). Update next. The joint of \((z_t, x_t)\) given \(x_{1:t-1}\) is Gaussian with means \((\bar\mu_t, H\bar\mu_t)\) and blocks \(\Sigma_{zz} = \bar P_t\), \(\Sigma_{zx} = \bar P_t H\T\), \(\Sigma_{xx} = H \bar P_t H\T + R\). The lemma, conditioning \(z_t\) on the observed \(x_t\), gives

The Kalman gain \(K_t\) weighs prior against measurement by their precisions. A one-line sanity check in one dimension, with prior \(\N(0, 1)\), \(H = 1\), \(R = 1\), and observation \(x = 1\), gives \(K = 1/(1+1) = 0.5\) and posterior \(\N(0.5, 0.5)\), the equal-precision average with halved variance. The Kalman filter is forward-backward with the discrete tables replaced by these two Gaussian recursions. The smoothing analogue (Rauch-Tung-Striebel) is the backward pass.

Nonlinear dynamics break the Gaussian closure, and the extensions are a menu of local approximations. The extended Kalman filter linearizes \(f\) and \(h\) at the current mean with Jacobians and runs the same equations. It is cheap, standard in navigation stacks, and can diverge when the linearization is poor. The unscented Kalman filter propagates \(2d + 1\) deterministically chosen sigma points through the exact nonlinearity and refits a Gaussian, capturing means and covariances to higher order without Jacobians. When the posterior is genuinely multimodal (a robot that may be in either of two corridors), no Gaussian fit is honest and the particle filter from the sampling section is the tool that remains.

Learning

Maximum likelihood in fully observed Bayesian networks

With complete data \(\D = \{x^{(1)}, \dots, x^{(m)}\}\), the log-likelihood of a Bayesian network splits by the factorization.

where \(N_i(v, u)\) counts examples with \(X_i = v\) and parents \(u\). The objective is a sum of terms, each touching only one CPD's parameters for one parent configuration, so it maximizes separately. For one such multinomial, maximize \(\sum_v N(v) \log \theta_v\) subject to \(\sum_v \theta_v = 1\). The Lagrangian \(\sum_v N(v)\log\theta_v + \lambda(1 - \sum_v \theta_v)\) has stationarity \(N(v)/\theta_v = \lambda\), so \(\theta_v \propto N(v)\), and normalizing,

Maximum likelihood in a fully observed Bayesian network is counting and dividing. This decomposition is the practical reason directed models are pleasant. There is no joint optimization and no normalization coupling, and the work is embarrassingly parallel across CPDs. Its weakness is also visible. A parent configuration observed three times yields estimates with three-sample noise, and one observed zero times yields \(0/0\). The Bayesian smoothing below exists for exactly this.

Why MRF learning is hard, the partition function again

Write an MRF in exponential form, \(p(x; \theta) = \exp\big( \theta\T f(x) - \log Z(\theta) \big)\) with feature vector \(f\). Differentiate the log-likelihood of one example.

So the gradient is empirical minus expected feature counts, and it vanishes when the model's moments match the data's, a clean condition (the objective is concave in \(\theta\)) with a heavy cost. Every gradient step requires an expectation under the current model, which is an inference problem, generally intractable, nested inside the learning loop. The workarounds define the field. Pseudo-likelihood (Besag, 1975) replaces the joint by the product of full conditionals, \(\sum_i \log p(x_i \mid x_{-i}; \theta)\), each of which normalizes over one variable and is cheap. The estimator is consistent as \(m \to \infty\) though less efficient, and it weights local dependence over global structure. Contrastive divergence (Hinton, 2002) approximates the model expectation with a handful of Gibbs steps started at the data. CD-\(k\) runs \(k\) sweeps and uses the resulting samples for \(\E_p[f]\). The estimator is biased (it does not follow the gradient of any fixed objective), but the bias is small when the chain stays near the data distribution, and CD-1 trained the restricted Boltzmann machines of the 2006-era deep learning revival. Persistent CD keeps the chains alive across parameter updates, trading staleness for better mixing.

EM for latent variables, derived

When some variables are never observed, the log-likelihood \(\log p(x; \theta) = \log \sum_z p(x, z; \theta)\) has a sum inside the log and no closed-form maximum. Introduce any distribution \(q(z)\) and use Jensen's inequality on the concave logarithm.

the same ELBO as the variational section, with the gap again \(\KL(q \,\|\, p(z \mid x; \theta))\). EM is coordinate ascent on \(\mathcal{F}\). The E-step maximizes over \(q\) with \(\theta\) fixed. The gap vanishes at \(q(z) = p(z \mid x; \theta^{(t)})\), making the bound tight. The M-step maximizes over \(\theta\) with \(q\) fixed, which is maximizing the expected complete-data log-likelihood \(\E_q[\log p(x, z; \theta)]\), typically as easy as the fully observed problem (for HMMs, the count-based updates above with expected counts). Monotonicity follows in two inequalities.

the first because \(\mathcal{F}\) lower-bounds the likelihood everywhere, the second by the M-step's definition, the equality because the E-step closed the gap. So EM never decreases the likelihood, converging to a stationary point. Variational EM substitutes an approximate E-step when the posterior is intractable, at which point the bound no longer closes and the algorithm ascends the ELBO instead. A VAE is exactly this with amortization and minibatch gradients.

Structure learning

Learning the graph itself has two traditions. Score-based search defines a score and optimizes over DAGs. The score cannot be raw likelihood, which always prefers more edges (more parameters never fit worse). The Bayesian information criterion penalizes dimension,

and is the \(O(1)\)-error Laplace approximation to the marginal likelihood \(\log p(\D \mid G)\), so it is consistent. With enough data it recovers the true equivalence class. Because both terms decompose over families \((X_i, \mathrm{pa}_i)\), local search with edge additions, deletions, and reversals only ever rescores the families a move touches. Hill climbing with random restarts, or greedy equivalence search (GES), which searches over equivalence classes directly, are the standard algorithms. The space is super-exponential (roughly \(4.2 \times 10^{18}\) DAGs at 10 nodes) and exact search is NP-hard, so all practical methods are heuristic beyond ~30 nodes. Constraint-based methods instead ask the data conditional-independence questions. The PC algorithm (Spirtes and Glymour) starts from the complete undirected graph. It removes the edge \(X - Y\) whenever some conditioning set among neighbors, of size 0, then 1, then 2, and so on, renders them independent by a statistical test. It orients every unshielded triple \(X - Z - Y\) as a collider \(X \to Z \leftarrow Y\) exactly when \(Z\) was not in the set that separated \(X\) and \(Y\), and it propagates remaining orientations by acyclicity. The output is the CPDAG of an equivalence class, as it must be. The method is quasi-linear on sparse graphs but inherits every failure of finite-sample independence testing, and both traditions lean on the faithfulness assumption, that the data's independences are exactly the graph's, with no coincidental cancellations.

Bayesian parameter estimation and conjugacy

The Bayesian repair for small counts places a prior on each CPD. For a Bernoulli parameter with a Beta\((\alpha, \beta)\) prior, observing \(h\) ones and \(t\) zeros gives

the same family back again. This is conjugacy, which makes the update a pair of additions and the predictive \(p(x_{\text{new}}{=}1 \mid \D) = (\alpha + h)/(\alpha + \beta + h + t)\), Laplace smoothing with the prior as pseudo-counts. The multivariate version, Dirichlet priors over each multinomial CPD row, behaves identically, and with parameter independence assumptions the whole-network posterior factorizes just as the MLE did, keeping learning local. The marginal likelihood this machinery produces (the BDe score) is what serious score-based structure learning actually optimizes. BIC is its large-sample shadow. Problem 7 works a numeric case.

Causality

The ladder, and why prediction is not intervention

Pearl's ladder of causation has three rungs. Rung one is association, \(p(y \mid x)\), what observing \(X\) tells about \(Y\). Everything so far on this page lives here. Rung two is intervention, \(p(y \mid \mathrm{do}(x))\), what happens to \(Y\) when \(X\) is set by an action that overrides its usual causes. Rung three is counterfactuals, what \(Y\) would have been for this individual had \(X\) been different, given what actually happened. The rungs are genuinely different objects. In the sprinkler network, observing the wet grass raises the probability that it rained, while turning the sprinkler on with a wrench does nothing to the rain. Formally, conditioning on \(S = 1\) changes beliefs about \(S\)'s causes upstream (through the fork at \(C\)), while intervening on \(S\) severs the edge \(C \to S\), leaving upstream beliefs untouched. A structural causal model makes this precise. Each variable is assigned a mechanism \(X_i := f_i(\mathrm{pa}_i, U_i)\) with independent exogenous noises \(U_i\). The model entails an observational distribution, and the intervention \(\mathrm{do}(X_j = x^*)\) replaces the \(j\)-th mechanism by the constant \(x^*\), producing the truncated factorization

Everything in causal inference is the question of when quantities on rung two can be computed from data gathered on rung one, and prediction and intervention part ways exactly when confounding is present. A machine-learned \(p(y \mid x)\) is the wrong object for any decision that changes \(x\), no matter how accurate the model, because it faithfully includes the back-door correlation the action would destroy.

Back-door and front-door identification

A set \(Z\) satisfies the back-door criterion relative to \((X, Y)\) if no member of \(Z\) is a descendant of \(X\) and \(Z\) blocks every path from \(X\) to \(Y\) that enters \(X\) through an arrow into it. Then the adjustment formula, derivable from the truncated factorization by summing it over the non-\((X,Y,Z)\) variables and comparing to the observational factorization, is

Note what it is not, namely \(\sum_z p(z \mid x)\, p(y \mid x, z) = p(y \mid x)\). The intervention weights the strata by their population frequency \(p(z)\), not by their frequency among the treated, and this difference is the entire arithmetic content of confounding. Problem 6 works a case where it flips the conclusion. The descendant clause matters just as much. Adjusting for a collider or a mediator opens or destroys paths and creates bias where none existed.

The front-door criterion rescues a case where back-door adjustment is impossible because the confounder is unobserved. Suppose \(X \to M \to Y\) with a latent \(U\) confounding \(X\) and \(Y\), and \(M\) satisfying three conditions. It intercepts all directed paths from \(X\) to \(Y\), there is no unblocked back-door path from \(X\) to \(M\), and every back-door path from \(M\) to \(Y\) is blocked by \(X\). Chaining two back-door adjustments (identify the effect of \(X\) on \(M\), then of \(M\) on \(Y\) adjusting for \(X\), then compose) gives

A worked identification with concrete numbers. The observational data are \(p(X{=}1) = 0.5\), then \(p(M{=}1 \mid X{=}1) = 0.8\) and \(p(M{=}1 \mid X{=}0) = 0.2\), and \(p(Y{=}1 \mid M, X)\) equal to 0.1, 0.3, 0.2, 0.4 for \((M,X) = (0,0), (0,1), (1,0), (1,1)\). The inner sums are \(\sum_{x'} p(x') p(Y{=}1 \mid M{=}0, x') = 0.5(0.1) + 0.5(0.3) = 0.2\) and, for \(M{=}1\), \(0.5(0.2) + 0.5(0.4) = 0.3\). Then \(p(Y{=}1 \mid \mathrm{do}(X{=}1)) = 0.2(0.2) + 0.8(0.3) = 0.28\) and \(p(Y{=}1 \mid \mathrm{do}(X{=}0)) = 0.8(0.2) + 0.2(0.3) = 0.22\), a causal effect of \(+0.06\). The naive conditional comparison gives \(p(Y{=}1 \mid X{=}1) = 0.2(0.3) + 0.8(0.4) = 0.38\) versus \(p(Y{=}1 \mid X{=}0) = 0.8(0.1) + 0.2(0.2) = 0.12\), a \(+0.26\) association, more than four times the causal effect. The excess is the unobserved confounder's contribution, and the front-door formula removed it using only observational quantities plus the graph.

Both criteria are corollaries of do-calculus, three rules for rewriting expressions containing \(\mathrm{do}\), each licensed by a d-separation condition in a mutilated graph (\(G_{\overline{X}}\) deletes arrows into \(X\), and \(G_{\underline{X}}\) deletes arrows out of it). Rule 1 inserts or deletes observations, \(p(y \mid \mathrm{do}(x), z, w) = p(y \mid \mathrm{do}(x), w)\) if \(Y \perp Z \mid X, W\) in \(G_{\overline{X}}\). Rule 2 exchanges intervention for observation, \(p(y \mid \mathrm{do}(x), \mathrm{do}(z), w) = p(y \mid \mathrm{do}(x), z, w)\) if \(Y \perp Z \mid X, W\) in \(G_{\overline{X}\underline{Z}}\). Rule 3 deletes interventions outright under a condition in \(G_{\overline{X}\overline{Z(W)}}\). The calculus is complete (Shpitser and Pearl, 2006, and Huang and Valtorta, 2006). A causal effect is identifiable from observational data and the graph if and only if some sequence of the three rules reduces it to a do-free expression, and the ID algorithm decides this mechanically.

Instrumental variables

When no observable set blocks the back door and no front-door mediator exists, an instrument can still identify a linear effect. \(Z\) is an instrument for \(X \to Y\) if \(Z\) affects \(X\), and \(Z\) affects \(Y\) only through \(X\), and \(Z\) shares no confounder with \(Y\). In the linear SCM \(X = \alpha Z + \delta U + \varepsilon_X\), \(Y = \beta X + \gamma U + \varepsilon_Y\) with \(\Cov(Z, U) = 0\), taking the covariance of \(Y\) with \(Z\) gives \(\Cov(Z, Y) = \beta\, \Cov(Z, X) + \gamma\, \Cov(Z, U) = \beta\, \Cov(Z, X)\), so

estimable entirely from observables even though \(U\) is latent and ordinary regression of \(Y\) on \(X\) is biased by \(\gamma \Cov(X, U)\). This is the Wald estimator. Two-stage least squares is its multi-instrument generalization. The nonparametric analogue buys less. Without linearity, instruments yield bounds or effects for the "complier" subpopulation only, a caution that applies to every natural experiment dressed as an instrument.

Decision making

Influence diagrams and maximum expected utility

Extending the network with rectangle decision nodes and diamond utility nodes yields an influence diagram, and the normative rule is maximum expected utility. Choose the policy, a mapping from each decision node's observed parents to actions, that maximizes \(\E[U]\) under the induced joint. Evaluation is the same sum-product machinery with a max over decision variables interleaved in the elimination ordering (order matters, because what is observed at decision time constrains where the max may go).

   Oil? (chance) ────────> [Utility]
        │                      ^
        v                      │
   (Test result) ···> ┌─────────┐
     optional obs     │  Drill?  │  (decision)
                      └─────────┘

A worked instance. Oil is present with probability 0.3. Drilling costs 30 if dry and nets +100 if oil is struck (utility already net of cost), and not drilling is worth 0. The expected utility of drilling is \(0.3(100) - 0.7(30) = 30 - 21 = 9 > 0\), so the MEU policy drills, with value 9.

Value of information

The same diagram prices the option of observing before deciding. With perfect information about the oil, the decision maker drills only when oil is present, for expected utility \(0.3(100) + 0.7(0) = 30\). The value of perfect information is therefore \(30 - 9 = 21\), the most a perfectly informative survey could be worth, and an upper bound on what any noisy survey is worth. For a partial signal \(T\), the calculation conditions on each outcome, re-solves the decision, and averages, \(\mathrm{VOI}(T) = \E_T\big[ \max_a \E[U \mid a, T] \big] - \max_a \E[U \mid a] \ge 0\), nonnegative because the no-information policy is always still available (information never hurts a rational agent, an inequality that is just Jensen applied to the max). VOI computations are the principled answer to "which measurement should this system buy next," from medical testing to active sensing to which unit test to run first, and they fall out of the graphical model for free.

Worked problems

Implementation, verified

Everything below was run on this machine with PyTorch 2.7 and JAX 0.6 in float64, and checked against an independent reference, brute-force enumeration of the exact joint where the model is small enough, and hmmlearn 0.3.3 (an implementation with no shared code) for the HMM algorithms. The agreement numbers quoted after each block are from those runs, not aspirations.

Variable elimination on the sprinkler network

The PyTorch version implements a general discrete factor with named axes and the elimination loop. The JAX version makes the point that for a fixed elimination order, variable elimination is a tensor contraction, expressible as einsum with the summation order chosen by the ordering. Both print the same intermediate factor \(\tau_1\) derived in the text.

import torch
torch.set_default_dtype(torch.float64)

class Factor:
    """Discrete factor: `vars` names the axes of `table`, in order."""
    def __init__(self, vars, table):
        self.vars, self.table = list(vars), table

    def align(self, out_vars):
        # permute own axes to out_vars order, size-1 axes for missing vars
        mine = [v for v in out_vars if v in self.vars]
        t = self.table.permute([self.vars.index(v) for v in mine])
        shape = [t.shape[mine.index(v)] if v in self.vars else 1
                 for v in out_vars]
        return t.reshape(shape)

    def __mul__(self, other):
        out = self.vars + [v for v in other.vars if v not in self.vars]
        return Factor(out, self.align(out) * other.align(out))

    def sum_out(self, var):
        i = self.vars.index(var)
        return Factor([v for v in self.vars if v != var],
                      self.table.sum(dim=i))

    def observe(self, var, val):
        i = self.vars.index(var)
        return Factor([v for v in self.vars if v != var],
                      self.table.select(i, val))

def variable_elimination(factors, order, query):
    """Sum out `order` one variable at a time; normalize over `query`."""
    factors = list(factors)
    for z in order:
        bucket = [f for f in factors if z in f.vars]
        factors = [f for f in factors if z not in f.vars]
        prod = bucket[0]
        for f in bucket[1:]:
            prod = prod * f
        tau = prod.sum_out(z)              # the intermediate factor
        print(f"eliminated {z}: tau over {tau.vars} = {tau.table.tolist()}")
        factors.append(tau)
    prod = factors[0]
    for f in factors[1:]:
        prod = prod * f
    t = prod.align(query)
    return t / t.sum()

# sprinkler network: C -> S, C -> R, S -> W <- R
fC = Factor(["C"], torch.tensor([0.5, 0.5]))                  # (C,)
fS = Factor(["C", "S"], torch.tensor([[0.5, 0.5], [0.9, 0.1]]))
fR = Factor(["C", "R"], torch.tensor([[0.8, 0.2], [0.2, 0.8]]))
fW = Factor(["S", "R", "W"], torch.tensor(                    # (S,R,W)
    [[[1.0, 0.0], [0.1, 0.9]],
     [[0.1, 0.9], [0.01, 0.99]]]))

posterior = variable_elimination(
    [fC, fS, fR, fW.observe("W", 1)], order=["C", "S"], query=["R"])
print(posterior.tolist())   # [0.29207..., 0.70792...]
# eliminated C: tau over ['S','R'] = [[0.29, 0.41], [0.21, 0.09]]
# eliminated S: tau over ['R']     = [0.189, 0.4581]

import jax.numpy as jnp
from jax import config
config.update("jax_enable_x64", True)

pC    = jnp.array([0.5, 0.5])                        # (C,)
pS_C  = jnp.array([[0.5, 0.5], [0.9, 0.1]])          # (C,S)
pR_C  = jnp.array([[0.8, 0.2], [0.2, 0.8]])          # (C,R)
pW_SR = jnp.array([[[1.0, 0.0], [0.1, 0.9]],
                   [[0.1, 0.9], [0.01, 0.99]]])      # (S,R,W)

# VE with order [C, S] is literally a pair of contractions:
psi1 = jnp.einsum("c,cs,cr->sr", pC, pS_C, pR_C)   # tau_1(S,R)
psi2 = jnp.einsum("sr,sr->r", psi1, pW_SR[:, :, 1])  # tau_2(R), W=1
posterior = psi2 / psi2.sum()
print(psi1)        # [[0.29 0.41] [0.21 0.09]]
print(posterior)   # [0.29207232 0.70792768]
# a single einsum over all factors is VE with the order chosen
# by the contraction planner; opt_einsum finds low-width orders

For verification, brute-force enumeration of the \(2^4\)-entry joint gives \(p(R{=}1 \mid W{=}1) = 0.7079276773296247\). Both implementations agree with it to \(2.2 \times 10^{-16}\), one ulp in float64, and their printed \(\tau_1\) matches the hand-computed table in Problem 2 exactly.

Forward-backward with scaling, and Viterbi

The weather HMM from the text. The PyTorch code follows the scaled recursions index by index. The JAX code expresses both passes as lax.scan, which is the natural JAX shape for any chain recursion and compiles to a single fused loop.

import torch
torch.set_default_dtype(torch.float64)

pi = torch.tensor([0.6, 0.4])                       # (K,)
A  = torch.tensor([[0.7, 0.3], [0.4, 0.6]])         # (K,K)
B  = torch.tensor([[0.1, 0.4, 0.5],                 # (K,V)
                   [0.6, 0.3, 0.1]])

def forward_backward(obs, pi, A, B):
    """obs: (T,) long. Returns gamma (T,K) and log p(obs)."""
    T, K = len(obs), len(pi)
    alpha = torch.zeros(T, K)          # scaled: alpha[t] = p(z_t|x_1:t)
    c = torch.zeros(T)                 # c[t] = p(x_t | x_1:t-1)
    a = pi * B[:, obs[0]]
    c[0] = a.sum(); alpha[0] = a / c[0]
    for t in range(1, T):
        a = (alpha[t-1] @ A) * B[:, obs[t]]
        c[t] = a.sum(); alpha[t] = a / c[t]
    beta = torch.zeros(T, K)           # scaled backward messages
    beta[-1] = 1.0
    for t in range(T - 2, -1, -1):
        beta[t] = (A @ (B[:, obs[t+1]] * beta[t+1])) / c[t+1]
    gamma = alpha * beta
    return gamma / gamma.sum(1, keepdim=True), c.log().sum()

def viterbi(obs, pi, A, B):
    T, K = len(obs), len(pi)
    delta = torch.zeros(T, K); back = torch.zeros(T, K, dtype=torch.long)
    delta[0] = pi.log() + B[:, obs[0]].log()
    for t in range(1, T):
        s = delta[t-1].unsqueeze(1) + A.log()    # (K_prev, K_next)
        delta[t] = s.max(0).values + B[:, obs[t]].log()
        back[t] = s.argmax(0)
    path = torch.zeros(T, dtype=torch.long)
    path[-1] = delta[-1].argmax()
    for t in range(T - 2, -1, -1):
        path[t] = back[t+1, path[t+1]]
    return path, delta[-1].max()

obs = torch.tensor([0, 1, 2])          # walk, shop, clean
gamma, ll = forward_backward(obs, pi, A, B)
path, lp = viterbi(obs, pi, A, B)
print(gamma)            # [[0.2317,0.7683],[0.6241,0.3759],[0.8640,0.1360]]
print(ll.exp())         # 0.033612
print(path, lp.exp())   # [1, 0, 0] (Sunny,Rainy,Rainy), 0.01344

import jax, jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

pi = jnp.array([0.6, 0.4])
A  = jnp.array([[0.7, 0.3], [0.4, 0.6]])
B  = jnp.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])

def forward_backward(obs, pi, A, B):
    like = B[:, obs].T                          # (T,K) likelihoods
    def fwd(a_prev, lk):
        a = (a_prev @ A) * lk
        c = a.sum()
        return a / c, (a / c, c)
    a0 = pi * like[0]; c0 = a0.sum()
    _, (a_rest, c_rest) = jax.lax.scan(fwd, a0 / c0, like[1:])
    alpha = jnp.vstack([a0[None] / c0, a_rest])
    c = jnp.concatenate([c0[None], c_rest])
    def bwd(b_next, xs):
        lk, cn = xs
        b = (A @ (lk * b_next)) / cn
        return b, b
    _, b_rev = jax.lax.scan(bwd, jnp.ones_like(pi),
                            (like[1:][::-1], c[1:][::-1]))
    beta = jnp.vstack([b_rev[::-1], jnp.ones_like(pi)[None]])
    gamma = alpha * beta
    return gamma / gamma.sum(1, keepdims=True), jnp.log(c).sum()

def viterbi(obs, pi, A, B):
    loglike, logA = jnp.log(B[:, obs].T), jnp.log(A)
    def step(d_prev, lk):
        s = d_prev[:, None] + logA               # (K_prev, K_next)
        return s.max(0) + lk, s.argmax(0)
    d0 = jnp.log(pi) + loglike[0]
    dT, back = jax.lax.scan(step, d0, loglike[1:])
    def backtrace(z_next, ptr):
        z = ptr[z_next]
        return z, z
    zT = dT.argmax()
    _, rev = jax.lax.scan(backtrace, zT, back[::-1])
    return jnp.concatenate([rev[::-1], zT[None]]), dT.max()

obs = jnp.array([0, 1, 2])
gamma, ll = forward_backward(obs, pi, A, B)
path, lp = viterbi(obs, pi, A, B)
print(gamma, jnp.exp(ll))     # same numbers as the PyTorch tab
print(path, jnp.exp(lp))      # [1 0 0], 0.01344

Verification, two ways. On the three-step sequence, brute-force enumeration of all \(2^3\) state paths gives \(p(x_{1:3}) = 0.033612\) and the same smoothed marginals to \(2.2 \times 10^{-16}\), and confirms (Sunny, Rainy, Rainy) at probability 0.01344 as the best path. On a 200-step sequence sampled from the model, against hmmlearn's CategoricalHMM, the log-likelihood is \(-219.7932803732\) from both, agreeing to \(4.8 \times 10^{-13}\) (PyTorch) and \(4.5 \times 10^{-13}\) (JAX). Posterior marginals agree to a maximum absolute difference of \(2.3 \times 10^{-14}\), and the 200-step Viterbi paths are identical element for element in both frameworks. Differences at the \(10^{-13}\) level are accumulation-order effects, not logic.

Gibbs sampling on the sprinkler posterior

The sampler draws \((C, S, R) \sim p(\cdot \mid W{=}1)\) using the Markov-blanket conditionals derived earlier. Each update multiplies only the CPDs that mention the variable. The JAX tab uses lax.scan over a pre-split key array, which is the idiomatic way to run a long sequential chain under JIT.

import torch
torch.set_default_dtype(torch.float64)

pC  = torch.tensor([0.5, 0.5])
pS  = torch.tensor([[0.5, 0.5], [0.9, 0.1]])     # (C,S)
pR  = torch.tensor([[0.8, 0.2], [0.2, 0.8]])     # (C,R)
pW1 = torch.tensor([[0.0, 0.9], [0.9, 0.99]])    # P(W=1|S,R), (S,R)

def gibbs(n_iter=200_000, burn=2_000, seed=0):
    g = torch.Generator().manual_seed(seed)
    u = torch.rand(n_iter, 3, generator=g)       # all uniforms up front
    c, s, r = 1, 1, 1
    counts = torch.zeros(3)
    for it in range(n_iter):
        w = pC * pS[:, s] * pR[:, r]             # C | s, r  (blanket)
        c = int(u[it, 0] < w[1] / w.sum())
        w = pS[c] * pW1[:, r]                    # S | c, r, W=1
        s = int(u[it, 1] < w[1] / w.sum())
        w = pR[c] * pW1[s]                       # R | c, s, W=1
        r = int(u[it, 2] < w[1] / w.sum())
        if it >= burn:
            counts += torch.tensor([c, s, r], dtype=torch.float64)
    return counts / (n_iter - burn)

print(gibbs())   # [0.5755, 0.4298, 0.7082]
# exact by enumeration: [0.5758, 0.4298, 0.7079]

import jax, jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

pC  = jnp.array([0.5, 0.5])
pS  = jnp.array([[0.5, 0.5], [0.9, 0.1]])
pR  = jnp.array([[0.8, 0.2], [0.2, 0.8]])
pW1 = jnp.array([[0.0, 0.9], [0.9, 0.99]])   # P(W=1|S,R)

def gibbs_step(state, key):
    c, s, r = state
    kc, ks, kr = jax.random.split(key, 3)
    w = pC * pS[:, s] * pR[:, r]                       # C | s, r
    c = jax.random.bernoulli(kc, w[1] / w.sum()).astype(jnp.int32)
    w = pS[c] * pW1[:, r]                              # S | c, r, W=1
    s = jax.random.bernoulli(ks, w[1] / w.sum()).astype(jnp.int32)
    w = pR[c] * pW1[s]                                 # R | c, s, W=1
    r = jax.random.bernoulli(kr, w[1] / w.sum()).astype(jnp.int32)
    return (c, s, r), jnp.array([c, s, r])

keys = jax.random.split(jax.random.PRNGKey(0), 200_000)
init = (jnp.int32(1), jnp.int32(1), jnp.int32(1))
_, samples = jax.lax.scan(gibbs_step, init, keys)   # (200000, 3)
print(samples[2_000:].mean(0))   # [C, S, R | W=1] estimates

Verify against the exact enumerated posterior \((0.5758, 0.4298, 0.70793)\) for \((C, S, R)\). The PyTorch chain's 198,000 kept samples give \((0.5755, 0.4298, 0.7082)\), errors of \(3 \times 10^{-4}\), \(1 \times 10^{-4}\), and \(2 \times 10^{-4}\). The JAX chain gives 0.7083 for \(R\) and 0.4277 for \(S\), errors of \(4 \times 10^{-4}\) and \(2 \times 10^{-3}\). Both are consistent with the \(O(1/\sqrt{n})\) Monte Carlo error scale for correlated draws, which is the point. The sampler is exact only in distribution, and its accuracy is budgeted in samples.

Mean-field VI on a 3×3 Ising model, against exact enumeration

The coordinate updates \(m_i \leftarrow \tanh(h_i + J \sum_j m_j)\) derived in the variational section, on the grid with fields \(h = (0.2, -0.1, 0.3, -0.2, 0.1, 0.4, -0.3, 0.2, -0.1)\), compared against exact marginals and \(\ln Z\) from enumeration of all 512 states.

import itertools, torch
torch.set_default_dtype(torch.float64)

J = 0.4
h = torch.tensor([0.2, -0.1, 0.3, -0.2, 0.1, 0.4, -0.3, 0.2, -0.1])
edges = [(3*i+j, 3*i+j+1) for i in range(3) for j in range(2)] + \
        [(3*i+j, 3*i+j+3) for i in range(2) for j in range(3)]
nbrs = [[] for _ in range(9)]
for a, b in edges:
    nbrs[a].append(b); nbrs[b].append(a)

# exact reference: enumerate all 2^9 = 512 spin configurations
S = torch.tensor(list(itertools.product([-1.0, 1.0], repeat=9)))  # (512,9)
E = S @ h + J * sum(S[:, a] * S[:, b] for a, b in edges)
logZ = torch.logsumexp(E, 0)
p = torch.exp(E - logZ)
exact_p1 = ((p.unsqueeze(1) * S).sum(0) + 1) / 2      # exact P(x_i=+1)

# mean-field coordinate ascent
m = torch.zeros(9)
for sweep in range(200):
    m_old = m.clone()
    for i in range(9):
        m[i] = torch.tanh(h[i] + J * sum(m[j] for j in nbrs[i]))
    if (m - m_old).abs().max() < 1e-12:
        break                                          # 44 sweeps at J=0.4
mf_p1 = (m + 1) / 2

# ELBO under q: E_q[energy] + H(q), independence gives E[x_i x_j]=m_i m_j
q1 = mf_p1
H = -(q1 * q1.log() + (1 - q1) * (1 - q1).log()).sum()
elbo = (h * m).sum() + J * sum(m[a] * m[b] for a, b in edges) + H
print((mf_p1 - exact_p1).abs().max())   # 0.2638 at J=0.4
print(elbo, logZ)                       # 6.7945 vs 7.4905: gap 0.6960

import itertools, numpy as np
import jax, jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

J = 0.4
h = jnp.array([0.2, -0.1, 0.3, -0.2, 0.1, 0.4, -0.3, 0.2, -0.1])
edges = [(3*i+j, 3*i+j+1) for i in range(3) for j in range(2)] + \
        [(3*i+j, 3*i+j+3) for i in range(2) for j in range(3)]
adj = np.zeros((9, 9))
for a, b in edges:
    adj[a, b] = adj[b, a] = 1.0
adj = jnp.array(adj)

S = jnp.array(list(itertools.product([-1.0, 1.0], repeat=9)))  # (512,9)
E = S @ h + J * jnp.einsum("si,ij,sj->s", S, adj, S) / 2
logZ = jax.scipy.special.logsumexp(E)
p = jnp.exp(E - logZ)
exact_p1 = ((p[:, None] * S).sum(0) + 1) / 2

m = jnp.zeros(9)
for _ in range(200):                    # parallel (Jacobi-style) updates
    m_new = jnp.tanh(h + J * adj @ m)
    if jnp.abs(m_new - m).max() < 1e-12:
        m = m_new; break
    m = m_new
mf_p1 = (m + 1) / 2

q1 = mf_p1
H = -(q1 * jnp.log(q1) + (1 - q1) * jnp.log(1 - q1)).sum()
elbo = (h * m).sum() + J * (m @ adj @ m) / 2 + H
print(jnp.abs(mf_p1 - exact_p1).max(), logZ - elbo)   # 0.2638, 0.6960

Verification and the honest reading. At \(J = 0.4\) mean field converges in 44 sweeps to marginals wrong by up to 0.2638 (mean error 0.1825), with ELBO 6.7945 against the true \(\ln Z = 7.4905\), a 0.6960-nat gap, which by the identity \(\ln Z - \mathcal{F}(q) = \KL(q \| p)\) is precisely the KL divergence of the converged approximation. The failure is structural. The exact marginals hover between 0.42 and 0.73, but the factorized \(q\) cannot represent the correlated fluctuations, breaks the near-symmetry, and commits to the all-up mode with marginals between 0.54 and 0.91. At the weaker coupling \(J = 0.15\) the same code is accurate to 0.0133 with a 0.1221-nat gap. Gibbs sampling on the identical \(J = 0.4\) model (49,000 kept sweeps) gets every marginal within 0.0031. One model, three inference engines, and the trade is exactly as the theory predicts. Exact enumeration is exponential, mean field is fast and biased, sampling is asymptotically exact and pays in compute.

How it is done in practice

The gap between the derivations and deployed systems is mostly numerical discipline and scale. Production inference lives in log space end to end. Log-sum-exp with the max subtracted replaces every sum of products, the forward algorithm's scaling constants become running log-normalizers, and CRF training computes \(\log Z(x)\) once per sequence with a vectorized forward pass over the whole batch, shaped (batch, time, labels, labels) so the recursion is a single batched matrix-plus-reduce per step. A neural tagger with a CRF head is trained exactly as derived here. The loss is \(\log Z(x) - \text{score}(y^{\text{gold}})\), its gradient is expected minus empirical transition counts via forward-backward, and decoding is Viterbi. The encoder changed from hand features to a transformer, while the structured layer did not change at all. Kalman filters run in every phone, drone, and car at kilohertz rates, usually in square-root form (propagating a Cholesky factor of \(P\)) because the naive covariance update loses positive-definiteness in float32. Belief propagation decodes LDPC codes in Wi-Fi, 5G, and deep space links, implemented in fixed point in silicon, running the exact sum-product updates from this page on graphs engineered to have few short cycles so that loopy BP behaves. Bayesian workflow at scale runs on Stan or NumPyro, with four NUTS chains, warmup with dual-averaging step-size adaptation, and an automatic gate on rank-normalized \(\hat R \le 1.01\) and per-quantity effective sample sizes before any number is reported. Particle filters localize vacuum robots and warehouse fleets with a few hundred particles and stratified resampling. The uniform engineering lesson is that the algorithms are the ones derived above, but every product became a sum of logs, every loop became a batched tensor op, and every convergence claim became a monitored statistic.

The current research frontier

Probabilistic programming. The inference engines of Stan (Carpenter et al., 2017), Pyro, NumPyro, and PyMC have converged on autodiff plus NUTS as the workhorse, and the frontier is throughput and automation, with GPU-vectorized chains (NumPyro can run hundreds of chains in parallel via vmap), normalizing-flow preconditioning of hard posteriors, and variational-MCMC hybrids such as Pathfinder (Zhang et al., 2022). Jun Zhu's group built ZhuSuan (Shi et al., 2017) to couple probabilistic programming with deep generative models, one of several signs the field's center of gravity spans industry and academic labs on multiple continents.

Graph neural networks as learned message passing. Gilmer et al. (2017) unified a family of graph networks as message-passing neural networks, with a learned message function of the endpoint states per edge and a learned update of the aggregated messages per node. The analogy to belief propagation is precise at the level of dataflow, node-local state updated by aggregated neighbor messages along the graph, and it should be stated with its limits. BP's messages are distributions with fixed algebra (multiply, marginalize) derived from a generative model, its fixed points approximate a specific free energy, and its beliefs have probability semantics. An MPNN's messages are unconstrained learned vectors with no normalization, no partition function, and no calibration guarantee, and standard MPNNs are bounded in expressive power by the one-dimensional Weisfeiler-Leman test, a limitation BP does not share because its messages carry potentials, not colorings. A productive line (neural-augmented BP, and hybrids that learn residual corrections to BP messages for decoding and inference) treats BP as the inductive bias and learning as the correction.

Latent-variable deep models. The ELBO derived above now trains far more than VAEs. Diffusion models are hierarchical latent-variable models whose objective is a reweighted ELBO over a fixed noising chain, a view made exact in variational diffusion models (Kingma et al., 2021), and discrete-latent models rely on relaxations or score-function estimators where reparameterization fails. Amortized inference has become its own subfield, simulation-based inference (Cranmer, Brehmer, and Louppe, 2020), where neural networks trained on simulator output approximate posteriors for models with intractable likelihoods in physics and epidemiology.

Causal machine learning. Score-based structure search was recast as continuous optimization by NOTEARS (Zheng et al., 2018), turning acyclicity into a differentiable constraint, followed by a wave of gradient-based discovery methods and an equally useful corrective literature showing benchmark artifacts (variance-ordering leakage) inflated early claims. Schölkopf and colleagues at Max Planck pushed causal representation learning, recovering causal variables from raw observations rather than assuming them given (Schölkopf et al., 2021), while applied causal effect estimation consolidated in libraries (DoWhy, EconML) that combine graphical identification with machine-learned nuisance models.

Bayesian deep learning, honestly. The aspiration is posterior uncertainty over network weights. The accounting after a decade is mixed. Deep ensembles (Lakshminarayanan, Pritzel, and Blundell, 2017) remain the strongest practical baseline for calibration and out-of-distribution behavior despite being only loosely Bayesian. MC dropout (Gal and Ghahramani, 2016) is cheap and widely used but approximates a narrow posterior family. Laplace and SWAG-style Gaussian fits around SGD solutions give useful uncertainty at modest cost. And the cold-posterior effect (Wenzel et al., 2020), where tempering the posterior beats the untempered Bayes rule, suggests the priors and likelihoods in use are misspecified, a finding probed by the full-batch HMC study of Izmailov et al. (2021), which found genuine Bayesian averaging helps accuracy while the cold-posterior effect largely traces to data augmentation distorting the likelihood. The graphical-model vocabulary, what is the latent, what is the approximating family, which KL direction, is exactly the vocabulary in which these debates are conducted.

Open source to read

Repositories worth reading rather than merely importing, with the file to open first.

pgmpy/pgmpy, discrete Bayesian networks with exact and approximate inference, the closest thing to this page in executable form. Open pgmpy/inference/ExactInference.py, which contains both variable elimination and junction-tree belief propagation over the same factor abstractions used here.

pyro-ppl/pyro, deep probabilistic programming on PyTorch, built around trace-based stochastic variational inference. Open pyro/infer/svi.py to see the ELBO loop stripped to its essentials, then pyro/infer/trace_elbo.py for how the bound is actually estimated from execution traces.

pyro-ppl/numpyro, the JAX sibling, with the fastest NUTS implementation in the Python world. Open numpyro/infer/hmc.py, where the leapfrog integrator, dual averaging, and the tree building of NUTS are all readable in one file.

stan-dev/stan, the reference Bayesian inference engine. Open src/stan/mcmc/hmc/nuts/base_nuts.hpp for the production NUTS recursion, multinomial sampling over the trajectory and all.

blackjax-devs/blackjax, samplers as pure JAX functions, no modeling language attached, ideal for seeing algorithms without framework scaffolding. Open blackjax/mcmc/nuts.py.

google-deepmind/distrax, distributions in JAX with clean log-prob/sample contracts. Open distrax/_src/utils/hmm.py, a complete HMM with forward-backward, Viterbi, and posterior sampling in a few hundred lines, directly comparable to the implementations above.

jmschrei/pomegranate, probabilistic models rebuilt on PyTorch with GPU support, offering HMMs, mixtures, and Bayesian networks with a scikit-like API. Open pomegranate/hmm/dense_hmm.py.

py-why/dowhy, causal inference organized around the four-step model-identify-estimate-refute workflow, with do-calculus-based identification. Open dowhy/causal_model.py and follow one identify_effect call through.

py-why/EconML, heterogeneous treatment effect estimation with double machine learning. Open econml/dml/dml.py to see nuisance-model cross-fitting, the estimator IV and back-door theory turn into at production quality.

Common misconceptions

"d-connected implies dependent." d-separation guarantees independence. Its absence guarantees nothing about a particular distribution. A d-connecting path can carry exactly cancelling effects, leaving the variables independent despite the graph. Only for generic (faithful) parameterizations do d-connection and dependence coincide, and faithfulness is an assumption, one that constraint-based structure learning relies on and that near-cancellations in finite data routinely violate.

"Conditioning on more variables can only help." Conditioning on a collider or its descendant opens a path and manufactures dependence. Select on hospital admission and independent diseases become correlated. Adjust a regression for a mediator and the causal effect being estimated silently changes. Filter a dataset on an outcome and every upstream cause becomes entangled. Which covariates to adjust is a graph question (the back-door criterion), not a more-is-better question.

"If loopy BP converges, its answers are exact." Convergence means a stationary point of the Bethe free energy was found, an approximation whose entropy term miscounts on every loop. The beliefs are locally consistent but may correspond to no joint distribution, and can be badly wrong on strongly coupled cyclic graphs while looking perfectly confident.

"After burn-in, MCMC samples are independent draws from the posterior." They are neither. The chain approaches stationarity only asymptotically, and successive states remain correlated forever. A million draws with strong autocorrelation can carry the information of a few dozen independent samples. The effective sample size, not the raw draw count, is the denominator in the Monte Carlo error, which is why \(\hat R\) and \(n_{\text{eff}}\) gate every serious analysis.

"Variational inference is conservative because it is a bound." The bound is on the marginal likelihood, not on the honesty of the posterior. The \(\KL(q \| p)\) direction is mode-seeking. It produces approximations that are too narrow and, under multimodality, confidently centered on one mode, as the Ising run above shows with marginals off by 0.26. Bounded objective, overconfident answers. Both are true at once.

"The partition function is a bookkeeping constant." For a fixed model it is a constant. For learning it is the entire problem. The gradient of an MRF's log-likelihood is empirical minus model-expected features, and the second term is an inference query against \(Z\)'s distribution at every gradient step. Pseudo-likelihood, contrastive divergence, and score-based objectives exist precisely to avoid computing it.

"The most probable path is the sequence of most probable states." The per-step argmaxes of the smoothed marginals ignore transition structure and can concatenate into a path of probability zero. Viterbi and forward-backward answer different questions (joint mode versus marginals), and the worked trellis above shows them disagreeing at \(t = 2\) on a three-step chain.

"Given enough observational data, the causal graph is identifiable." Observational data identifies at best the Markov equivalence class, same skeleton, same v-structures. \(A \to B\) and \(A \leftarrow B\) are observationally indistinguishable without interventions, functional-form assumptions, or auxiliary variables, which is why honest structure-learning output is a CPDAG with some edges undirected.

Self-check

References

1. Koller, D. and Friedman, N. (2009). Probabilistic Graphical Models: Principles and Techniques. The comprehensive reference for everything in the representation, inference, and learning sections.
2. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Chapters 8-13 cover graphical models, EM, variational inference, and state-space models with full derivations. book page
3. Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective, and (2022, 2023) Probabilistic Machine Learning, Vols. 1-2. probml.github.io/pml-book
4. Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems. The source of d-separation and belief propagation.
5. Pearl, J. (2009). Causality: Models, Reasoning, and Inference, 2nd ed. Do-calculus, back-door and front-door criteria, structural causal models.
6. Pearl, J. and Mackenzie, D. (2018). The Book of Why. The ladder of causation, for intuition.
7. Spirtes, P., Glymour, C., and Scheines, R. (2000). Causation, Prediction, and Search, 2nd ed. The PC algorithm and constraint-based discovery.
8. Peters, J., Janzing, D., and Schölkopf, B. (2017). Elements of Causal Inference. open access
9. Wainwright, M. J. and Jordan, M. I. (2008). Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning 1(1-2). doi:10.1561/2200000001
10. Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999). An introduction to variational methods for graphical models. Machine Learning 37. doi:10.1023/A:1007665907178
11. Yedidia, J. S., Freeman, W. T., and Weiss, Y. (2005). Constructing free-energy approximations and generalized belief propagation algorithms. IEEE Transactions on Information Theory 51(7). doi:10.1109/TIT.2005.850085
12. Minka, T. P. (2001). Expectation propagation for approximate Bayesian inference. UAI. arXiv:1301.2294
13. Lafferty, J., McCallum, A., and Pereira, F. (2001). Conditional random fields: probabilistic models for segmenting and labeling sequence data. ICML. ACM DL
14. Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE 77(2). doi:10.1109/5.18626
15. Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS B 39(1). doi:10.1111/j.2517-6161.1977.tb01600.x
16. Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering 82(1). doi:10.1115/1.3662552
17. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo. arXiv:1206.1901
18. Hoffman, M. D. and Gelman, A. (2014). The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. JMLR 15. jmlr.org/papers/v15/hoffman14a
19. Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC. Bayesian Analysis 16(2). arXiv:1903.08008
20. Kingma, D. P. and Welling, M. (2014). Auto-encoding variational Bayes. ICLR. arXiv:1312.6114
21. Rezende, D. J., Mohamed, S., and Wierstra, D. (2014). Stochastic backpropagation and approximate inference in deep generative models. ICML. arXiv:1401.4082
22. Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017). Variational inference: a review for statisticians. JASA 112(518). arXiv:1601.00670
23. Hinton, G. E. (2002). Training products of experts by minimizing contrastive divergence. Neural Computation 14(8). doi:10.1162/089976602760128018
24. Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. (2017). Neural message passing for quantum chemistry. ICML. arXiv:1704.01212
25. Zheng, X., Aragam, B., Ravikumar, P., and Xing, E. P. (2018). DAGs with NO TEARS: continuous optimization for structure learning. NeurIPS. arXiv:1803.01422