"""Bishop (2006), Chapter 8: Graphical Models. Reproduces the figures of Chapter 8 of C. M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006), pp. 359-422, written from scratch with numpy, scipy and matplotlib. Most figures of this chapter are diagrams; we draw them with the small helper class Diagram below and, wherever a diagram makes a numerical claim, compute the numbers with the algorithms of this file and print them on the figure. Every algorithm is a short function meant to be read next to the lecture notes (ch08-graphical-models.html): * Section 8.1, Bayesian networks with discrete variables: joint_table (the factorisation (8.5)), ancestral_sample (Section 8.1.2), random_cpts; linear-Gaussian models: linear_gaussian_moments (the recursions (8.15) and (8.16)) and linear_gaussian_sample (8.14). * Section 8.2, conditional independence: cond_mutual_info (a numerical test of a | b given c on a joint table), d_separated (the reachability or "Bayes-ball" rule of Section 8.2.2), trails (every path between two nodes, each with the node that blocks it and why), markov_blanket (Section 8.2.2, Figure 8.26), fuel_gauge (the numbers (8.30) to (8.33)). * Section 8.3, Markov random fields: moralize (Section 8.3.4), ising_energy (8.42), icm (iterated conditional modes), gibbs_ising (Gibbs sampling at temperature 1, for the demo) and graph_cut_map (the exact MAP of (8.42) by a minimum s-t cut, Greig, Porteous and Seheult 1989), with make_binary_image (our own test image) and flip_noise. * Section 8.4, inference: chain_sum_product ((8.54) to (8.58)), FactorGraph with brute-force enumeration, sum_product ((8.61) to (8.72), observed variables by indicator functions as in Section 8.4.4), max_sum ((8.87) to (8.98) with back-tracking on a tree), chain_max_sum ((8.99) to (8.102), with the erratum: mu_{x_{n-1} -> f_{n-1,n}}(x_{n-1}) inside (8.101)), loopy_bp (the flooding schedule of Section 8.4.7). Figures: 8.2, 8.6, 8.15 (with 8.16 to 8.20), 8.21, 8.22, 8.26, 8.28, 8.30, 8.32 (with 8.33), 8.40, 8.51, 8.52, 8.53. Run `python3 ch08_graphical_models.py --list`, or with no option for all figures (about ten seconds). Files go to ../img/ch08/ as PNG, SVG and PDF. What differs from the book, and the choices made: * Figure 8.30: the book's image (hand-written "Bayes' Theorem") is not available. We render our own 150 x 200 binary image (the words "Graphical models" in a rounded box, drawn with matplotlib's DejaVu Sans Bold font; 21.2% of the pixels are ink), flip every pixel independently with probability 0.1 (numpy seed 8: 9.88% flipped), and use the book's beta = 1.0, eta = 2.1, h = 0. ICM visits the pixels in raster order, starting from x = y, until a whole sweep changes nothing (here the second sweep). Pixels that differ from the original: ICM 4.52% (book: 4%), graph cut 2.05% (book: 1%). Energies (8.42): noisy image -98150, ICM -104208, graph cut -105351 (the global minimum; the original image has -104821, so the exact MAP is not the original: the remaining 2% are corners and thin strokes that the Ising prior smooths away). The graph cut is checked against brute force on small images in checks_ch08.py. * Figures 8.15 to 8.22: the probabilities in the three 3-node graphs and in Figure 8.22 come from conditional probability tables drawn from a Dirichlet(1) distribution (numpy seed 3); the numbers written on the figures are conditional mutual informations, which are zero exactly when the independence holds. The fuel-gauge numbers are the book's: 0.315, 0.81, 0.257, 0.111. * Figure 8.52: binary variables with the integer factor tables FIG51_TABLES below, so that the messages are whole numbers; Figure 8.53: a chain of N = 5 variables with K = 3 states and the pairwise tables CHAIN_PSI below (the same chain as in the sum-product demo). """ import itertools import math import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle, FancyArrowPatch, FancyBboxPatch, Rectangle from scipy.sparse import csr_matrix from scipy.sparse.csgraph import breadth_first_order, maximum_flow import common NODE = common.RED # node outlines and links, as in the book OBS = "#8c9cf5" # fill of observed (shaded) nodes PLATE = common.BLUE GREY = "#b0b0b0" OK = common.GREEN # =========================================================================== shared examples # Figure 8.51: p(x) proportional to f_a(x1, x2) f_b(x2, x3) f_c(x2, x4), binary variables # (variables numbered 0..3 for x1..x4). Integer tables keep the messages whole numbers. FIG51_VARS = [(0, 1), (1, 2), (1, 3)] FIG51_TABLES = [np.array([[2.0, 1.0], [1.0, 3.0]]), # f_a(x1, x2) np.array([[1.0, 2.0], [3.0, 1.0]]), # f_b(x2, x3) np.array([[2.0, 1.0], [1.0, 1.0]])] # f_c(x2, x4) # Figure 8.53 and the sum-product demo: a chain x1 - x2 - x3 - x4 - x5 with K = 3 states and # pairwise potentials psi_{n,n+1}(x_n, x_{n+1}) (rows: x_n, columns: x_{n+1}). CHAIN_PSI = [np.array([[1.0, 0.4, 2.2], [2.0, 0.6, 0.3], [0.5, 1.5, 1.0]]), np.array([[0.3, 1.8, 0.7], [1.2, 0.5, 2.4], [0.9, 2.0, 0.4]]), np.array([[2.5, 0.6, 0.8], [0.4, 1.1, 1.6], [1.0, 0.3, 2.1]]), np.array([[0.7, 1.9, 0.5], [2.2, 0.4, 1.3], [0.6, 1.0, 1.7]])] # =========================================================================== drawing helper class Diagram: """Draw a graphical model on a matplotlib axes, in the conventions of Section 8.1. node(key, (x, y), label, observed=..., kind=...) adds a random variable (an open circle, shaded when observed), a deterministic parameter (kind="param": a small solid dot) or a factor (kind="factor": a small solid square). edge(a, b) adds a directed link a -> b (directed=False for an undirected or factor-graph link), plate(...) a rounded box with its replication count, and done() fixes the limits and hides the axes. """ LPOS = {"left": (-1, 0, "right", "center"), "right": (1, 0, "left", "center"), "above": (0, 1, "center", "bottom"), "below": (0, -1, "center", "top"), "upper left": (-0.72, 0.72, "right", "bottom"), "upper right": (0.72, 0.72, "left", "bottom"), "lower left": (-0.72, -0.72, "right", "top"), "lower right": (0.72, -0.72, "left", "top")} def __init__(self, ax, r=0.28, lw=2.0, fs=15, color=NODE): self.ax, self.r, self.lw, self.fs, self.color = ax, r, lw, fs, color self.nodes = {} self.box = [np.inf, np.inf, -np.inf, -np.inf] def _grow(self, x0, y0, x1, y1): b = self.box self.box = [min(b[0], x0), min(b[1], y0), max(b[2], x1), max(b[3], y1)] def node(self, key, xy, label=None, observed=False, kind="var", lpos="left", color=None, fill=None, ls="-", fs=None, r=None): x, y = xy col = color or self.color if kind == "var": rad = r or self.r face = fill if fill is not None else (OBS if observed else "white") self.ax.add_patch(Circle(xy, rad, fc=face, ec=col, lw=self.lw, ls=ls, zorder=3)) elif kind == "param": rad = 0.06 self.ax.add_patch(Circle(xy, rad, fc=col, ec=col, lw=0, zorder=3)) elif kind == "point": # an invisible anchor for links that end in dots rad = 0.0 else: # factor rad = 0.12 self.ax.add_patch(Rectangle((x - rad, y - rad), 2 * rad, 2 * rad, fc=col, ec=col, lw=0, zorder=3)) self.nodes[key] = (np.array(xy, float), rad, kind) self._grow(x - rad, y - rad, x + rad, y + rad) if label: dx, dy, ha, va = self.LPOS[lpos] if isinstance(lpos, str) else lpos off = rad + (0.07 if kind == "var" else 0.08) tx, ty = x + dx * off, y + dy * off self.ax.text(tx, ty, label, ha=ha, va=va, fontsize=fs or self.fs, zorder=4) w = 0.16 * len(label.replace("$", "").replace("\\", "").replace("_", "")[:4]) + 0.1 self._grow(tx - (w if ha != "left" else 0), ty - 0.3, tx + (w if ha != "right" else 0), ty + 0.3) def _cut(self, key, u): """Distance from the centre of node key to its boundary in the direction u.""" _, rad, kind = self.nodes[key] if kind == "factor": return rad / max(abs(u[0]), abs(u[1])) return rad def edge(self, a, b, directed=True, color=None, lw=None, ls="-", gap=0.0): pa, pb = self.nodes[a][0], self.nodes[b][0] u = (pb - pa) / np.linalg.norm(pb - pa) p0 = pa + u * (self._cut(a, u) + gap) p1 = pb - u * (self._cut(b, u) + gap) col = color or self.color if directed: self.ax.add_patch(FancyArrowPatch(p0, p1, arrowstyle="-|>", mutation_scale=15, lw=lw or self.lw, color=col, ls=ls, shrinkA=0, shrinkB=0, zorder=2)) else: self.ax.plot([p0[0], p1[0]], [p0[1], p1[1]], color=col, lw=lw or self.lw, ls=ls, zorder=2, solid_capstyle="butt") def plate(self, x0, y0, x1, y1, label, color=PLATE): self.ax.add_patch(FancyBboxPatch((x0, y0), x1 - x0, y1 - y0, boxstyle="round,pad=0,rounding_size=0.15", fc="none", ec=color, lw=self.lw, zorder=1)) self.ax.text(x1 - 0.12, y0 + 0.1, label, ha="right", va="bottom", fontsize=self.fs) self._grow(x0, y0, x1, y1) def text(self, x, y, s, **kw): kw.setdefault("fontsize", 10) kw.setdefault("ha", "center") kw.setdefault("va", "center") self.ax.text(x, y, s, **kw) self._grow(x, y, x, y) def done(self, pad=0.15): b = self.box self.ax.set_xlim(b[0] - pad, b[2] + pad) self.ax.set_ylim(b[1] - pad, b[3] + pad) self.ax.set_aspect("equal") self.ax.axis("off") def draw_dag(ax, pos, parents, observed=(), labels=None, lpos=None, **kw): """Draw a directed graph given node positions and a parents dict; returns the Diagram.""" d = Diagram(ax, **kw) for v, xy in pos.items(): d.node(v, xy, (labels or {}).get(v, "$%s$" % v), observed=v in observed, lpos=(lpos or {}).get(v, "left")) for v, ps in parents.items(): for p in ps: d.edge(p, v) return d # =========================================================================== Section 8.1: discrete Bayesian networks def topological_order(parents): """Nodes ordered so that every parent comes before its children (Exercise 8.2: a DAG has one).""" order, done = [], set() while len(order) < len(parents): ready = [v for v in parents if v not in done and all(p in done for p in parents[v])] if not ready: raise ValueError("the graph has a directed cycle") for v in ready: order.append(v) done.add(v) return order def children_of(parents): ch = {v: [] for v in parents} for v, ps in parents.items(): for p in ps: ch[p].append(v) return ch def random_cpts(parents, card, rng, alpha=1.0): """Conditional probability tables p(x_k | pa_k) with rows drawn from Dirichlet(alpha). cpt[k] has one axis per parent (in the order of parents[k]) and a last axis for x_k. """ cpts = {} for v, ps in parents.items(): shape = tuple(card[p] for p in ps) cpts[v] = rng.dirichlet(alpha * np.ones(card[v]), size=shape) if shape else rng.dirichlet(alpha * np.ones(card[v])) return cpts def joint_table(order, parents, cpts): """The joint distribution p(x) = prod_k p(x_k | pa_k) (8.5) as an array with one axis per node of order.""" letters = {v: chr(ord("a") + i) for i, v in enumerate(order)} terms, ops = [], [] for v in order: terms.append("".join(letters[p] for p in parents[v]) + letters[v]) ops.append(cpts[v]) return np.einsum(",".join(terms) + "->" + "".join(letters[v] for v in order), *ops) def ancestral_sample(parents, cpts, rng, n): """Ancestral sampling (Section 8.1.2): draw each node in topological order from p(x_k | pa_k). Returns a dict node -> array of n integer states. """ out = {} for v in topological_order(parents): table = cpts[v] rows = table[tuple(out[p] for p in parents[v])] if parents[v] else np.broadcast_to(table, (n, table.shape[-1])) cum = np.cumsum(rows, axis=-1) u = rng.random(n)[:, None] out[v] = np.minimum((u > cum).sum(-1), table.shape[-1] - 1) return out def cond_mutual_info(P, order, A, B, C=()): """Conditional mutual information I(A; B | C) in nats from a joint table P over the nodes order. It is zero exactly when A and B are conditionally independent given C (Section 8.2), so it is the numerical test we compare d-separation with. """ ax = {v: i for i, v in enumerate(order)} keep = [ax[v] for v in list(A) + list(B) + list(C)] other = tuple(i for i in range(len(order)) if i not in keep) Q = P.sum(axis=other) if other else P # reorder the axes of Q as A, B, C cur = sorted(keep) Q = np.transpose(Q, [cur.index(k) for k in keep]) nA, nB = len(A), len(B) pABC = Q pAC = Q.sum(axis=tuple(range(nA, nA + nB)), keepdims=True) pBC = Q.sum(axis=tuple(range(nA)), keepdims=True) pC = Q.sum(axis=tuple(range(nA + nB)), keepdims=True) with np.errstate(divide="ignore", invalid="ignore"): t = pABC * (np.log(pABC) + np.log(pC) - np.log(pAC) - np.log(pBC)) return float(np.nansum(np.where(pABC > 0, t, 0.0))) # ----------------------------------------------------------------------------- Section 8.2.2: d-separation def ancestors_of(parents, S): """S together with every ancestor of a node of S.""" out, stack = set(), list(S) while stack: v = stack.pop() if v not in out: out.add(v) stack.extend(parents[v]) return out def d_separated(parents, A, B, C=()): """True when the node sets A and B are d-separated by C (Section 8.2.2). Reachability ("Bayes-ball") form of the rule: walk along links from A, remembering whether we arrived at a node from a child ("up") or from a parent ("down"). A path may continue * through an unobserved node in any direction except down-then-up (a head-to-tail or tail-to-tail node that is not in C does not block, and it is not a collider), and * from down to up (a head-to-head node, the collider) only when the node or one of its descendants is in C, which are exactly the blocking conditions (i) and (ii) of Section 8.2.2. B is d-separated from A when no node of B can be reached. """ A, B, C = set(A), set(B), set(C) ch = children_of(parents) anc_C = ancestors_of(parents, C) # nodes with a descendant (or themselves) in C seen, reach = set(), set() stack = [(a, "up") for a in A] while stack: v, d = stack.pop() if (v, d) in seen: continue seen.add((v, d)) if v not in C: reach.add(v) if d == "up" and v not in C: # came from a child: go on to parents and children stack.extend((p, "up") for p in parents[v]) stack.extend((c, "down") for c in ch[v]) elif d == "down": # came from a parent if v not in C: # head-to-tail through v stack.extend((c, "down") for c in ch[v]) if v in anc_C: # head-to-head at v, unblocked by C stack.extend((p, "up") for p in parents[v]) return not (reach & (B - C)) def trails(parents, a, b, C=()): """Every simple path between a and b in the graph with directions ignored. Returns a list of (path, blockers): blockers lists the nodes (node, reason) at which C blocks the path, with reason 'tail-to-tail, in C', 'head-to-tail, in C' or 'head-to-head, neither it nor a descendant in C' (Section 8.2.2); an empty list means the path is open. a and b are d-separated by C exactly when every path has at least one blocker. """ C = set(C) ch = children_of(parents) nb = {v: set(parents[v]) | set(ch[v]) for v in parents} desc_ok = ancestors_of(parents, C) out = [] def walk(path): v = path[-1] if v == b: out.append(list(path)) return for w in sorted(nb[v], key=str): if w not in path: walk(path + [w]) walk([a]) res = [] for p in out: blockers = [] for i in range(1, len(p) - 1): u, m, w = p[i - 1], p[i], p[i + 1] if u in parents[m] and w in parents[m]: if m not in desc_ok: blockers.append((m, "head-to-head, neither it nor a descendant in C")) elif m in C: kind = "tail-to-tail" if (m in parents[u] and m in parents[w]) else "head-to-tail" blockers.append((m, kind + ", in C")) res.append((p, blockers)) return res def markov_blanket(parents, v): """Parents, children and co-parents of v (Figure 8.26).""" ch = children_of(parents) mb = set(parents[v]) | set(ch[v]) for c in ch[v]: mb |= set(parents[c]) mb.discard(v) return mb def fuel_gauge(): """The fuel-gauge example of Section 8.2.1 ((8.30) to (8.33)): returns the four probabilities. B battery (1 = charged), F fuel tank (1 = full), G gauge (1 = reads full). """ pB = np.array([0.1, 0.9]) pF = np.array([0.1, 0.9]) pG1 = np.array([[0.1, 0.2], [0.2, 0.8]]) # p(G = 1 | B, F), indexed [B, F] pG = np.stack([1 - pG1, pG1], -1) # p(G | B, F), axes B, F, G P = pB[:, None, None] * pF[None, :, None] * pG # joint, axes B, F, G p_g0 = P[:, :, 0].sum() # (8.30) p_g0_f0 = (pG[:, 0, 0] * pB).sum() # (8.31) p_f0_g0 = P[:, 0, 0].sum() / p_g0 # (8.32) p_f0_g0_b0 = P[0, 0, 0] / P[0, :, 0].sum() # (8.33) return {"p(G=0)": p_g0, "p(G=0|F=0)": p_g0_f0, "p(F=0|G=0)": p_f0_g0, "p(F=0|G=0,B=0)": p_f0_g0_b0, "p(F=0)": pF[0], "joint": P} # ----------------------------------------------------------------------------- Section 8.1.4: linear-Gaussian models def linear_gaussian_moments(W, b, v): """Mean and covariance of a linear-Gaussian DAG by the recursions (8.15) and (8.16). Node i has p(x_i | pa_i) = N(x_i | sum_j W[i, j] x_j + b_i, v_i) (8.11), with W strictly lower triangular (the nodes are numbered in topological order). E[x_i] = sum_j w_ij E[x_j] + b_i (8.15) and cov[x_i, x_j] = sum_k w_jk cov[x_i, x_k] + I_ij v_j (8.16), filled in for j = 1, ..., D and i <= j, so that every term on the right is already known. """ D = len(b) mu = np.zeros(D) S = np.zeros((D, D)) for j in range(D): mu[j] = W[j, :j] @ mu[:j] + b[j] for i in range(j + 1): S[i, j] = sum(W[j, k] * S[i, k] for k in range(j)) + (v[j] if i == j else 0.0) S[j, i] = S[i, j] return mu, S def linear_gaussian_sample(W, b, v, rng, n): """Ancestral sampling of (8.14): x_i = sum_j w_ij x_j + b_i + sqrt(v_i) eps_i.""" D = len(b) X = np.zeros((n, D)) for i in range(D): X[:, i] = X[:, :i] @ W[i, :i] + b[i] + np.sqrt(v[i]) * rng.standard_normal(n) return X # =========================================================================== Section 8.3: Markov random fields def moralize(parents): """Moral graph of a DAG (Section 8.3.4): drop directions and marry the parents of every node. Returns the set of undirected links as frozensets {u, v}. """ E = set() for v, ps in parents.items(): for p in ps: E.add(frozenset((p, v))) for p, q in itertools.combinations(ps, 2): E.add(frozenset((p, q))) return E def ising_energy(x, y, h=0.0, beta=1.0, eta=2.1): """E(x, y) = h sum_i x_i - beta sum_{i~j} x_i x_j - eta sum_i x_i y_i (8.42), 4-neighbour grid.""" x = np.asarray(x, float) pair = (x[1:, :] * x[:-1, :]).sum() + (x[:, 1:] * x[:, :-1]).sum() return float(h * x.sum() - beta * pair - eta * (x * y).sum()) def _nbr_sum(x, i, j): H, W = x.shape s = 0.0 if i > 0: s += x[i - 1, j] if i < H - 1: s += x[i + 1, j] if j > 0: s += x[i, j - 1] if j < W - 1: s += x[i, j + 1] return s def icm(y, h=0.0, beta=1.0, eta=2.1, max_sweeps=50, x0=None): """Iterated conditional modes (Kittler and Foglein 1984; Section 8.3.3) for the energy (8.42). Start from x = y (or x0), visit the pixels in raster order and set each x_j to the state of lower energy with all other pixels fixed: the terms of (8.42) that involve x_j are x_j (h - beta sum_{i~j} x_i - eta y_j), so x_j = +1 when beta sum x_i + eta y_j - h > 0, and x_j keeps its value on a tie. Stop after a sweep with no change (a local minimum). Returns x and the energy after 0, 1, 2, ... sweeps (never increasing, Exercise 8.13). """ x = (np.array(y, float) if x0 is None else np.array(x0, float)).copy() H, W = x.shape trace = [ising_energy(x, y, h, beta, eta)] for _ in range(max_sweeps): changed = 0 for i in range(H): for j in range(W): a = beta * _nbr_sum(x, i, j) + eta * y[i, j] - h new = 1.0 if a > 0 else (-1.0 if a < 0 else x[i, j]) if new != x[i, j]: x[i, j] = new changed += 1 trace.append(ising_energy(x, y, h, beta, eta)) if changed == 0: break return x, np.array(trace) def gibbs_ising(y, rng, h=0.0, beta=1.0, eta=2.1, sweeps=20, x0=None): """Gibbs sampling from p(x | y) proportional to exp(-E(x, y)) (8.43), raster order. p(x_j = +1 | rest) = sigma(2 (beta sum_{i~j} x_i + eta y_j - h)). Returns the last sample and the average of the samples (an estimate of E[x_j | y]). """ x = (np.array(y, float) if x0 is None else np.array(x0, float)).copy() H, W = x.shape mean = np.zeros_like(x) for _ in range(sweeps): u = rng.random((H, W)) for i in range(H): for j in range(W): a = beta * _nbr_sum(x, i, j) + eta * y[i, j] - h x[i, j] = 1.0 if u[i, j] < common.sigmoid(2 * a) else -1.0 mean += x / sweeps return x, mean def graph_cut_map(y, h=0.0, beta=1.0, eta=2.1, scale=1000): """Exact minimiser of the Ising energy (8.42) by a minimum s-t cut (Greig et al. 1989). With labels x_i in {-1, +1} the energy is, up to a constant, sum_i u_i(x_i) + sum_{i~j} 2 beta [x_i != x_j], u_i(x) = x (h - eta y_i), because -beta x_i x_j = 2 beta [x_i != x_j] - beta. Build a graph with a source s (label +1) and a sink t (label -1): an edge i -> t of capacity u_i(+1) - m_i (cut when i is on the source side), an edge s -> i of capacity u_i(-1) - m_i (cut when i is on the sink side), m_i = min(u_i(+1), u_i(-1)), and edges i -> j, j -> i of capacity 2 beta between neighbours. Every s-t cut is a labelling and its capacity is the energy minus a constant, so the minimum cut (= maximum flow) is the global minimum when beta >= 0. Capacities are multiplied by `scale` and rounded, which is exact for parameters with at most three decimals. Returns x (+1 on the source side of the minimum cut). """ y = np.asarray(y, float) H, W = y.shape N = H * W s, t = N, N + 1 idx = np.arange(N).reshape(H, W) a = (h - eta * y).ravel() u_plus, u_minus = a, -a m = np.minimum(u_plus, u_minus) cap_t = np.rint((u_plus - m) * scale).astype(np.int64) cap_s = np.rint((u_minus - m) * scale).astype(np.int64) pw = int(round(2 * beta * scale)) rows, cols, caps = [], [], [] for p, q in ((idx[:, :-1], idx[:, 1:]), (idx[:-1, :], idx[1:, :])): p, q = p.ravel(), q.ravel() rows += [p, q] cols += [q, p] caps += [np.full(p.size, pw), np.full(p.size, pw)] rows += [np.full(N, s), np.arange(N)] cols += [np.arange(N), np.full(N, t)] caps += [cap_s, cap_t] r, c, k = np.concatenate(rows), np.concatenate(cols), np.concatenate(caps) keep = k > 0 G = csr_matrix((k[keep].astype(np.int32), (r[keep], c[keep])), shape=(N + 2, N + 2)) res = maximum_flow(G, s, t) R = (G - res.flow).tocsr() R.data = (R.data > 0).astype(np.int8) R.eliminate_zeros() reach = breadth_first_order(R, s, directed=True, return_predecessors=False) x = -np.ones(N + 2) x[reach] = 1.0 return x[:N].reshape(H, W) def brute_force_map(y, h=0.0, beta=1.0, eta=2.1): """Global minimiser of (8.42) by trying all 2^D images (only for tiny images).""" H, W = y.shape best, bx = np.inf, None for bits in itertools.product((-1.0, 1.0), repeat=H * W): x = np.array(bits).reshape(H, W) e = ising_energy(x, y, h, beta, eta) if e < best - 1e-12: best, bx = e, x return bx, best def make_binary_image(H=150, W=200): """Our test image for Figure 8.30: 'Graphical models' in a rounded box, as a +-1 array. Rendered with matplotlib (DejaVu Sans Bold, which ships with matplotlib) and thresholded; +1 is the ink (shown in blue, as in the book) and -1 the background. """ from matplotlib.backends.backend_agg import FigureCanvasAgg from matplotlib.figure import Figure fig = Figure(figsize=(W / 100, H / 100), dpi=100) canvas = FigureCanvasAgg(fig) ax = fig.add_axes([0, 0, 1, 1]) ax.set_xlim(0, W) ax.set_ylim(0, H) ax.axis("off") ax.add_patch(FancyBboxPatch((9, 9), W - 18, H - 18, boxstyle="round,pad=0,rounding_size=16", fc="none", ec="black", lw=4.2)) ax.text(W / 2, H * 0.64, "Graphical", ha="center", va="center", fontsize=21, fontweight="bold", family="DejaVu Sans") ax.text(W / 2, H * 0.33, "models", ha="center", va="center", fontsize=21, fontweight="bold", family="DejaVu Sans") canvas.draw() img = np.asarray(canvas.buffer_rgba())[..., 0].astype(float) return np.where(img < 128, 1.0, -1.0) def flip_noise(x, p, rng): """Flip the sign of every pixel independently with probability p (Section 8.3.3).""" return np.where(rng.random(x.shape) < p, -x, x) def image_to_hex(x): """Pack a +-1 image into a hex string, row by row, 1 bit per pixel (most significant first). This is how the demo image-denoising-demo.html stores the same image. """ bits = (np.asarray(x).ravel() > 0).astype(np.uint8) return np.packbits(bits).tobytes().hex() def fig30_setup(seed=8, p=0.1): x = make_binary_image() y = flip_noise(x, p, np.random.default_rng(seed)) return x, y # =========================================================================== Section 8.4: inference def chain_sum_product(psis): """Marginals of a chain p(x) = (1/Z) prod_n psi_{n,n+1}(x_n, x_{n+1}) (8.49) by (8.54) to (8.58). psis[n] is the K_n x K_{n+1} table of psi_{n,n+1}. The forward messages mu_alpha(x_n) = sum_{x_{n-1}} psi_{n-1,n}(x_{n-1}, x_n) mu_alpha(x_{n-1}) (8.55), mu_alpha(x_1) = 1, and backward messages mu_beta(x_n) = sum_{x_{n+1}} psi_{n,n+1}(x_n, x_{n+1}) mu_beta(x_{n+1}) (8.57, in the corrected index order of the errata), mu_beta(x_N) = 1, give p(x_n) = mu_alpha mu_beta / Z (8.54) and p(x_{n-1}, x_n) = mu_alpha(x_{n-1}) psi mu_beta(x_n) / Z (8.58). """ N = len(psis) + 1 alpha = [np.ones(psis[0].shape[0])] for n in range(1, N): alpha.append(psis[n - 1].T @ alpha[-1]) beta = [np.ones(psis[-1].shape[1])] for n in range(N - 2, -1, -1): beta.insert(0, psis[n] @ beta[0]) Z = float(alpha[-1].sum()) marg = [a * b / Z for a, b in zip(alpha, beta)] pair = [alpha[n][:, None] * psis[n] * beta[n + 1][None, :] / Z for n in range(N - 1)] return {"alpha": alpha, "beta": beta, "Z": Z, "marginals": marg, "pairwise": pair} class FactorGraph: """A factor graph p(x) = (1/Z) prod_s f_s(x_s) (8.59) over discrete variables 0, ..., n-1. card[i] is the number of states of variable i; add(vars, table) adds a factor whose table has one axis per variable in vars (in that order). """ def __init__(self, card): self.card = list(card) self.factors = [] def add(self, vars_, table): table = np.asarray(table, float) assert table.shape == tuple(self.card[v] for v in vars_) self.factors.append((tuple(vars_), table)) return len(self.factors) - 1 def var_nbrs(self, i): return [s for s, (vs, _) in enumerate(self.factors) if i in vs] def joint(self, evidence=None): """The unnormalised joint prod_s f_s(x_s), times indicators of the evidence, by enumeration.""" n = len(self.card) letters = [chr(ord("a") + i) for i in range(n)] terms = [",".join("".join(letters[v] for v in vs) for vs, _ in self.factors)] ops = [t for _, t in self.factors] P = np.einsum(terms[0] + "->" + "".join(letters), *ops) if n > 0 else None # variables in no factor get a constant factor 1 for i, ind in (evidence or {}).items(): mask = np.zeros(self.card[i]) mask[ind] = 1.0 shape = [1] * n shape[i] = self.card[i] P = P * mask.reshape(shape) return P def brute_marginals(self, evidence=None): P = self.joint(evidence) Z = P.sum() n = len(self.card) return [P.sum(axis=tuple(j for j in range(n) if j != i)) / Z for i in range(n)], float(Z) def is_tree(self): nodes = len(self.card) + len(self.factors) edges = sum(len(vs) for vs, _ in self.factors) return edges == nodes - 1 and self._connected() def _connected(self): seen, stack = {("v", 0)}, [("v", 0)] while stack: for w in self.nbrs(stack.pop()): if w not in seen: seen.add(w) stack.append(w) return len(seen) == len(self.card) + len(self.factors) def nbrs(self, node): kind, k = node if kind == "v": return [("f", s) for s in self.var_nbrs(k)] return [("v", v) for v in self.factors[k][0]] def tree_schedule(self, root=0): """Message order on a tree: leaves to the root ("collect"), then the root back out ("distribute").""" order, parent = [], {("v", root): None} stack = [("v", root)] while stack: u = stack.pop() order.append(u) for w in self.nbrs(u): if w not in parent: parent[w] = u stack.append(w) collect = [(u, parent[u]) for u in reversed(order) if parent[u] is not None] distribute = [(parent[u], u) for u in order if parent[u] is not None] return collect, distribute def _log(a): """ln a with ln 0 = -inf and no warning (zeros come from observed variables and zero table entries).""" with np.errstate(divide="ignore"): return np.log(a) def _evidence_vec(fg, i, evidence): e = np.ones(fg.card[i]) if evidence and i in evidence: e = np.zeros(fg.card[i]) e[evidence[i]] = 1.0 return e def _factor_to_var(fg, s, i, incoming, op=np.sum, log=False): """Message f_s -> x_i (8.66): multiply f_s by the messages from its other variables, sum them out. With op=np.max and log=True this is the max-sum message (8.93). """ vs, table = fg.factors[s] T = _log(table) if log else table.copy() for k, v in enumerate(vs): if v == i: continue shape = [1] * len(vs) shape[k] = fg.card[v] m = incoming[(("v", v), ("f", s))].reshape(shape) T = T + m if log else T * m k = vs.index(i) return op(np.moveaxis(T, k, 0).reshape(fg.card[i], -1), axis=1) def _var_to_factor(fg, i, s, incoming, evidence, log=False): """Message x_i -> f_s (8.69): the product of the messages from the other factors of x_i (1 at a leaf (8.70)), times the indicator of the evidence on x_i (Section 8.4.4).""" e = _evidence_vec(fg, i, evidence) m = _log(e) if log else e.copy() for t in fg.var_nbrs(i): if t != s: m = m + incoming[(("f", t), ("v", i))] if log else m * incoming[(("f", t), ("v", i))] return m def sum_product(fg, root=0, evidence=None): """The sum-product algorithm on a tree-structured factor graph ((8.61) to (8.72)). Messages go from the leaves to the root and back (fg.tree_schedule). Observed variables (evidence: dict variable -> observed state) are clamped by multiplying in indicator functions, so that Z becomes the probability of the evidence when the factors are normalised conditionals (Section 8.4.4). Returns the marginals p(x_i | evidence) (8.63), the factor marginals p(x_s | evidence) (8.72), Z and every message in the order sent. """ msgs, sent = {}, [] collect, distribute = fg.tree_schedule(root) for src, dst in collect + distribute: if src[0] == "v": msgs[(src, dst)] = _var_to_factor(fg, src[1], dst[1], msgs, evidence) else: msgs[(src, dst)] = _factor_to_var(fg, src[1], dst[1], msgs) sent.append((src, dst, msgs[(src, dst)].copy())) unnorm = [] for i in range(len(fg.card)): m = _evidence_vec(fg, i, evidence) for s in fg.var_nbrs(i): m = m * msgs[(("f", s), ("v", i))] unnorm.append(m) Z = float(unnorm[root].sum()) fmarg = [] for s, (vs, table) in enumerate(fg.factors): T = table.copy() for k, v in enumerate(vs): shape = [1] * len(vs) shape[k] = fg.card[v] T = T * msgs[(("v", v), ("f", s))].reshape(shape) fmarg.append(T / Z) return {"marginals": [u / u.sum() for u in unnorm], "Z": Z, "factor_marginals": fmarg, "messages": sent} def max_sum(fg, root=0, evidence=None): """The max-sum algorithm with back-tracking on a tree ((8.87) to (8.98)). Log-domain messages from the leaves to the root: factor to variable (8.93), with the arg max of the other variables stored for every value of the receiving variable, and variable to factor (8.94), starting from 0 at leaf variables (8.95) and ln f at leaf factors (8.96). At the root, ln p_max = max of the sum of incoming messages (8.97) and x_root = its arg max (8.98). Back-tracking then walks out from the root and reads the stored arg max of each factor given the value of the variable it sent its message to, which is (8.101) and (8.102) on a general tree. Returns (x_max, ln of the unnormalised max). """ msgs, back = {}, {} collect, _ = fg.tree_schedule(root) for src, dst in collect: if src[0] == "v": msgs[(src, dst)] = _var_to_factor(fg, src[1], dst[1], msgs, evidence, log=True) else: s, i = src[1], dst[1] vs, table = fg.factors[s] T = _log(table) for k, v in enumerate(vs): if v != i: shape = [1] * len(vs) shape[k] = fg.card[v] T = T + msgs[(("v", v), ("f", s))].reshape(shape) k = vs.index(i) flat = np.moveaxis(T, k, 0).reshape(fg.card[i], -1) msgs[(src, dst)] = flat.max(1) others = [v for v in vs if v != i] shp = [fg.card[v] for v in others] back[s] = (i, others, [np.unravel_index(a, shp) for a in flat.argmax(1)]) tot = _log(_evidence_vec(fg, root, evidence)) for s in fg.var_nbrs(root): tot = tot + msgs[(("f", s), ("v", root))] x = {root: int(np.argmax(tot))} stack = [root] while stack: i = stack.pop() for s in fg.var_nbrs(i): if s in back and back[s][0] == i: _, others, arg = back[s] for v, val in zip(others, arg[x[i]]): if v not in x: x[v] = int(val) stack.append(v) return [x[i] for i in range(len(fg.card))], float(tot.max()) def chain_max_sum(psis): """Max-sum on the chain of (8.49), root x_N, with back-tracking ((8.99) to (8.102)). mu(x_1) = 0 (8.99); mu(x_n) = max_{x_{n-1}} [ln psi_{n-1,n}(x_{n-1}, x_n) + mu(x_{n-1})] and the arg max phi(x_n) (8.101, with the erratum: the message is a function of x_{n-1}); x_N = arg max mu(x_N) (8.100) and x_{n-1} = phi(x_n) (8.102). Returns the path, the messages mu and back-pointers phi (phi[n][k] = best state of x_{n-1} when x_n = k), and ln max prod psi. """ mu = [np.zeros(psis[0].shape[0])] phi = [None] for P in psis: S = np.log(P) + mu[-1][:, None] phi.append(S.argmax(0)) mu.append(S.max(0)) path = [int(np.argmax(mu[-1]))] for n in range(len(psis), 0, -1): path.insert(0, int(phi[n][path[0]])) return {"path": path, "mu": mu, "phi": phi, "logmax": float(mu[-1].max())} def loopy_bp(fg, iters=50, evidence=None, damping=0.0, tol=0.0): """Loopy belief propagation with the flooding schedule (Section 8.4.7). All variable-to-factor messages start uniform; each iteration sends every factor-to-variable message (8.66) from the current variable messages, then every variable-to-factor message (8.69), all normalised to sum to one. Beliefs are the normalised products (8.63). On a tree the beliefs are exact after a number of iterations equal to the diameter; on a graph with loops they are an approximation (and may fail to converge). Returns the history of beliefs. """ msgs = {} for s, (vs, _) in enumerate(fg.factors): for v in vs: msgs[(("v", v), ("f", s))] = np.ones(fg.card[v]) / fg.card[v] history = [] for _ in range(iters): new = {} for s, (vs, _) in enumerate(fg.factors): for v in vs: m = _factor_to_var(fg, s, v, msgs) new[(("f", s), ("v", v))] = m / m.sum() if damping and history: for k in new: new[k] = (1 - damping) * new[k] + damping * msgs[k] msgs.update(new) for s, (vs, _) in enumerate(fg.factors): for v in vs: m = _var_to_factor(fg, v, s, msgs, evidence) msgs[(("v", v), ("f", s))] = m / m.sum() bel = [] for i in range(len(fg.card)): b = _evidence_vec(fg, i, evidence) for s in fg.var_nbrs(i): b = b * msgs[(("f", s), ("v", i))] bel.append(b / b.sum()) history.append(bel) if tol and len(history) > 1 and max(np.abs(a - b).max() for a, b in zip(history[-1], history[-2])) < tol: break return history def fig51_graph(): fg = FactorGraph([2, 2, 2, 2]) for vs, t in zip(FIG51_VARS, FIG51_TABLES): fg.add(vs, t) return fg def chain_graph(psis=CHAIN_PSI): fg = FactorGraph([psis[0].shape[0]] + [P.shape[1] for P in psis]) for n, P in enumerate(psis): fg.add((n, n + 1), P) return fg def cycle_graph(J=1.0, fields=(0.6, -0.2, 0.3, -0.5, 0.1), couplings=(1.0, 0.8, 1.2, 0.9, 1.1)): """A cycle x1 - x2 - ... - x5 - x1 of binary (+-1) variables, used for loopy BP: pairwise factors exp(J c_k x_k x_{k+1}) and unary factors exp(h_k x_k) (an Ising ring).""" n = len(fields) fg = FactorGraph([2] * n) s = np.array([-1.0, 1.0]) for k in range(n): fg.add((k,), np.exp(fields[k] * s)) for k in range(n): fg.add((k, (k + 1) % n), np.exp(J * couplings[k] * np.outer(s, s))) return fg # =========================================================================== figures def _cmi_example(kind, seed=8): """Random CPTs for the three 3-node graphs: returns I(a;b) and I(a;b|c).""" G = {"tail": {"c": (), "a": ("c",), "b": ("c",)}, "chain": {"a": (), "c": ("a",), "b": ("c",)}, "head": {"a": (), "b": (), "c": ("a", "b")}}[kind] order = topological_order(G) cp = random_cpts(G, {v: 2 for v in G}, np.random.default_rng(seed)) P = joint_table(order, G, cp) return cond_mutual_info(P, order, ["a"], ["b"]), cond_mutual_info(P, order, ["a"], ["b"], ["c"]) def _fmt_mi(v): return "0" if abs(v) < 1e-12 else "%.4f" % v def fig8_2(out_dir=None): """Figure 8.2: the directed acyclic graph of (8.4) on seven variables.""" fig, ax = plt.subplots(figsize=(3.4, 4.0)) pos = {"x1": (1.0, 3.0), "x2": (-0.3, 2.0), "x3": (2.3, 2.0), "x4": (0.35, 0.75), "x5": (1.65, 0.75), "x6": (0.0, -0.5), "x7": (1.95, -0.5)} par = {"x1": (), "x2": (), "x3": (), "x4": ("x1", "x2", "x3"), "x5": ("x1", "x3"), "x6": ("x4",), "x7": ("x4", "x5")} lab = {v: "$x_%s$" % v[1] for v in pos} lp = {"x1": "left", "x2": "above", "x3": "above", "x4": "left", "x5": "right", "x6": "left", "x7": "right"} d = draw_dag(ax, pos, par, labels=lab, lpos=lp) d.done() common.savefig(fig, 8, "fig8-2", out_dir) def fig8_6(out_dir=None): """Figure 8.6: Bayesian polynomial regression with a plate, observed t_n and deterministic parameters.""" fig, ax = plt.subplots(figsize=(4.4, 3.0)) d = Diagram(ax) d.node("t", (0, 0), "$t_n$", observed=True, lpos="lower left") d.node("x", (0, 1.2), "$x_n$", kind="param", lpos="above") d.node("w", (1.9, 0), r"$\mathbf{w}$", lpos="upper right") d.node("alpha", (1.9, 1.2), r"$\alpha$", kind="param", lpos="above") d.node("s2", (-1.6, 0), r"$\sigma^2$", kind="param", lpos="left") d.edge("x", "t") d.edge("w", "t") d.edge("alpha", "w") d.edge("s2", "t") d.plate(-0.9, -0.85, 0.9, 1.75, "$N$") d.done() common.savefig(fig, 8, "fig8-6", out_dir) THREE = {"tail": ({"a": (-0.8, 0), "b": (0.8, 0), "c": (0, 1.0)}, {"c": (), "a": ("c",), "b": ("c",)}), "chain": ({"a": (-1.1, 0.5), "c": (0, 0.5), "b": (1.1, 0.5)}, {"a": (), "c": ("a",), "b": ("c",)}), "head": ({"a": (-0.8, 1.0), "b": (0.8, 1.0), "c": (0, 0)}, {"a": (), "b": (), "c": ("a", "b")})} def fig8_15(out_dir=None): """Figures 8.15 to 8.20: the three 3-node graphs, without and with c observed, with the conditional mutual information of a and b computed from random tables.""" fig, axes = plt.subplots(2, 3, figsize=(9.6, 4.6)) titles = {"tail": "tail-to-tail at $c$\n(Figs 8.15, 8.16)", "chain": "head-to-tail at $c$\n(Figs 8.17, 8.18)", "head": "head-to-head at $c$\n(Figs 8.19, 8.20)"} lp = {"tail": {"a": "left", "b": "right", "c": "left"}, "chain": {"a": "above", "b": "above", "c": "above"}, "head": {"a": "left", "b": "right", "c": "left"}} for j, kind in enumerate(("tail", "chain", "head")): pos, par = THREE[kind] mi0, mi1 = _cmi_example(kind) for i, obs in enumerate(((), ("c",))): ax = axes[i, j] d = draw_dag(ax, pos, par, observed=obs, lpos=lp[kind]) sep = d_separated(par, ["a"], ["b"], obs) cond = r" \mid c" if obs else r" \mid \emptyset" mi = mi1 if obs else mi0 stmt = (r"$a \perp b%s$" % cond) if sep else (r"$a, b$ dependent" + (" given $c$" if obs else "")) d.text(0, -0.75, stmt, fontsize=12, color=OK if sep else NODE) d.text(0, -1.15, "I(a; b%s) = %s" % (" | c" if obs else "", _fmt_mi(mi)), fontsize=9, color="0.25") d.done(pad=0.1) ax.set_xlim(-1.75, 1.75) if i == 0: ax.set_title(titles[kind], fontsize=10.5) fig.text(0.005, 0.73, "c unobserved", rotation=90, va="center", fontsize=10) fig.text(0.005, 0.27, "c observed", rotation=90, va="center", fontsize=10) fig.tight_layout(rect=(0.02, 0, 1, 1)) common.savefig(fig, 8, "fig8-15", out_dir) def fig8_21(out_dir=None): """Figure 8.21: explaining away in the fuel-gauge network, with the probabilities (8.30) to (8.33).""" fg = fuel_gauge() fig, axes = plt.subplots(1, 3, figsize=(9.0, 3.0)) pos = {"B": (-0.75, 1), "F": (0.75, 1), "G": (0, 0)} par = {"B": (), "F": (), "G": ("B", "F")} obs = [(), ("G",), ("G", "B")] lines = [r"$p(F=0) = %.3f$" % fg["p(F=0)"], r"$p(F=0\mid G=0) = %.3f$" % fg["p(F=0|G=0)"], r"$p(F=0\mid G=0, B=0) = %.3f$" % fg["p(F=0|G=0,B=0)"]] subs = ["prior", "gauge reads empty", "gauge empty, battery flat"] for k, ax in enumerate(axes): d = draw_dag(ax, pos, par, observed=obs[k], labels={v: "$%s$" % v for v in pos}, lpos={"B": "left", "F": "right", "G": "left"}) d.text(0, -0.65, lines[k], fontsize=11) d.text(0, -1.0, subs[k], fontsize=9, color="0.35") d.done(pad=0.1) ax.set_xlim(-1.7, 1.7) ax.set_title("(%s)" % "abc"[k], fontsize=10) fig.tight_layout() common.savefig(fig, 8, "fig8-21", out_dir) FIG22_POS = {"a": (0, 2), "f": (2, 2), "e": (1, 1), "b": (3, 1), "c": (1, 0)} FIG22_PAR = {"a": (), "f": (), "e": ("a", "f"), "b": ("f",), "c": ("e",)} def fig22_cpts(q=0.9): """Tables for Figure 8.22 with clearly visible dependences: a and f fair coins, e = (a != f) with probability 0.9, b a noisy copy of f and c a noisy copy of e (each copied with probability q).""" copy = np.array([[q, 1 - q], [1 - q, q]]) e = np.zeros((2, 2, 2)) for a in range(2): for f in range(2): p1 = 0.9 if a != f else 0.1 e[a, f] = [1 - p1, p1] return {"a": np.array([0.5, 0.5]), "f": np.array([0.5, 0.5]), "e": e, "b": copy, "c": copy} def fig8_22(out_dir=None): """Figure 8.22: d-separation in a five-node graph; (a) c observed, (b) f observed. With the tables of fig22_cpts, I(a; b | c) = 0.1375 nats and I(a; b | f) = 0. """ order = topological_order(FIG22_PAR) P = joint_table(order, FIG22_PAR, fig22_cpts()) fig, axes = plt.subplots(1, 2, figsize=(8.4, 3.6)) lp = {"a": "upper left", "f": "upper left", "e": "left", "b": "left", "c": "left"} for k, obs in enumerate((("c",), ("f",))): ax = axes[k] d = draw_dag(ax, FIG22_POS, FIG22_PAR, observed=obs, lpos=lp) sep = d_separated(FIG22_PAR, ["a"], ["b"], obs) mi = cond_mutual_info(P, order, ["a"], ["b"], obs) (path, blk), = trails(FIG22_PAR, "a", "b", obs) if blk: msg = "path a-e-f-b blocked\n" + "\n".join("at %s (%s)" % (v, why.split(",")[0]) for v, why in blk) else: msg = "path a-e-f-b open:\ne head-to-head, its\ndescendant c observed" stmt = (r"$a \perp b \mid %s$" % obs[0]) if sep else (r"$a, b$ dependent given $%s$" % obs[0]) d.text(1.5, -0.85, stmt, fontsize=12, color=OK if sep else NODE) d.text(1.5, -1.25, "I(a; b | %s) = %s" % (obs[0], _fmt_mi(mi)), fontsize=9, color="0.25") d.text(2.6, 0.05, msg, fontsize=8.5, color="0.3", ha="center") d.done(pad=0.1) ax.set_title("(%s) %s observed" % ("ab"[k], obs[0]), fontsize=10) fig.tight_layout() common.savefig(fig, 8, "fig8-22", out_dir) # Figure 8.26, extended by a grandparent g and a grandchild d that are outside the blanket. MB_POS = {"p1": (-0.75, 2), "p2": (0.75, 2), "xi": (0, 1), "q1": (-1.5, 1), "q2": (1.5, 1), "c1": (-0.75, 0), "c2": (0.75, 0), "g": (-1.5, 3), "d": (1.5, -1)} MB_PAR = {"g": (), "p1": ("g",), "p2": (), "q1": (), "q2": (), "xi": ("p1", "p2"), "c1": ("q1", "xi"), "c2": ("xi", "q2"), "d": ("c2",)} def fig8_26(out_dir=None): """Figure 8.26: the Markov blanket of x_i in a directed graph (parents, children, co-parents).""" mb = markov_blanket(MB_PAR, "xi") fig, ax = plt.subplots(figsize=(4.0, 4.3)) d = Diagram(ax) for v, xy in MB_POS.items(): inside = v in mb or v == "xi" d.node(v, xy, r"$\mathbf{x}_i$" if v == "xi" else None, fill=OBS if v in mb else "white", color=NODE if inside else GREY, ls="-" if inside else "--") for v, ps in MB_PAR.items(): for p in ps: inside = (v in mb or v == "xi") and (p in mb or p == "xi") d.edge(p, v, color=NODE if inside else GREY) d.text(-1.5, 3.42, "outside", fontsize=8, color="0.45") d.text(1.5, -1.42, "outside", fontsize=8, color="0.45") d.done() common.savefig(fig, 8, "fig8-26", out_dir) def fig8_28(out_dir=None): """Figure 8.28: the Markov blanket in an undirected graph is the set of neighbours.""" fig, ax = plt.subplots(figsize=(3.6, 3.6)) d = Diagram(ax) pos = {(i, j): (j, -i) for i in range(3) for j in range(3)} for (i, j), xy in pos.items(): centre = (i, j) == (1, 1) nbr = abs(i - 1) + abs(j - 1) == 1 d.node((i, j), xy, r"$x_i$" if centre else None, fill=OBS if nbr else "white", color=NODE if (centre or nbr) else GREY, ls="-" if (centre or nbr) else "--", lpos="upper left") for (i, j) in pos: for (di, dj) in ((0, 1), (1, 0)): q = (i + di, j + dj) if q in pos: inside = (1, 1) in ((i, j), q) d.edge((i, j), q, directed=False, color=NODE if inside else GREY) d.done() common.savefig(fig, 8, "fig8-28", out_dir) def fig8_30(out_dir=None): """Figure 8.30: image de-noising with the Ising model (8.42): ICM and the exact MAP by graph cut. beta = 1.0, eta = 2.1, h = 0 as in the book; each pixel flipped with probability 0.1 (numpy seed 8). Pixels that differ from the original: noisy 9.88%, ICM 4.52% after 2 sweeps, graph cut 2.05% (book: 4% and 1%). """ x, y = fig30_setup() xi, tr = icm(y) xg = graph_cut_map(y) err = lambda z: 100.0 * np.mean(z != x) fig, axes = plt.subplots(2, 2, figsize=(7.4, 5.9)) cmap = plt.matplotlib.colors.ListedColormap(["#f5d27a", "#1f1fe0"]) panels = [(x, "original"), (y, "noisy: %.1f%% of pixels flipped" % err(y)), (xi, "ICM (%d sweeps): %.2f%% wrong, E = %.0f" % (len(tr) - 1, err(xi), tr[-1])), (xg, "graph cut (exact MAP): %.2f%% wrong, E = %.0f" % (err(xg), ising_energy(xg, y)))] for ax, (img, t) in zip(axes.ravel(), panels): ax.imshow(img, cmap=cmap, vmin=-1, vmax=1, interpolation="nearest", rasterized=True) ax.set_title(t, fontsize=9.5) ax.set_xticks([]) ax.set_yticks([]) fig.tight_layout() common.savefig(fig, 8, "fig8-30", out_dir) return {"noisy": err(y), "icm": err(xi), "cut": err(xg), "E_icm": tr[-1], "E_cut": ising_energy(xg, y), "E_noisy": tr[0], "sweeps": len(tr) - 1} def fig8_32(out_dir=None): """Figures 8.32 and 8.33: a directed chain and its undirected version; a node with three parents and its moral graph, whose links are computed by moralize().""" fig, axes = plt.subplots(1, 4, figsize=(11.5, 2.4), gridspec_kw={"width_ratios": [1.7, 1.7, 1, 1]}) chain = {"x1": (), "x2": ("x1",), "x3": ("x2",), "x4": ("x3",)} lab = {"x1": "$x_1$", "x2": "$x_2$", "x3": "$x_{N-1}$", "x4": "$x_N$"} pos = {"x1": (0, 0), "x2": (1.2, 0), "x3": (3.0, 0), "x4": (4.2, 0)} for k, directed in enumerate((True, False)): ax = axes[k] d = Diagram(ax) for v, xy in pos.items(): d.node(v, xy, lab[v], lpos="above") d.node("d1", (1.95, 0), kind="point") d.node("d2", (2.25, 0), kind="point") d.edge("x1", "x2", directed) d.edge("x2", "d1", directed) d.edge("d2", "x3", directed) d.edge("x3", "x4", directed) d.text(2.1, 0.02, r"$\cdots$", fontsize=13, color=NODE) d.done() ax.set_title("(%s) %s chain" % ("ab"[k], "directed" if directed else "undirected"), fontsize=10) par = {"x1": (), "x2": (), "x3": (), "x4": ("x1", "x2", "x3")} pos2 = {"x1": (-1, 1.2), "x2": (0, 1.2), "x3": (1, 1.2), "x4": (0, 0)} lp = {"x1": "above", "x2": "above", "x3": "above", "x4": "below"} d = draw_dag(axes[2], pos2, par, labels={v: "$x_%s$" % v[1] for v in pos2}, lpos=lp) d.done() axes[2].set_title("(c) directed", fontsize=10) d = Diagram(axes[3]) for v, xy in pos2.items(): d.node(v, xy, "$x_%s$" % v[1], lpos=lp[v]) for e in sorted(moralize(par), key=lambda e: sorted(e)): u, v = sorted(e) married = not (u in par[v] or v in par[u]) if married and {u, v} == {"x1", "x3"}: # draw the long link as an arc above x2 axes[3].annotate("", xy=(1 - 0.2, 1.4), xytext=(-1 + 0.2, 1.4), arrowprops=dict(arrowstyle="-", color=common.BLUE, lw=2, connectionstyle="arc3,rad=-0.45")) d._grow(-1, 1.2, 1, 2.0) else: d.edge(u, v, directed=False, color=common.BLUE if married else NODE) d.done() axes[3].set_title("(d) moral graph (new links blue)", fontsize=10) fig.tight_layout() common.savefig(fig, 8, "fig8-32", out_dir) def fig8_40(out_dir=None): """Figure 8.40: the factor graph of p(x) = f_a(x1, x2) f_b(x1, x2) f_c(x2, x3) f_d(x3) (8.60).""" fig, ax = plt.subplots(figsize=(4.6, 2.6)) d = Diagram(ax) for k, xy in enumerate([(-1, 1), (0, 1), (1, 1)]): d.node("x%d" % (k + 1), xy, "$x_%d$" % (k + 1), lpos="above") for k, xy in enumerate([(-1.5, 0), (-0.5, 0), (0.5, 0), (1.5, 0)]): d.node("f" + "abcd"[k], xy, "$f_%s$" % "abcd"[k], kind="factor", lpos="below") for f, vs in {"fa": ("x1", "x2"), "fb": ("x1", "x2"), "fc": ("x2", "x3"), "fd": ("x3",)}.items(): for v in vs: d.edge(v, f, directed=False) d.done() common.savefig(fig, 8, "fig8-40", out_dir) FIG51_POS = {"x1": (0, 1.6), "x2": (2, 1.6), "x3": (4, 1.6), "x4": (2, 0), "fa": (1, 1.6), "fb": (3, 1.6), "fc": (2, 0.8)} def _fig51_diagram(ax): d = Diagram(ax) lp = {"x1": "above", "x2": "above", "x3": "above", "x4": "below", "fa": "below", "fb": "below", "fc": "right"} for v in ("x1", "x2", "x3", "x4"): d.node(v, FIG51_POS[v], "$x_%s$" % v[1], lpos=lp[v]) for f in ("fa", "fb", "fc"): d.node(f, FIG51_POS[f], "$f_%s$" % f[1], kind="factor", lpos=lp[f]) for f, vs in {"fa": ("x1", "x2"), "fb": ("x2", "x3"), "fc": ("x2", "x4")}.items(): for v in vs: d.edge(v, f, directed=False) return d def fig8_51(out_dir=None): """Figure 8.51: the factor graph p~(x) = f_a(x1, x2) f_b(x2, x3) f_c(x2, x4) used for sum-product.""" fig, ax = plt.subplots(figsize=(4.4, 2.8)) d = _fig51_diagram(ax) d.done() common.savefig(fig, 8, "fig8-51", out_dir) def _key(node): return ("x%d" % (node[1] + 1)) if node[0] == "v" else "f" + "abc"[node[1]] def fig8_52(out_dir=None): """Figure 8.52: sum-product on the graph of Figure 8.51 with root x3: (a) messages from the leaves to the root, (b) from the root back to the leaves, with their values for the integer tables FIG51_TABLES (binary variables; vectors list the states 0 and 1).""" fg = fig51_graph() res = sum_product(fg, root=2) sent = res["messages"] fig, axes = plt.subplots(1, 2, figsize=(10.0, 3.6)) for k, ax in enumerate(axes): d = _fig51_diagram(ax) part = sent[:6] if k == 0 else sent[6:] for j, (src, dst, m) in enumerate(part): a, b = FIG51_POS[_key(src)], FIG51_POS[_key(dst)] a, b = np.array(a), np.array(b) u = (b - a) / np.linalg.norm(b - a) nrm = np.array([-u[1], u[0]]) off = 0.32 if abs(u[1]) < 0.5 else 0.3 side = nrm * off * (1 if (u[0] > 0 or u[1] < 0) else -1) if abs(u[1]) > 0.5: side = np.array([-0.3, 0]) p0 = a + 0.3 * (b - a) + side p1 = a + 0.7 * (b - a) + side ax.annotate("", xy=p1, xytext=p0, arrowprops=dict(arrowstyle="-|>", color=common.BLUE, lw=1.8, mutation_scale=13)) mid = (p0 + p1) / 2 + (side * 0.75 if abs(u[1]) < 0.5 else np.array([-0.05, 0])) ax.text(mid[0], mid[1], "(%g, %g)" % tuple(m), fontsize=8, color=common.BLUE, ha="right" if abs(u[1]) > 0.5 else "center", va="center") d._grow(mid[0] - 0.6, mid[1] - 0.2, mid[0] + 0.3, mid[1] + 0.2) d.done(pad=0.05) ax.set_title("(%s) %s" % ("ab"[k], "leaves to the root $x_3$" if k == 0 else "root $x_3$ back to the leaves"), fontsize=10) Pm = res["marginals"] fig.text(0.5, 0.01, "Z = %g; p(x1) = (%.3f, %.3f), p(x2) = (%.3f, %.3f), p(x3) = (%.3f, %.3f), p(x4) = (%.3f, %.3f)" % ((res["Z"],) + tuple(np.concatenate(Pm))), ha="center", fontsize=9) fig.tight_layout(rect=(0, 0.05, 1, 1)) common.savefig(fig, 8, "fig8-52", out_dir) def fig8_53(out_dir=None): """Figure 8.53: the max-sum lattice (trellis) of the chain CHAIN_PSI (N = 5, K = 3): messages mu(x_n), back-pointers phi(x_n) (8.101) as black lines, and the back-tracked path (8.102).""" r = chain_max_sum(CHAIN_PSI) N, K = len(CHAIN_PSI) + 1, 3 cols = [common.RED, common.GREEN, common.BLUE] fig, ax = plt.subplots(figsize=(8.0, 3.9)) X = lambda n: 1.6 * n Y = lambda k: -1.1 * k for n in range(1, N): for k in range(K): j = r["phi"][n][k] on = r["path"][n] == k and r["path"][n - 1] == j ax.plot([X(n - 1) + 0.2, X(n) - 0.2], [Y(j), Y(k)], color="k" if on else "0.55", lw=3.2 if on else 1.2, zorder=1) for n in range(N): for k in range(K): on = r["path"][n] == k ax.add_patch(FancyBboxPatch((X(n) - 0.2, Y(k) - 0.2), 0.4, 0.4, boxstyle="round,pad=0,rounding_size=0.06", fc="white", ec=cols[k], lw=4 if on else 2.2, zorder=2)) ax.text(X(n), Y(k) - 0.36, "%.2f" % r["mu"][n][k], ha="center", va="top", fontsize=8, color="0.2") for k in range(K): ax.text(X(0) - 0.45, Y(k), "$k = %d$" % (k + 1), color=cols[k], ha="right", va="center", fontsize=12) for n in range(N): ax.text(X(n), Y(K - 1) - 0.8, "$x_%d$" % (n + 1), ha="center", va="top", fontsize=12) ax.annotate("", xy=(X(N - 1), 0.75), xytext=(X(0), 0.75), arrowprops=dict(arrowstyle="-|>", color="k", lw=2)) ax.text(X(N - 1) / 2, 0.85, r"max-sum messages $\mu(x_n)$ and back-pointers $\phi(x_n)$ (8.101), root $x_5$", ha="center", va="bottom", fontsize=9) ax.annotate("", xy=(X(0), -3.55), xytext=(X(N - 1), -3.55), arrowprops=dict(arrowstyle="-|>", color="0.3", lw=1.5, ls="--")) ax.text(X(N - 1) / 2, -3.65, "back-tracking (8.102): path x = (%s), ln max = %.3f" % (", ".join(str(p + 1) for p in r["path"]), r["logmax"]), ha="center", va="top", fontsize=9) ax.set_xlim(X(0) - 1.3, X(N - 1) + 0.5) ax.set_ylim(-4.0, 1.2) ax.set_aspect("equal") ax.axis("off") common.savefig(fig, 8, "fig8-53", out_dir) FIGURES = { "8.2": fig8_2, "8.6": fig8_6, "8.15": fig8_15, "8.21": fig8_21, "8.22": fig8_22, "8.26": fig8_26, "8.28": fig8_28, "8.30": fig8_30, "8.32": fig8_32, "8.40": fig8_40, "8.51": fig8_51, "8.52": fig8_52, "8.53": fig8_53, } if __name__ == "__main__": common.run(FIGURES, 8)