Probabilistic Graphical Models

A probabilistic graphical model (PGM) is a joint probability distribution over many random variables together with a graph that says which of them interact directly. Each node is a variable, and a missing edge is an assumption of conditional independence. The assumptions let the distribution be written as a product of small factors, each touching a few variables, and that product is what makes the model possible to specify, to store and to compute with. Murphy (2023, §4.1) calls PGMs a convenient formalism for defining joint distributions and remarks that “independence diagrams” would have been a better name. Torralba, Isola and Freeman (2024, §29.1) stress the modularity: vision needs to reason about complicated situations, and factorization keeps the computation tractable.

Vision uses several kinds. Markov random fields are undirected models of pixels and labels on the image grid. Conditional random fields condition such a model on the image. Factor graphs write SLAM and bundle adjustment as products of measurement factors. State-space models, behind the Bayes filter, are chains in time. Gaussian mixtures and other latent-variable models are small directed graphs, and pictorial structures are trees of object parts. Each of those has its own article. This one covers what they share: directed models (Bayesian networks) and how to read independence from them, the comparison with undirected models, what inference is and why it is hard, belief propagation on trees and on graphs with loops, the two other workhorses of approximate inference (mean field and Gibbs sampling) in outline, and learning. An example computes exact marginals of a small image model by enumeration and measures how far loopy belief propagation, mean field and Gibbs sampling stray from them as the coupling between pixels grows.

Definition

A probabilistic graphical model is a family of joint distributions $p(x_1, \dots, x_N)$ described by a graph $G$ whose nodes are the variables: every distribution in the family factorizes according to $G$, and therefore satisfies the conditional independences that $G$ encodes. In a directed model, or Bayesian network, $G$ is a directed acyclic graph and each factor is the conditional distribution of a node given its parents. In an undirected model, or Markov random field, each factor is a non-negative potential on a clique of $G$ and the product is normalized by a constant $Z$. Murphy (2023, §4.2) notes that there is nothing inherently Bayesian about Bayesian networks: they are simply a way of defining probability distributions.

Intuition

A grey pixel on a wall can be dark paint in full light or light paint in shadow. Reflectance and illumination are causes; the pixel’s intensity is their effect. Drawn as a graph, the two causes each have an arrow into the intensity and no edge between them: the paint does not know about the light. Before the pixel is seen, learning that the wall is in shadow says nothing about the paint. Once the grey pixel has been seen, it does: if the wall is in shadow, shadow already accounts for the low intensity, so the paint is probably light. Each cause, once established, explains away the evidence for the other. The graph makes this visible before any number is assigned: two arrows meeting at an observed node connect their parents.

The same graph tells a computer how to do the arithmetic. A joint table over the reflectance of every pixel, the illumination of every pixel and every intensity would have astronomically many entries, but the factorization says that each intensity depends only on its own two causes, and sums over the table can be pushed inside the product. Message passing organizes those partial sums so that each is computed once (Torralba et al., 2024, §29.5).

Formal Definition

Directed models (Bayesian networks)

Order the nodes of a directed acyclic graph so that parents come before children. If each node is independent of its other predecessors given its parents $\mathrm{pa}(i)$, the chain rule of probability collapses to

$$p(x_1, \dots, x_N) = \prod_{i=1}^{N} p\big(x_i \mid x_{\mathrm{pa}(i)}\big),$$

a product of one conditional probability distribution (CPD) per node (Murphy, 2023, Eqs. 4.1–4.2). With $K$ states per variable an unstructured joint needs $O(K^N)$ numbers; with at most $P$ parents per node the factorized form needs $O(N K^{P+1})$, which can be exponentially fewer if the graph is sparse (§4.2.1). Each factor is already normalized, so no global constant is needed, and drawing a sample is easy: visit the nodes parents first and sample each from its CPD, which Murphy calls ancestral sampling (§4.2.5). Torralba et al. (2024, §29.3) note that directed models are often used to describe causal processes, as in the paint-and-light example.

Pearl (1986, §1.1) describes belief networks as directed acyclic graphs whose nodes are propositions or variables, whose arcs are direct dependencies and whose strengths are conditional probabilities, and uses the names belief networks, Bayesian networks and influence networks interchangeably.

Plates. When a model repeats the same sub-graph for $N$ independent data points sharing parameters $\theta$, the repeated nodes are drawn once inside a box, a plate, labelled $N$; the box is notation, and the model it stands for is the unrolled graph $p(\theta)\prod_n p(x_n \mid \theta)$ (Murphy, 2023, §4.2.8).

Reading independence from a directed graph: d-separation

Which conditional independences does a directed graph imply? Murphy (2023, §4.2.4.1) gives the d-separation criterion. A path between two nodes (ignoring arrow directions) is blocked by an observed set $C$ if it contains

1. a chain $s \to m \to t$ (or $s \leftarrow m \leftarrow t$) with $m \in C$;
2. a fork $s \leftarrow m \to t$ with $m \in C$; or
3. a collider, or v-structure, $s \to m \leftarrow t$ in which neither $m$ nor any descendant of $m$ is in $C$.

Sets $A$ and $B$ are d-separated by $C$ if every path between them is blocked, and then $x_A \perp x_B \mid x_C$ in every distribution that factorizes over the graph. The first two rules say that observing the middle of a chain or a common cause cuts the flow of information. The third is the reverse: two causes of a common effect are independent until the effect, or anything downstream of it, is observed. Murphy derives each rule from the factorization (Eqs. 4.18–4.24) and calls the third effect explaining away, inter-causal reasoning or Berkson’s paradox. Because conditioning can create dependence, the independences of a directed graph can be non-monotonic: adding variables to the conditioning set can destroy an independence (§4.5.1).

The Markov blanket of a node, the smallest set that makes it independent of everything else, is its parents, its children and its children’s other parents (the co-parents), and the node’s full conditional is proportional to its own CPD times its children’s CPDs (Murphy, 2023, Eqs. 4.26 and 4.32). Gibbs sampling and mean field both work through this local quantity.

Undirected models and factor graphs

An undirected model writes $p(x) = \frac{1}{Z}\prod_c \psi_c(x_c)$ over the cliques of the graph, and the Hammersley–Clifford theorem says that for strictly positive distributions this factorization and the graph-separation independences go together. Markov Random Field covers the definition, the Gibbs distribution and the theorem, the Ising and Potts models and Gaussian MRFs; Torralba et al. (2024, §29.2) introduce graphical models through the same undirected examples. A factor graph is bipartite, with variable nodes and factor nodes, and represents any of these products explicitly, including factors on more than two variables; Factor Graph covers it with the elimination algorithm and the Gaussian case used in SLAM.

Directed against undirected

Neither kind is more expressive than the other. A graph is a perfect map of a distribution when it encodes exactly the distribution’s conditional independences, and directed and undirected graphs are perfect maps for different sets of distributions (Murphy, 2023, §4.5.1). The v-structure $A \to C \leftarrow B$ asserts $A \perp B$ but not $A \perp B \mid C$; dropping the arrows gives $A - C - B$, which asserts the opposite, and no undirected graph says exactly what the v-structure says. Conversely, no directed graph says exactly what an undirected 4-cycle says. Distributions that both can represent perfectly have chordal (decomposable) graphs, trees and chains among them (§4.5.1).

To convert a directed model to an undirected one, moralize it: connect (“marry”) every pair of parents that share a child, then drop the arrows (Murphy, 2023, §4.5.2.1). The married parents give each CPD a clique to live in, and the moral graph asserts only independences that the directed graph also implies, though it may lose some. In the paint-and-light example, moralization joins reflectance and illumination, which is exactly the dependence that observing the intensity creates. The reverse conversion adds an observed dummy child to each potential (§4.5.2.2).

The two kinds also differ in normalization. Directed models are locally normalized, each CPD summing to one; undirected models are globally normalized through $Z$, a sum over all configurations. Murphy (2023, §4.5.3) traces the label bias problem of directed conditional models of sequences to local normalization, and notes the price of global normalization: a conditional random field gives no valid distribution over a sequence’s labels until the whole sequence is seen, and $Z$ couples all its parameters, which makes training slower.

A model can mix the two. Murphy (2023, Fig. 4.41) draws a grid model with hidden labels this way, an undirected Ising prior on the labels with a directed likelihood from each label to its pixel, which is the model of the example below.

Properties

Inference and what it costs

In a graphical model, inference means computing the distribution of query variables $Q$ given observed variables $V$, summing out the rest $R$:

$$p(x_Q \mid x_V) = \frac{\sum_{x_R} p(x_Q, x_V, x_R)}{\sum_{x_Q', x_R} p(x_Q', x_V, x_R)}$$

(Murphy, 2023, Eq. 4.35). For a single query node the answer is a posterior marginal. The other common query is the jointly most probable configuration, which Bayesian Inference and Markov Random Field discuss as the maximum a posteriori (MAP) estimate. Murphy (2023, §9.3.3.3) warns about a clash of names: in the graphical-models literature, for which he cites Pearl (1988, not read), the jointly most probable configuration is called the most probable explanation (MPE), and “MAP” means maximizing some variables while summing out others.

The cost of exact inference depends on the graph. Variable elimination sums out one variable at a time; each step multiplies the factors that mention the variable and creates a new factor on its neighbours. The order matters, because a bad order creates large intermediate factors. The size of the largest factor minus one is the induced width of the order, the smallest width over all orders is the graph’s treewidth $w$, and elimination then takes $O(N K^{w+1})$ time for $N$ variables of $K$ states (Murphy, 2023, §9.5.2 and §9.5.4). Finding the best order is NP-complete, and greedy heuristics such as min-fill are used in practice (§9.5.3). In general, exact inference in discrete models is NP-hard, and in fact #P-hard, intuitively because computing a normalizing constant amounts to counting satisfying assignments of a logic formula (§9.5.4). The junction tree algorithm turns the graph into a tree of cliques by running elimination symbolically and then passes messages on that tree, so that all posterior marginals are computed without repeating redundant work (§9.6). Torralba et al. (2024, §29.7) point out that the treewidth of the 2D Markov random fields used in vision can be intractably large. Factor Graph shows how elimination order and fill-in play out in a visual SLAM problem.

Belief propagation on trees

On a chain or tree, elimination can be organized so that every marginal is computed at once. The sum-product algorithm, or belief propagation (BP), passes a message from each node to each neighbour. For a pairwise model $p(x) \propto \prod_i \psi_i(x_i) \prod_{(i,j)} \psi_{ij}(x_i, x_j)$, the message from $i$ to $j$ and the belief at $i$ are

$$m_{ij}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j)\, \psi_i(x_i) \prod_{k \in N(i) \setminus j} m_{ki}(x_i), \qquad b_i(x_i) \propto \psi_i(x_i) \prod_{k \in N(i)} m_{ki}(x_i)$$

(Yedidia, Freeman and Weiss, 2000, Eqs. 2–3; Torralba et al., 2024, Eqs. 29.7–29.8). A message summarizes everything on the far side of the edge, which on a tree is a separate sub-tree, and the belief is the exact posterior marginal. Torralba et al. (2024, §29.5–29.6) show that this turns the cost of marginalizing an $N$-node chain from order $|x|^N$ into order $(N-1)|x|^2$, and that two sweeps, leaves to root and back, give every marginal for twice the cost of one. Replacing the sum by a maximum gives max-product BP, whose max-marginals recover the MAP configuration on a tree if they are unique (Murphy, 2023, §9.3.3.1).

Familiar algorithms are special cases. Murphy (2023, §9.3) presents the forward–backward algorithm for hidden Markov models and the Kalman smoother as message passing on a chain, notes that chain message passing serves inference in state-space models (§9.1), and relates the traceback of max-product on trees to the Viterbi algorithm for the most probable sequence (§9.2.6, §9.3.3.1); Factor Graph reports Kschischang, Frey and Loeliger’s (2001) account of the Viterbi algorithm and the Kalman filter as special cases of sum-product. The recursions of the Bayes filter, the Kalman filter and connectionist temporal classification, and the dynamic programme of pictorial structures (which its authors call a generalization of the Viterbi algorithm) all belong to this family. Torralba et al. (2024, §29.6.4) give a vision example: BP along one scanline of a stereo pair with four disparity states and a Potts prior turns noisy local evidence into a posterior that places the depth discontinuity at the object boundary.

Murphy (2023, §9.4) calls Pearl the inventor of belief propagation for trees. In 1982 he proposed that each variable be a processor that keeps its own belief parameters and exchanges messages with its neighbours, confining the analysis to trees in which every node has one parent (Pearl, 1982). In 1986 he showed that for singly connected networks, which allow several parents but no loops, beliefs can be updated by local propagation among parallel, autonomous processors in time proportional to the longest path, and warned that if loops are ignored, messages may circulate indefinitely and fail to settle in the correct equilibrium (Pearl, 1986, abstract and §2.4).

Loopy belief propagation

Running the same updates on a graph with cycles is loopy belief propagation (LBP): initialize all messages, update them in parallel or in some order, and repeat until they stop changing (Murphy, 2023, §9.4.1). It may not converge, and when it does, the beliefs are approximations. Murphy reports that Pearl recommended it as a heuristic despite his reservations and that Murphy, Weiss and Jordan (1999) found empirically that it works on various models. Damping, replacing each message by a weighted average of its new and old values, makes convergence more likely, with a weight of about 0.5 standard practice, and Murphy describes asynchronous and adaptive (residual) update orders as further ways to help convergence (§9.4.4). Convergence analysis uses the computation tree: $T$ iterations of LBP are exact inference on a tree of height $T + 1$ unrolled from the graph, and with weak enough edges the influence of the leaves fades (§9.4.4.1).

Yedidia, Freeman and Weiss (2000) explained what LBP approximates. Writing $\phi_{ij}(x_i, x_j) = \psi_{ij}(x_i, x_j)\psi_i(x_i)\psi_j(x_j)$ and $q_i$ for the number of neighbours of node $i$, the Bethe free energy of pair beliefs $b_{ij}$ and node beliefs $b_i$ is

$$F_\beta = \sum_{(i,j)} \sum_{x_i, x_j} b_{ij}(x_i, x_j)\big[\ln b_{ij}(x_i, x_j) - \ln \phi_{ij}(x_i, x_j)\big] - \sum_i (q_i - 1) \sum_{x_i} b_i(x_i)\big[\ln b_i(x_i) - \ln \psi_i(x_i)\big],$$

and their Claim 1 states that BP’s fixed points are exactly the zero-gradient points of $F_\beta$ under the normalization and marginalization constraints (their Eq. 4, §2). Empirically, they found stable fixed points to be local minima; they also show that $F_\beta$ is bounded below, so BP equations always have a fixed point, and that the Bethe free energy is exact for trees while the mean-field free energy is not (§2.1). Torralba et al. (2024, §29.7) summarize this as fixed points corresponding to minima of the Bethe free energy. At a fixed point, $-F_\beta$ serves as an approximation to $\log Z$; this reading of the free energy is ours, and the example checks it on a chain, where it must be exact. Clustering variables into larger regions gives the Kikuchi approximations and generalized belief propagation, which can be significantly more accurate at an adjustable extra cost (Yedidia et al., 2000, §3–4; Murphy, 2023, §9.4.6).

On grids, min-sum LBP is one of the optimizers that Markov Random Field compares for MAP labelling, and semi-global matching has been interpreted as min-sum BP along scanlines.

Mean field and sampling

Two other approximations are used when exact inference is too costly.

• Mean field. Replace the posterior by a fully factorized $q(x) = \prod_i q_i(x_i)$ and update each factor using the expected values of its neighbours. For the Ising-prior denoising model, with $\pm 1$ labels, coupling weights $W_{ij}$ and local evidence $L_i(\pm 1) = \log p(y_i \mid x_i = \pm 1)$, the update is $\mu_i \leftarrow \tanh\big(\sum_{j \in N(i)} W_{ij}\mu_j + \tfrac{1}{2}(L_i^+ - L_i^-)\big)$, with $\mu_i = \mathbb{E}_q[x_i]$, optionally damped to improve convergence (Murphy, 2023, Eqs. 10.100–10.101). Murphy describes the over-confidence of factorized approximations as a well-known flaw of variational inference (§10.1.2, of a Gaussian example) and notes that mean field often underestimates posterior uncertainty (§10.3.4.5). Evidence Lower Bound owns variational inference and the bound mean field maximizes.
• Gibbs sampling. Resample each variable in turn from its full conditional, which depends only on its Markov blanket (Murphy, 2023, Eq. 12.19). On a grid, the pixels of one checkerboard colour are independent given the other colour and can be sampled together (§12.3.3). Averages over the samples estimate marginals with a standard error of order $1/\sqrt{S}$ for $S$ independent samples (§11.6); successive Gibbs samples are correlated, so more are needed. Pearl (1986, §2.4) listed this “stochastic relaxation” as an alternative to propagation in networks with loops, guaranteed to converge but usually needing very long relaxation times. Markov Random Field covers its use for annealing and MPM estimates on images.

Learning

With every variable observed, the likelihood and the posterior of a Bayesian network’s parameters factorize over the CPDs, so each node’s parameters are fitted separately, by counting for tables (Murphy, 2023, §4.2.7.1–4.2.7.2). With hidden variables the log does not pass through the sum over them and the objective no longer decomposes; Murphy describes expectation–maximization, which alternates inference of the hidden variables with refitting, as a popular method for this case, alongside gradient methods (§4.2.7.4–4.2.7.6). The Gaussian mixture model is the simplest example, a two-node network from the component label to the observation. Undirected models are harder to fit because the gradient of $\log Z$ is an expectation under the model; Markov Random Field covers this and pseudolikelihood. The graph itself can also be learned from data (Murphy, 2023, §4.2.7).

Examples

Exact, loopy and sampled marginals on a 4 × 4 patch

The example builds a model small enough to solve exactly and compares three approximate inference methods with the exact answer as the coupling between pixels grows. A second part checks that belief propagation is exact on a chain, and a third shows explaining away.

• Model. A $4 \times 4$ patch of labels $x_s \in \{-1, +1\}$ with the true labels in rows 0000, 0110, 0111 and 0011 (1 meaning $+1$; 7 of the 16 pixels). Each pixel is observed as $y_s = [x_s = +1] + \varepsilon_s$ with Gaussian noise of standard deviation 0.6 (noise seeds 0–4). With an Ising prior of coupling $J$ on the 4-connected grid (24 edges), the posterior is $p(x \mid y) \propto \exp\big(J \sum_{s \sim t} x_s x_t + \sum_s h_s x_s\big)$, where $h_s = \tfrac{1}{2}[\log p(y_s \mid +1) - \log p(y_s \mid -1)] = (2y_s - 1)/(4 \cdot 0.6^2)$. In the cost form of Markov Random Field, $J$ is half a Potts penalty. “Log Z” below is the logarithm of this posterior’s normalizing sum, not of the evidence $p(y)$. The couplings are $J = 0.1, 0.25, 0.44, 0.7$ and $1.0$; 0.44 is about the critical coupling of the infinite Ising lattice that Markov Random Field quotes, which a $4 \times 4$ patch does not have. A base block at $J = 0$ has a closed-form answer, $p(x_s = +1 \mid y) = \sigma(2h_s)$.
• Exact. Enumeration of all $2^{16} = 65{,}536$ labellings, with log-sum-exp.
• Loopy BP. The sum-product updates above, all messages at once, started uniform, normalized and damped by 0.5, until no new message differs from the current one by $10^{-14}$, for at most 1000 iterations. The Bethe estimate of $\log Z$ is $-F_\beta$ at the final messages.
• Mean field. Murphy’s damped parallel update (Eq. 10.101, damping 0.5, as in his Ising example), started from $\mu = 0$, with the same stopping rule; $q(x_s = +1) = (1 + \mu_s)/2$.
• Gibbs. Checkerboard updates from $p(x_s = +1 \mid \text{rest}) = \sigma\big(2(J\eta_s + h_s)\big)$, with $\eta_s$ the sum of the neighbours’ labels (Murphy’s Eqs. 12.27 and 12.29, the latter with the local evidence), four chains from random starts, 100 discarded sweeps, and marginal estimates after 100, 1000 and 10,000 recorded sweeps.
• Chain. The same model on a $1 \times 16$ chain (15 edges, true labels 0000011111100000, noise seeds 100–104), where loopy BP is ordinary BP on a tree.
• Directed block. Reflectance $R$ (dark 0.2 or light 0.8) and illumination $L$ (shadow 0.25 or lit 1.0), each with prior 0.5, cause a pixel intensity $I \sim \mathcal{N}(\rho \ell, 0.05^2)$, with $\rho$ the reflectance and $\ell$ the illumination, observed at $I = 0.2$; all four configurations are enumerated.

The measures are the mean and maximum absolute error of each method’s marginals $P(x_s = +1 \mid y)$ against the exact ones, the Bethe $\log Z$ minus the exact $\log Z$, and a “confidence”, the mean over pixels of $|P(x_s = +1) - 0.5|$, for the exact and the mean-field marginals.

Seven checks test the code and data: the counts of edges, states and foreground pixels (C1); at $J = 0$, enumeration against the closed form, $\sum p = 1$ and $\log Z = \sum_s \log 2\cosh h_s$, within $10^{-12}$ (C2); on the chain, loopy BP’s marginals and Bethe $\log Z$ against enumeration within $10^{-12}$, with every run converged (C3); at $J = 0$, Gibbs within 0.015 of the exact marginals after 10,000 sweeps (C4); at $J = 0$, loopy BP, mean field and the Bethe $\log Z$ exact within $10^{-12}$ (C5); two runs giving identical output (C6); and in the directed block, $P(\text{dark} \mid \text{shadow}) = P(\text{dark})$ before the pixel is seen (C7). The hypotheses were written before any script for the example ran: (H1) loopy BP’s mean absolute error is larger at $J = 1.0$ than at $J = 0.1$ on every seed (H1a), and its mean over seeds does not decrease from one coupling to the next (H1b); (H2) the Bethe $\log Z$ error is larger in magnitude at $J = 1.0$ than at $J = 0.1$ on every seed; (H3) at $J = 0.7$ and $J = 1.0$ mean field’s confidence exceeds the exact confidence on every seed; (H4) at every coupling the seed-mean of Gibbs’s mean absolute error falls from 100 to 1000 to 10,000 sweeps; (H5) observing the shadow lowers $P(\text{dark} \mid I)$. H1 and H2 are our hypotheses: Murphy (2023, §9.4.4.1) ties LBP’s convergence to weak edges but does not say that its error grows with coupling. H3 is ours too: Murphy calls mean field over-confident in general, not at strong coupling in particular.

The script uses NumPy 2.3.5 and SciPy 1.16.3 and runs in about 4 seconds.

# Exact, loopy and sampled marginals on a 4 x 4 binary label patch (Ising prior, Gaussian pixel noise),
# a 1 x 16 chain where belief propagation must be exact, and explaining away in a 3-node Bayes net.
import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"; os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy.special import logsumexp, expit
from scipy.stats import norm

SIGMA = 0.6                               # pixel noise; intensities 0 (label -1) and 1 (label +1)
JS = [0.1, 0.25, 0.44, 0.7, 1.0]          # Ising couplings swept
SEEDS = range(5)
DAMP = 0.5                                # damping of LBP messages and of mean-field updates
TOL, MAXIT = 1e-14, 1000
CHAINS, BURN, SWEEPS = 4, 100, [100, 1000, 10000]
GRID_TRUTH = np.array([[0,0,0,0],[0,1,1,0],[0,1,1,1],[0,0,1,1]])
CHAIN_TRUTH = np.array([[0]*5 + [1]*6 + [0]*5])

def grid_edges(h, w):
    idx = np.arange(h*w).reshape(h, w)
    return np.array([(idx[r,c], idx[r,c+1]) for r in range(h) for c in range(w-1)] +
                    [(idx[r,c], idx[r+1,c]) for r in range(h-1) for c in range(w)])

def field(truth, rng):                    # h_s = (log p(y|+1) - log p(y|-1)) / 2
    y = truth.ravel() + SIGMA * rng.standard_normal(truth.size)
    return (2*y - 1) / (4*SIGMA**2)

N = 16
S = 1 - 2*((np.arange(2**N)[:, None] >> np.arange(N)[None, :]) & 1)    # all 65,536 states, spins +-1

def exact(J, h, E):                       # p(x|y) ∝ exp(J sum_{s~t} x_s x_t + sum_s h_s x_s)
    lw = J*(S[:, E[:,0]] * S[:, E[:,1]]).sum(1) + S @ h
    logZ = logsumexp(lw); p = np.exp(lw - logZ)
    return p @ (S == 1), logZ, p.sum()

def lbp(J, h, E):                         # synchronous damped sum-product; Bethe log Z (Yedidia et al. Eq. 4)
    s = np.array([-1., 1.]); psi = np.exp(J*np.outer(s, s)); phi = np.exp(np.outer(h, s))
    src = np.concatenate([E[:,0], E[:,1]]); dst = np.concatenate([E[:,1], E[:,0]])
    k = len(E); rev = np.concatenate([np.arange(k, 2*k), np.arange(k)])
    m = np.full((2*k, 2), 0.5)
    for it in range(1, MAXIT + 1):
        login = np.zeros((N, 2)); np.add.at(login, dst, np.log(m))
        cav = phi[src] * np.exp(login[src] - np.log(m[rev]))           # phi_i times messages into i except from j
        new = cav @ psi; new /= new.sum(1, keepdims=True)
        change = np.max(np.abs(new - m))
        m = DAMP*new + (1 - DAMP)*m
        if change < TOL: break
    login = np.zeros((N, 2)); np.add.at(login, dst, np.log(m))
    b = phi * np.exp(login); b /= b.sum(1, keepdims=True)
    i, j = E[:,0], E[:,1]
    cav_i = phi[i] * np.exp(login[i] - np.log(m[k:]))                    # m[k:] are the messages j -> i
    cav_j = phi[j] * np.exp(login[j] - np.log(m[:k]))
    bij = cav_i[:, :, None] * psi[None] * cav_j[:, None, :]; bij /= bij.sum((1, 2), keepdims=True)
    lphij = np.log(psi)[None] + np.log(phi[i])[:, :, None] + np.log(phi[j])[:, None, :]
    q = np.bincount(E.ravel(), minlength=N)
    F = np.sum(bij*(np.log(bij) - lphij)) - np.sum((q - 1)[:, None] * b*(np.log(b) - np.log(phi)))
    return b[:, 1], -F, it, change < TOL

def mean_field(J, h, A):                  # Murphy PML2 Eq. 10.101, parallel, damped
    mu = np.zeros(N)
    for it in range(1, MAXIT + 1):
        new = (1 - DAMP)*mu + DAMP*np.tanh(J*(A @ mu) + h)
        change = np.max(np.abs(new - mu)); mu = new
        if change < TOL: break
    return (1 + mu)/2, it, change < TOL

def gibbs(J, h, A, colour, rng):          # checkerboard Gibbs, CHAINS chains; running estimates of P(x_s = +1)
    X = np.where(rng.random((CHAINS, N)) < 0.5, 1, -1)
    counts, out, done = np.zeros(N), {}, 0
    for sweep in range(BURN + SWEEPS[-1]):
        for c in (colour, ~colour):
            p1 = expit(2*(J*(X @ A) + h))
            X = np.where(c, np.where(rng.random((CHAINS, N)) < p1, 1, -1), X)
        if sweep >= BURN:
            counts += (X == 1).sum(0); done += 1
            if done in SWEEPS: out[done] = counts / (CHAINS*done)
    return out

Eg, Ec = grid_edges(4, 4), grid_edges(1, 16)
A = np.zeros((N, N)); A[Eg[:,0], Eg[:,1]] = A[Eg[:,1], Eg[:,0]] = 1
colour = (np.add.outer(np.arange(4), np.arange(4)) % 2 == 0).ravel()
print(f"C1 data: {len(Eg)} grid edges, {len(Ec)} chain edges, {len(np.unique(S, axis=0))} distinct states, "
      f"foreground {GRID_TRUTH.sum()} of 16 (grid) and {CHAIN_TRUTH.sum()} of 16 (chain)")

# ---- base block, J = 0: every method must reach the exact marginals (Gibbs within Monte Carlo error)
c2 = c5 = c4 = 0.0
for seed in SEEDS:
    h = field(GRID_TRUTH, np.random.default_rng(seed))
    pe, logZ, tot = exact(0.0, h, Eg)
    c2 = max(c2, abs(tot - 1), np.max(np.abs(pe - expit(2*h))), abs(logZ - np.sum(np.log(2*np.cosh(h)))))
    pb, zb, _, _ = lbp(0.0, h, Eg); pm, _, _ = mean_field(0.0, h, A)
    c5 = max(c5, np.max(np.abs(pb - pe)), np.max(np.abs(pm - pe)), abs(zb - logZ))
    g = gibbs(0.0, h, A, colour, np.random.default_rng([seed, 0]))
    c4 = max(c4, np.max(np.abs(g[SWEEPS[-1]] - pe)))
print(f"C2 J=0 enumeration vs closed form: worst {c2:.1e} (tol 1e-12)")
print(f"C5 J=0 LBP, mean field, Bethe log Z vs exact: worst {c5:.1e} (tol 1e-12)")
print(f"C4 J=0 Gibbs, {SWEEPS[-1]} sweeps: worst max|err| {c4:.5f} (tol 0.015)")

# ---- chain: belief propagation on a tree must equal enumeration
c3m = c3z = 0.0; c3conv = True
for J in JS:
    for seed in SEEDS:
        h = field(CHAIN_TRUTH, np.random.default_rng(100 + seed))
        pe, logZ, _ = exact(J, h, Ec); pb, zb, _, conv = lbp(J, h, Ec)
        c3m = max(c3m, np.max(np.abs(pb - pe))); c3z = max(c3z, abs(zb - logZ)); c3conv &= conv
print(f"C3 chain, 25 runs: max|LBP - exact| {c3m:.1e}, max|Bethe log Z - log Z| {c3z:.1e} (tol 1e-12), "
      f"all converged {c3conv}")

# ---- 4 x 4 grid: sweep of the coupling J
R = {}
for jj, J in enumerate(JS, start=1):
    for seed in SEEDS:
        h = field(GRID_TRUTH, np.random.default_rng(seed))
        pe, logZ, _ = exact(J, h, Eg)
        pb, zb, itb, convb = lbp(J, h, Eg); pm, itm, convm = mean_field(J, h, A)
        g = gibbs(J, h, A, colour, np.random.default_rng([seed, jj]))
        r = {"lbp": pb, "mf": pm, **{f"g{n}": g[n] for n in SWEEPS}}
        R[J, seed] = dict(maxerr={k: np.max(np.abs(v - pe)) for k, v in r.items()},
                          meanerr={k: np.mean(np.abs(v - pe)) for k, v in r.items()},
                          conf_exact=np.mean(np.abs(pe - 0.5)), conf_mf=np.mean(np.abs(pm - 0.5)),
                          dz=zb - logZ, itb=itb, convb=convb, itm=itm, convm=convm)

def row(J, f, fmt="9.6f"):
    v = [f(R[J, s]) for s in SEEDS]
    return " ".join(f"{x:{fmt}}" for x in v) + f"   mean {np.mean(v):{fmt}}"
print("\nper seed 0-4, then the seed mean")
for J in JS:
    print(f"J = {J}:  LBP converged {sum(R[J, s]['convb'] for s in SEEDS)}/5 "
          f"(iterations {', '.join(str(R[J, s]['itb']) for s in SEEDS)}),  "
          f"mean field converged {sum(R[J, s]['convm'] for s in SEEDS)}/5")
    print("  LBP   mean|err| " + row(J, lambda d: d["meanerr"]["lbp"]))
    print("  LBP   max|err|  " + row(J, lambda d: d["maxerr"]["lbp"]))
    print("  Bethe logZ err  " + row(J, lambda d: d["dz"], "9.2e"))
    print("  MF    mean|err| " + row(J, lambda d: d["meanerr"]["mf"]))
    print("  MF    max|err|  " + row(J, lambda d: d["maxerr"]["mf"]))
    print("  conf. exact     " + row(J, lambda d: d["conf_exact"]))
    print("  conf. MF        " + row(J, lambda d: d["conf_mf"]))
    for n in SWEEPS:
        print(f"  Gibbs {n:>5} mean " + row(J, lambda d: d["meanerr"][f"g{n}"]))
    print(f"  Gibbs 10000 max  " + row(J, lambda d: d["maxerr"][f"g{SWEEPS[-1]}"]))

# ---- hypotheses
lo, hi = JS[0], JS[-1]
h1a = [R[hi, s]["meanerr"]["lbp"] > R[lo, s]["meanerr"]["lbp"] for s in SEEDS]
m1 = [np.mean([R[J, s]["meanerr"]["lbp"] for s in SEEDS]) for J in JS]
h2 = [abs(R[hi, s]["dz"]) > abs(R[lo, s]["dz"]) for s in SEEDS]
h3 = [R[J, s]["conf_mf"] > R[J, s]["conf_exact"] for J in (0.7, 1.0) for s in SEEDS]
h4 = [np.mean([R[J, s]["meanerr"]["g100"] for s in SEEDS]) > np.mean([R[J, s]["meanerr"]["g1000"] for s in SEEDS])
      > np.mean([R[J, s]["meanerr"]["g10000"] for s in SEEDS]) for J in JS]
print(f"\nH1a LBP mean|err| larger at J={hi} than at J={lo}: {sum(h1a)}/5 seeds")
print(f"H1b seed-mean LBP mean|err| non-decreasing in J: {all(np.diff(m1) >= 0)}")
print(f"H2  |Bethe log Z error| larger at J={hi} than at J={lo}: {sum(h2)}/5 seeds")
print(f"H3  mean-field confidence above exact at J=0.7 and 1.0: {sum(h3)}/10 (J, seed) pairs")
print(f"H4  Gibbs seed-mean error falls 100 > 1000 > 10000 sweeps: {sum(h4)}/5 couplings")

# ---- directed block: reflectance R and illumination L -> pixel intensity I (explaining away)
rho = {"dark": 0.2, "light": 0.8}; lum = {"shadow": 0.25, "lit": 1.0}; I_OBS, NOISE = 0.2, 0.05
joint = {(r, l): 0.5 * 0.5 * norm.pdf(I_OBS, rho[r]*lum[l], NOISE) for r in rho for l in lum}
prior = {(r, l): 0.25 for r in rho for l in lum}
def p_dark(tab, l=None):
    keep = {k: v for k, v in tab.items() if l is None or k[1] == l}
    return sum(v for k, v in keep.items() if k[0] == "dark") / sum(keep.values())
print(f"\nC7 before the pixel is seen: P(dark) {p_dark(prior):.4f}, P(dark | shadow) {p_dark(prior, 'shadow'):.4f}")
print(f"after I = {I_OBS}: P(dark | I) {p_dark(joint):.4f}, P(dark | I, shadow) {p_dark(joint, 'shadow'):.4f}, "
      f"P(dark | I, lit) {p_dark(joint, 'lit'):.4f}")
print(f"H5 observing the shadow lowers P(dark | I): {p_dark(joint, 'shadow') < p_dark(joint)}")

Output:

C1 data: 24 grid edges, 15 chain edges, 65536 distinct states, foreground 7 of 16 (grid) and 6 of 16 (chain)
C2 J=0 enumeration vs closed form: worst 1.1e-14 (tol 1e-12)
C5 J=0 LBP, mean field, Bethe log Z vs exact: worst 1.3e-14 (tol 1e-12)
C4 J=0 Gibbs, 10000 sweeps: worst max|err| 0.00571 (tol 0.015)
C3 chain, 25 runs: max|LBP - exact| 1.2e-14, max|Bethe log Z - log Z| 1.4e-14 (tol 1e-12), all converged True

per seed 0-4, then the seed mean
J = 0.1:  LBP converged 5/5 (iterations 49, 50, 48, 49, 49),  mean field converged 5/5
  LBP   mean|err|  0.000011  0.000014  0.000005  0.000006  0.000012   mean  0.000009
  LBP   max|err|   0.000035  0.000060  0.000013  0.000030  0.000040   mean  0.000036
  Bethe logZ err  -1.01e-04 -1.11e-04 -3.06e-05 -4.74e-05 -1.08e-04   mean -7.96e-05
  MF    mean|err|  0.002447  0.002425  0.001862  0.001739  0.002234   mean  0.002141
  MF    max|err|   0.005742  0.004834  0.003547  0.003940  0.004767   mean  0.004566
  conf. exact      0.285516  0.270111  0.328163  0.340731  0.314168   mean  0.307738
  conf. MF         0.287963  0.272535  0.330019  0.342470  0.316402   mean  0.309878
  Gibbs   100 mean  0.012308  0.013379  0.010910  0.011790  0.010218   mean  0.011721
  Gibbs  1000 mean  0.005075  0.004940  0.005738  0.003525  0.005679   mean  0.004991
  Gibbs 10000 mean  0.001505  0.001839  0.001257  0.001269  0.001029   mean  0.001380
  Gibbs 10000 max   0.004067  0.005158  0.003836  0.003666  0.004748   mean  0.004295
J = 0.25:  LBP converged 5/5 (iterations 66, 67, 59, 61, 66),  mean field converged 5/5
  LBP   mean|err|  0.000396  0.000462  0.000166  0.000192  0.000467   mean  0.000337
  LBP   max|err|   0.001254  0.002006  0.000377  0.001082  0.001809   mean  0.001306
  Bethe logZ err  -3.19e-03 -3.27e-03 -9.88e-04 -1.67e-03 -4.19e-03   mean -2.66e-03
  MF    mean|err|  0.015344  0.013410  0.010846  0.009962  0.013234   mean  0.012559
  MF    max|err|   0.036840  0.031766  0.020441  0.023479  0.045318   mean  0.031569
  conf. exact      0.299192  0.279521  0.341301  0.339544  0.303461   mean  0.312604
  conf. MF         0.314536  0.292505  0.352146  0.348847  0.316177   mean  0.324842
  Gibbs   100 mean  0.016565  0.016266  0.011923  0.010544  0.015824   mean  0.014224
  Gibbs  1000 mean  0.004387  0.006534  0.003687  0.003624  0.003804   mean  0.004407
  Gibbs 10000 mean  0.001509  0.001758  0.001512  0.000846  0.001120   mean  0.001349
  Gibbs 10000 max   0.004425  0.005087  0.004483  0.003194  0.003860   mean  0.004210
J = 0.44:  LBP converged 5/5 (iterations 100, 99, 78, 82, 100),  mean field converged 5/5
  LBP   mean|err|  0.004251  0.004596  0.001670  0.001511  0.005706   mean  0.003547
  LBP   max|err|   0.012889  0.017231  0.004355  0.009327  0.020387   mean  0.012838
  Bethe logZ err  -2.88e-02 -2.53e-02 -9.25e-03 -1.50e-02 -3.46e-02   mean -2.26e-02
  MF    mean|err|  0.044734  0.035595  0.029855  0.033901  0.056003   mean  0.040018
  MF    max|err|   0.114534  0.099532  0.058054  0.101585  0.151487   mean  0.105038
  conf. exact      0.296527  0.278333  0.351118  0.337351  0.301923   mean  0.313050
  conf. MF         0.341261  0.313928  0.380974  0.367754  0.355361   mean  0.351855
  Gibbs   100 mean  0.032020  0.018893  0.014943  0.011268  0.016647   mean  0.018754
  Gibbs  1000 mean  0.012998  0.006322  0.005005  0.004217  0.005869   mean  0.006882
  Gibbs 10000 mean  0.001680  0.002159  0.001438  0.001464  0.001457   mean  0.001640
  Gibbs 10000 max   0.003597  0.006053  0.003272  0.003589  0.003713   mean  0.004045
J = 0.7:  LBP converged 5/5 (iterations 351, 163, 131, 136, 163),  mean field converged 5/5
  LBP   mean|err|  0.021321  0.046650  0.016015  0.006524  0.041109   mean  0.026324
  LBP   max|err|   0.067545  0.102417  0.045159  0.035217  0.107123   mean  0.071492
  Bethe logZ err  -2.69e-01 -1.33e-01 -7.56e-02 -9.44e-02 -1.47e-01   mean -1.44e-01
  MF    mean|err|  0.155674  0.144815  0.076454  0.095863  0.129514   mean  0.120464
  MF    max|err|   0.411025  0.335709  0.182949  0.350648  0.288592   mean  0.313785
  conf. exact      0.272050  0.250911  0.331390  0.338273  0.307481   mean  0.300021
  conf. MF         0.366833  0.380884  0.404906  0.426657  0.424471   mean  0.400750
  Gibbs   100 mean  0.052781  0.037318  0.024000  0.009220  0.021990   mean  0.029062
  Gibbs  1000 mean  0.010726  0.009297  0.013844  0.004597  0.010576   mean  0.009808
  Gibbs 10000 mean  0.002218  0.002795  0.001424  0.001887  0.001575   mean  0.001980
  Gibbs 10000 max   0.006861  0.006376  0.005037  0.008129  0.004064   mean  0.006093
J = 1.0:  LBP converged 5/5 (iterations 110, 122, 213, 154, 208),  mean field converged 5/5
  LBP   mean|err|  0.025633  0.182160  0.201876  0.035665  0.137108   mean  0.116488
  LBP   max|err|   0.053124  0.231727  0.336303  0.108630  0.278930   mean  0.201743
  Bethe logZ err  -5.77e-02 -2.79e-01 -4.58e-01 -1.27e-01 -3.61e-01   mean -2.56e-01
  MF    mean|err|  0.432465  0.248235  0.273394  0.253512  0.212425   mean  0.284006
  MF    max|err|   0.926632  0.352571  0.500774  0.802118  0.444833   mean  0.605385
  conf. exact      0.464544  0.223216  0.206207  0.412501  0.265179   mean  0.314329
  conf. MF         0.444524  0.448188  0.454047  0.457461  0.446337   mean  0.450111
  Gibbs   100 mean  0.021919  0.351203  0.116871  0.049219  0.010296   mean  0.109902
  Gibbs  1000 mean  0.001873  0.081031  0.049965  0.008149  0.009621   mean  0.030128
  Gibbs 10000 mean  0.001889  0.038433  0.005391  0.002066  0.004078   mean  0.010372
  Gibbs 10000 max   0.005420  0.046885  0.012468  0.005453  0.011061   mean  0.016257

H1a LBP mean|err| larger at J=1.0 than at J=0.1: 5/5 seeds
H1b seed-mean LBP mean|err| non-decreasing in J: True
H2  |Bethe log Z error| larger at J=1.0 than at J=0.1: 5/5 seeds
H3  mean-field confidence above exact at J=0.7 and 1.0: 9/10 (J, seed) pairs
H4  Gibbs seed-mean error falls 100 > 1000 > 10000 sweeps: 5/5 couplings

C7 before the pixel is seen: P(dark) 0.5000, P(dark | shadow) 0.5000
after I = 0.2: P(dark | I) 0.5028, P(dark | I, shadow) 0.0110, P(dark | I, lit) 1.0000
H5 observing the shadow lowers P(dark | I): True

All seven checks held. The data have the planned 24 and 15 edges, 65,536 distinct states and 7 and 6 foreground pixels (C1). At $J = 0$ enumeration matches the closed form to $1.1 \times 10^{-14}$ (C2), and loopy BP, mean field and the Bethe $\log Z$ match it to $1.3 \times 10^{-14}$ (C5): without coupling every pixel is on its own and all three methods are exact by construction. Gibbs is within 0.00571 of the exact marginals after 10,000 sweeps (C4). On the chain, loopy BP converged in all 25 runs and its marginals and Bethe $\log Z$ match enumeration to $1.2 \times 10^{-14}$ and $1.4 \times 10^{-14}$ (C3): on a tree, BP is exact and so is the Bethe free energy, as Yedidia et al. state. The two runs printed identical output (C6), and before the pixel is seen $P(\text{dark} \mid \text{shadow})$ equals $P(\text{dark})$, 0.5000 (C7), as d-separation predicts for two parents of an unobserved child.

Loopy BP. H1 held. H1a held on all five seeds: on each, loopy BP’s mean absolute error at $J = 1.0$ is above its value at $J = 0.1$. H1b held: the seed means rise at every step, from 0.000009 at $J = 0.1$ to 0.116488 at $J = 1.0$. Loopy BP converged on all 25 grid runs, in 48 to 351 iterations, so its errors are those of its fixed points, not of unfinished runs. At every coupling and seed its mean error is below mean field’s.

Bethe log Z. H2 held: on every seed the magnitude of the Bethe error at $J = 1.0$ (from $5.77 \times 10^{-2}$ to $4.58 \times 10^{-1}$) is above its magnitude at $J = 0.1$ (from $3.06 \times 10^{-5}$ to $1.11 \times 10^{-4}$). On every printed grid run the Bethe estimate is below the exact $\log Z$; the sources read make no claim about its sign, and the example offers no reason for it.

Mean field. H3 failed: it held in nine of its ten (coupling, seed) cases and failed in one. Mean field’s confidence exceeds the exact one on all five seeds at $J = 0.7$ and on seeds 1–4 at $J = 1.0$; on seed 0 at $J = 1.0$ it is below, 0.444524 against the exact 0.464544. Mean field’s mean error is above loopy BP’s on every seed at every coupling, with seed means from 0.002141 at $J = 0.1$ to 0.284006 at $J = 1.0$, and its largest single-pixel error reaches 0.926632 (seed 0, $J = 1.0$).

Gibbs. H4 held: at every coupling the seed-mean error falls from 100 to 1000 to 10,000 sweeps, from 0.011721 to 0.004991 to 0.001380 at $J = 0.1$ and from 0.109902 to 0.030128 to 0.010372 at $J = 1.0$. With 10,000 sweeps Gibbs is less accurate than loopy BP on every seed at $J = 0.1$ and $J = 0.25$, and more accurate on every seed at $J = 0.44$ and above. It converges more slowly on some runs than others: seed 1 at $J = 1.0$ still has a mean error of 0.038433 after 10,000 sweeps, against 0.001889–0.005391 on the other seeds at that coupling.

Explaining away. H5 held. After the grey pixel is seen, the surface is dark with probability 0.5028; learning that it is in shadow lowers this to 0.0110, and learning that it is lit raises it to 1.0000. By construction, dark paint in shadow predicts an intensity of 0.05, three noise standard deviations below the observed 0.2, while light paint in shadow predicts 0.2 exactly; in light, dark paint predicts 0.2 and light paint 0.8. Reflectance and illumination were independent before the observation (C7) and are strongly dependent after it, the collider rule of d-separation in numbers.

Process note. The plan listed eleven departures from the batch plan’s description of this example, before the pilot: (1) the data were fixed (the truth patch and chain, intensities 0 and 1, $\sigma = 0.6$, seeds 0–4); (2) the couplings were fixed at 0.1, 0.25, 0.44, 0.7 and 1.0, with a $J = 0$ base block in which every method must reach the exact answer; (3) loopy BP’s settings were fixed (synchronous, damping 0.5, uniform start, tolerance $10^{-14}$, at most 1000 iterations); (4) mean field uses damped parallel updates from $\mu = 0$, as Murphy’s Ising example; (5) Gibbs uses four chains, 100 discarded sweeps and nested estimates at 100, 1000 and 10,000 sweeps; (6) the “confidence” measure was added, because the batch plan’s over-confidence hypothesis had no measure; (7) the chain check also compares the Bethe $\log Z$ with the exact one; (8) the Bethe free energy is taken from Yedidia, Freeman and Weiss (2000, Eq. 4), because Murphy’s main text treats it only in a supplement we did not read; (9) checks C1–C7 were added; (10) the optional directed block was included, with C7 and H5; (11) the batch plan’s hypotheses were made precise as H1–H5, with loopy BP’s growth of error tested both per seed and as a trend of the seed means. The tolerances of four checks were set after a planning script had printed the reachable values they were set against: C2, C3 and C5 ($10^{-12}$ each), after that script printed floating-point agreement of about $10^{-14}$ between enumeration and closed forms or an independent forward–backward pass; and C4 (0.015), after it printed 0.00797 as the largest error of a direct simulation of 40,000 independent draws per pixel. The counts in C1 (24 and 15 edges) and the directed block’s numbers (reflectances 0.2 and 0.8, illuminations 0.25 and 1.0, noise 0.05, observation 0.2) and the Gibbs seeds were first written after that script ran; the script printed nothing about the directed block. It ran no loopy BP, mean field or Gibbs sampler and printed none of the measures H1–H5 compare, so H1–H5 are tests, not illustrations. The pilot printed the same output as the first run; no change was made after it.

Common Misconceptions

• “Bayesian networks are Bayesian.” There is nothing inherently Bayesian about them; they are just a way of defining probability distributions (Murphy, 2023, §4.2). A Bayesian neural network, a network with a posterior over its weights, is a different thing.
• “An edge means dependence, a missing edge means independence.” A missing edge means conditional independence given other variables, and in a directed graph conditioning can also create dependence: two causes with no edge between them become dependent once their common effect is observed (explaining away; C7 and H5 in the example).
• “Directed and undirected models are interchangeable.” Each can express independences the other cannot (the v-structure, the 4-cycle); moralizing a directed graph loses its v-structures (Murphy, 2023, §4.5).
• “Belief propagation is exact.” Only on trees and chains (C3 in the example). On graphs with loops it is an approximation whose fixed points are stationary points of the Bethe free energy (Yedidia et al., 2000), and it need not converge (Murphy, 2023, §9.4). In the example its error grew from about $10^{-5}$ to about $10^{-1}$ as the coupling rose.
• “Sampling is exact if run long enough, so it is always the better choice.” Its error falls only as $1/\sqrt{S}$ for independent samples, and successive Gibbs samples are correlated; in the example loopy BP was more accurate than 10,000 Gibbs sweeps at weak coupling, and Gibbs more accurate at strong coupling.
• “The graph of a neural network is a graphical model.” Torralba et al. (2024, §29.8) contrast the two: a PGM’s nodes are random variables and its edges statistical dependencies, while a network’s nodes are activations and its edges terms of weighted sums; PGMs compute probabilities, networks minimize a loss. In our reading, computational graphs, scene graphs and graph neural networks are not graphical models in this sense either.

Where It Is Used

• Pixel and label models. Markov random fields for restoration, segmentation and stereo, minimized by graph cuts, loopy BP and other optimizers, and in semi-global matching along scanlines. Conditional random fields condition the potentials on the image, and DeepLab v1 and v2 refine network outputs with a fully connected one.
• Estimation in geometry. Factor graphs for SLAM, pose graph optimization and bundle adjustment, where elimination is sparse matrix factorization.
• Time series and tracking. State-space models and their recursive inference: the Bayes filter, Kalman filter and particle filter. Hidden Markov models in handwriting recognition and in early temporal action segmentation, and forward–backward recursions in connectionist temporal classification.
• Part-based models. Pictorial structures and deformable part models are tree-structured models matched exactly by max-product dynamic programming.
• Latent-variable models. The Gaussian mixture model and others fitted with expectation–maximization; deep latent-variable models such as the variational autoencoder, fitted through the evidence lower bound.
• Bayesian reasoning about images. Bayesian inference and image formation as inference, in which explanations such as reflectance and illumination compete for the same evidence.