"""Bishop (2006), Chapter 3: Linear Models for Regression. Reproduces the numerical figures of Chapter 3 of C. M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006), pp. 137-173, written from scratch with numpy, scipy and matplotlib. Every algorithm is a small function meant to be read next to the lecture notes (ch03-linear-regression.html): * Section 3.1, basis functions: poly_basis, gauss_basis (3.4), sigmoid_basis (3.5), design_matrix. Section 3.1.1: ml_weights (the normal equations (3.15) solved by least squares, with noise precision (3.21)). Section 3.1.4: ridge_weights (3.28). * Section 3.1.3, sequential learning: lms_step (3.23). * Section 3.2, bias-variance: bias_variance (3.37)-(3.46) over L data sets. * Section 3.3, Bayesian linear regression: posterior (3.50)-(3.54), predictive (3.58)-(3.59), equivalent_kernel (3.60)-(3.62). * Sections 3.4 and 3.5, evidence: log_evidence (3.86), gamma_eff (3.91), evidence_update (3.92), (3.95), (3.98), evidence_iterate (the alpha/beta fixed point), evidence_path (alpha scanned). Data. Everything is synthetic, as in the book: sin(2 pi x) plus Gaussian noise, with x uniform on (0, 1) (Section 1.1; Figures 3.5, 3.6, 3.8, 3.9, 3.14, 3.16, 3.17), or the line a0 + a1 x with a0 = -0.3, a1 = 0.5, noise standard deviation 0.2 and x uniform on (-1, 1) (Figure 3.7, as in the text). Random numbers come from JSRandom (the seeded generator DK.rng of the browser demos), so the demos bayesian-regression-demo.html, bias-variance-demo.html and evidence-demo.html draw the same data number for number. Figures: 3.1 to 3.11, 3.14 to 3.17 (3.12 and 3.13 are schematic sketches, not reproduced). Run `python3 ch03_linear_regression.py --list`, `python3 ch03_linear_regression.py` (all figures, about half a minute) or `--fig 3.7 3.8`. Files go to ../img/ch03/. Settings chosen where the book gives none, and what differs: * Figure 3.1: polynomials x^1..x^10, Gaussians (3.4) with s = 0.2 centred at -1, -0.8, ..., 1, and sigmoids (3.5) with s = 0.1 and the same centres. * Figure 3.2: N = 3 observations, M = 2 basis functions, in 3-d (the book's sketch, with numbers). * Figure 3.3: the contours |w1|^q + |w2|^q = 1. Figure 3.4: the contours of a quadratic error centred at w_ML = (1.1, 0.55) with the constraint region |w1|^q + |w2|^q <= 0.8; the optimum w* is computed numerically (the lasso optimum has w1 = 0). * Figures 3.5 and 3.6: L = 100 data sets of N = 25 points, noise standard deviation 0.3, 24 Gaussian basis functions (s = 0.1, centres equally spaced on [0, 1]) plus a bias, ridge regression (3.28), ln lambda = 2.6, -0.31, -2.4 for Figure 3.5. The bias and variance integrals are taken over a grid of 200 equally spaced x on [0, 1], the density p(x) being uniform. The book's curves are lower than ours for the test error: its test error is about 0.065 whereas bias^2 + variance + sigma^2 here is about 0.1 for sigma = 0.3; the book's value is consistent with a noise variance of about 0.05, so we cannot recover its exact test noise. Our test error is bias^2 + variance + sigma^2, computed on 1000 fresh points per data set (it agrees with that sum to Monte Carlo accuracy). * Figures 3.8 and 3.9: 9 Gaussian basis functions (s = 0.1, centres equally spaced on [0, 1]) plus a bias, alpha = 2, beta = 25 (noise standard deviation 0.2), N = 1, 2, 4, 25. * Figures 3.10 and 3.11: alpha = 0.005, beta = 11.1; 200 equally spaced x on (-1, 1); the basis functions of Figure 3.1. * Figure 3.14: N = 10 points x = 0, 1/9, ..., 1, noise standard deviation 0.3 (beta = 1/0.09), alpha = 5e-3 as in the text; polynomial orders M = 0..8. * Figures 3.16 and 3.17: N = 25 points, noise standard deviation 0.3, 9 Gaussian basis functions (s = 0.1) plus a bias (10 parameters), beta fixed at 1/0.09. Figure 3.16(c) is an addition: the iteration (3.92)/(3.95)/(3.98) of alpha and beta from alpha = beta = 1. * Figure 3.15 is computed for a two-parameter example (a Hessian with eigenvalues 1 and 9). """ import math import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse from scipy import optimize import common from ch11_sampling import JSRandom NOISE_SD = 0.3 SEED = 3 # --------------------------------------------------------------------------- data def sine_data(rng, N, sd=NOISE_SD, lo=0.0, hi=1.0, grid=False): """N pairs (x, t), t = sin(2 pi x) + N(0, sd^2); x uniform on (lo, hi) or on an equally spaced grid.""" x = np.empty(N) t = np.empty(N) for n in range(N): x[n] = (lo + (hi - lo) * n / (N - 1)) if grid else lo + (hi - lo) * rng.u() t[n] = math.sin(2 * math.pi * x[n]) + sd * rng.norm() return x, t def line_data(rng, N, a0=-0.3, a1=0.5, sd=0.2): """Data for Figure 3.7: x ~ U(-1, 1), t = a0 + a1 x + N(0, sd^2).""" x = np.empty(N) t = np.empty(N) for n in range(N): x[n] = -1 + 2 * rng.u() t[n] = a0 + a1 * x[n] + sd * rng.norm() return x, t # --------------------------------------------------------------------------- Section 3.1: basis functions def poly_basis(x, powers): """phi_j(x) = x^j for the listed powers (Fig 3.1, left).""" return np.asarray(x, float)[:, None] ** np.asarray(powers)[None, :] def gauss_basis(x, mu, s): """Gaussian basis functions (3.4): exp(-(x - mu_j)^2 / (2 s^2)).""" return np.exp(-(np.asarray(x, float)[:, None] - np.asarray(mu)[None, :]) ** 2 / (2 * s * s)) def sigmoid_basis(x, mu, s): """Sigmoidal basis functions (3.5): sigma((x - mu_j) / s).""" return common.sigmoid((np.asarray(x, float)[:, None] - np.asarray(mu)[None, :]) / s) def design_matrix(x, kind, **kw): """Design matrix Phi (3.16), first column phi_0 = 1 (the bias).""" if kind == "poly": return poly_basis(x, np.arange(kw["M"] + 1)) if kind == "gauss": B = gauss_basis(x, kw["mu"], kw["s"]) elif kind == "sigmoid": B = sigmoid_basis(x, kw["mu"], kw["s"]) else: raise ValueError(kind) return np.hstack([np.ones((len(x), 1)), B]) def gauss_design(x, M, s=0.1): """Bias plus M Gaussian basis functions with centres equally spaced on [0, 1] (M + 1 columns).""" return design_matrix(x, "gauss", mu=np.linspace(0, 1, M), s=s) # --------------------------------------------------------------------------- Section 3.1: least squares def ml_weights(Phi, t): """Maximum likelihood weights: the solution of the normal equations (3.15), via the pseudo-inverse (Phi^T Phi) w = Phi^T t. Also returns the noise precision 1/beta_ML = mean squared residual (3.21).""" w = np.linalg.lstsq(Phi, t, rcond=None)[0] r = t - Phi @ w return w, 1.0 / np.mean(r ** 2) def ridge_weights(Phi, t, lam): """Regularised least squares (3.28): w = (lambda I + Phi^T Phi)^(-1) Phi^T t.""" M = Phi.shape[1] return np.linalg.solve(lam * np.eye(M) + Phi.T @ Phi, Phi.T @ t) def lms_step(w, phi, t, eta): """One step of sequential learning (3.23): w <- w + eta (t_n - w^T phi_n) phi_n.""" return w + eta * (t - w @ phi) * phi # --------------------------------------------------------------------------- Section 3.2: bias-variance def bias_variance(fits, h): """Decomposition (3.46) from L fitted curves. fits: (L, X) array of y^(l)(x_i), h: (X,) the regression function h(x_i), with the average over x taken as the mean over the grid (uniform p(x)). Returns (bias^2, variance) with bias^2 = mean_i (ybar_i - h_i)^2 (3.45) and variance = mean_i mean_l (y_i^(l) - ybar_i)^2 (3.46).""" ybar = fits.mean(0) return float(np.mean((ybar - h) ** 2)), float(np.mean((fits - ybar) ** 2)) def bias_variance_experiment(lnlams, L=100, N=25, M=24, sd=NOISE_SD, s=0.1, Ntest=1000, seed=1000, xgrid=None): """The experiment of Figures 3.5 and 3.6 for every ln lambda in lnlams. Returns dict(xg, h, fits[ln lambda] -> (L, len(xg)), bias2, var, test, curves) with the test error measured on Ntest fresh noisy points per data set.""" xg = np.linspace(0, 1, 200) if xgrid is None else xgrid h = np.sin(2 * np.pi * xg) Pg = gauss_design(xg, M, s) datasets = [] for l in range(L): rng = JSRandom(seed + l) datasets.append(sine_data(rng, N, sd)) rt = JSRandom(seed + 10 ** 6) xt, tt = sine_data(rt, Ntest, sd) Pt = gauss_design(xt, M, s) out = dict(xg=xg, h=h, lnlam=np.asarray(lnlams), fits={}, bias2=[], var=[], test=[]) for ll in lnlams: lam = math.exp(ll) F = np.array([Pg @ ridge_weights(gauss_design(x, M, s), t, lam) for x, t in datasets]) b2, v = bias_variance(F, h) te = np.mean([np.mean((Pt @ ridge_weights(gauss_design(x, M, s), t, lam) - tt) ** 2) for x, t in datasets]) out["fits"][ll] = F out["bias2"].append(b2) out["var"].append(v) out["test"].append(float(te)) for k in ("bias2", "var", "test"): out[k] = np.array(out[k]) return out # --------------------------------------------------------------------------- Section 3.3: Bayesian regression def posterior(Phi, t, alpha, beta): """Posterior over w for the prior N(0, alpha^-1 I): S_N^-1 = alpha I + beta Phi^T Phi (3.54), m_N = beta S_N Phi^T t (3.53). Returns (m_N, S_N).""" M = Phi.shape[1] SN = np.linalg.inv(alpha * np.eye(M) + beta * Phi.T @ Phi) return beta * SN @ Phi.T @ t, SN def predictive(Phi_new, mN, SN, beta): """Predictive distribution (3.58)-(3.59): mean m_N^T phi(x), variance 1/beta + phi(x)^T S_N phi(x).""" return Phi_new @ mN, 1.0 / beta + np.einsum("ij,jk,ik->i", Phi_new, SN, Phi_new) def equivalent_kernel(Phi_x, Phi_y, SN, beta): """Equivalent kernel (3.61): k(x, x') = beta phi(x)^T S_N phi(x'); the posterior mean (3.60) is y(x, m_N) = sum_n k(x, x_n) t_n.""" return beta * Phi_x @ SN @ Phi_y.T # --------------------------------------------------------------------------- Sections 3.4, 3.5: evidence def log_evidence(Phi, t, alpha, beta): """ln p(t | alpha, beta) (3.86) = M/2 ln alpha + N/2 ln beta - E(m_N) - 1/2 ln |A| - N/2 ln 2 pi, with A = alpha I + beta Phi^T Phi (3.81) and E(m_N) = beta/2 ||t - Phi m_N||^2 + alpha/2 m_N^T m_N (3.82).""" N, M = Phi.shape A = alpha * np.eye(M) + beta * Phi.T @ Phi mN = beta * np.linalg.solve(A, Phi.T @ t) E = beta / 2 * np.sum((t - Phi @ mN) ** 2) + alpha / 2 * mN @ mN return (M / 2 * math.log(alpha) + N / 2 * math.log(beta) - E - 0.5 * np.linalg.slogdet(A)[1] - N / 2 * math.log(2 * math.pi)) def gamma_eff(Phi, alpha, beta): """Effective number of parameters (3.91): gamma = sum_i lambda_i / (alpha + lambda_i), where the lambda_i are the eigenvalues of beta Phi^T Phi (3.87).""" lam = np.linalg.eigvalsh(beta * Phi.T @ Phi) return float(np.sum(lam / (alpha + lam))) def evidence_update(Phi, t, alpha, beta): """One pass of the re-estimation equations: m_N from (3.53), gamma (3.91), then alpha <- gamma / (m_N^T m_N) (3.92) and 1/beta <- sum_n (t_n - m_N^T phi_n)^2 / (N - gamma) (3.95).""" N = len(t) mN, _ = posterior(Phi, t, alpha, beta) g = gamma_eff(Phi, alpha, beta) return g / (mN @ mN), (N - g) / np.sum((t - Phi @ mN) ** 2), g def evidence_iterate(Phi, t, alpha=1.0, beta=1.0, iters=200, tol=1e-10, fix_beta=False): """Iterate (3.92) and (3.95) (the evidence approximation, Sections 3.5.2 and 3.5.3). Returns the final (alpha, beta, gamma) and the trace [(alpha, beta, gamma, ln evidence)].""" trace = [] for _ in range(iters): a2, b2, g = evidence_update(Phi, t, alpha, beta) trace.append((alpha, beta, g, log_evidence(Phi, t, alpha, beta))) if fix_beta: b2 = beta done = abs(a2 - alpha) < tol * alpha and abs(b2 - beta) < tol * beta alpha, beta = a2, b2 if done: break return alpha, beta, gamma_eff(Phi, alpha, beta), trace def evidence_path(Phi, t, beta, lnalphas): """For each ln alpha: m_N, gamma, ln evidence, 2 alpha E_W(m_N) = alpha m_N^T m_N (Figures 3.16, 3.17).""" rows = [] for la in lnalphas: a = math.exp(la) mN, _ = posterior(Phi, t, a, beta) rows.append((mN, gamma_eff(Phi, a, beta), log_evidence(Phi, t, a, beta), a * mN @ mN)) return rows # --------------------------------------------------------------------------- figures COLS = ["#0000ff", "#ff0000", "#00c000", "#b0a800", "#c000c0", "#00a8b8"] def _box(ax, ticks=None): for sp in ax.spines.values(): sp.set_linewidth(1) def fig3_1(out_dir=None): """Figure 3.1: polynomial, Gaussian and sigmoidal basis functions.""" x = np.linspace(-1, 1, 400) mu = np.linspace(-1, 1, 11) fig, axes = common.panel_grid(1, 3, size=(2.9, 2.7)) for j in range(1, 11): axes[0, 0].plot(x, x ** j, color=COLS[j % 6], lw=1.3) G, S = gauss_basis(x, mu, 0.2), sigmoid_basis(x, mu, 0.1) for j in range(11): axes[0, 1].plot(x, G[:, j], color=COLS[j % 6], lw=1.3) axes[0, 2].plot(x, S[:, j], color=COLS[j % 6], lw=1.3) for a, lab, yl in zip(axes[0], "abc", [(-1, 1), (0, 1), (0, 1)]): a.set_xlim(-1, 1) a.set_ylim(*yl) a.set_xlabel("$x$") common.panel_label(a, lab, "upper left" if lab != "a" else "lower right") axes[0, 0].set_title("polynomial $x^j$", fontsize=9.5) axes[0, 1].set_title("Gaussian (3.4), $s = 0.2$", fontsize=9.5) axes[0, 2].set_title("sigmoidal (3.5), $s = 0.1$", fontsize=9.5) fig.tight_layout() common.savefig(fig, 3, "fig3-1", out_dir) def fig3_2(out_dir=None): """Figure 3.2: least squares as orthogonal projection of t onto the subspace spanned by the columns of Phi.""" Phi = np.array([[1.0, 0.2], [1.0, 1.0], [1.0, 1.8]]) * np.array([1, 1.2]) Phi = np.array([[1.0, 0.1], [0.4, 1.0], [0.2, 0.3]]) t = np.array([0.6, 0.5, 1.1]) w, _ = ml_weights(Phi, t) y = Phi @ w resid = t - y fig = plt.figure(figsize=(4.4, 4.0)) ax = fig.add_subplot(111, projection="3d") u = np.linspace(-0.3, 1.3, 2) U, V = np.meshgrid(u, u) P = U[..., None] * Phi[:, 0] + V[..., None] * Phi[:, 1] ax.plot_surface(P[..., 0], P[..., 1], P[..., 2], color=common.BLUE, alpha=0.15, rasterized=True) for v, lab, c in [(Phi[:, 0], r"$\varphi_1$", common.GREEN), (Phi[:, 1], r"$\varphi_2$", common.GREEN), (t, r"$\mathbf{t}$", common.RED), (y, r"$\mathbf{y}$", common.BLUE)]: ax.quiver(0, 0, 0, *v, color=c, arrow_length_ratio=0.08, lw=1.6) ax.text(*(v * 1.07), lab, fontsize=12) ax.plot(*zip(t, y), "k--", lw=1) ax.set_title("$\\mathbf{y} = \\Phi\\mathbf{w}_{ML}$ is the orthogonal projection of $\\mathbf{t}$\n" "(residual $\\cdot$ $\\varphi_j$ = %.1e, %.1e)" % tuple(Phi.T @ resid), fontsize=8.5) ax.set_axis_off() ax.set_box_aspect((1, 1, 1)) common.savefig(fig, 3, "fig3-2", out_dir) return dict(orth=float(np.abs(Phi.T @ resid).max())) def fig3_3(out_dir=None): """Figure 3.3: contours of the regulariser sum_j |w_j|^q for q = 0.5, 1, 2, 4.""" fig, axes = common.panel_grid(1, 4, size=(2.3, 2.3)) th = np.linspace(0, 2 * np.pi, 2000) for a, q, lab in zip(axes[0], [0.5, 1, 2, 4], "abcd"): c, s = np.cos(th), np.sin(th) r = (np.abs(c) ** q + np.abs(s) ** q) ** (-1.0 / q) a.plot(r * c, r * s, color=common.BLUE, lw=1.6) a.set_xlim(-1.5, 1.5) a.set_ylim(-1.5, 1.5) a.set_aspect("equal") a.set_xticks([]) a.set_yticks([]) a.set_title("$q = %g$" % q, fontsize=10) common.panel_label(a, lab, "upper left") fig.tight_layout() common.savefig(fig, 3, "fig3-3", out_dir) def fig3_4(out_dir=None): """Figure 3.4: error contours and the constraint region, quadratic (q = 2) versus lasso (q = 1).""" wml = np.array([1.1, 0.55]) H = np.array([[1.0, 0.7], [0.7, 2.2]]) fig, axes = common.panel_grid(1, 2, size=(3.2, 3.2)) g = np.linspace(-1.6, 1.8, 300) X, Y = np.meshgrid(g, g) E = 0.5 * (H[0, 0] * (X - wml[0]) ** 2 + 2 * H[0, 1] * (X - wml[0]) * (Y - wml[1]) + H[1, 1] * (Y - wml[1]) ** 2) out = {} for a, q, lab in zip(axes[0], [2, 1], "ab"): eta = 0.8 cons = [dict(type="ineq", fun=lambda w, q=q: eta - np.sum(np.abs(w) ** q))] res = optimize.minimize(lambda w: 0.5 * (w - wml) @ H @ (w - wml), x0=[0.2, 0.2], constraints=cons, method="SLSQP", options=dict(ftol=1e-14, maxiter=500)) ws = res.x a.contour(X, Y, E, levels=[0.1, 0.3, 0.55, 0.9], colors=common.BLUE, linewidths=1.2) th = np.linspace(0, 2 * np.pi, 2000) c, s = np.cos(th), np.sin(th) r = (eta / (np.abs(c) ** q + np.abs(s) ** q)) ** (1.0 / q) a.fill(r * c, r * s, color=common.CYAN, alpha=0.3) a.plot(r * c, r * s, color=common.CYAN, lw=1.2) a.plot(*wml, "x", color="k") a.plot(*ws, "o", color=common.RED, ms=5) a.annotate("$\\mathbf{w}^\\star$", ws, textcoords="offset points", xytext=(-14, -14), fontsize=11) a.axhline(0, color="k", lw=0.5) a.axvline(0, color="k", lw=0.5) a.set_xlabel("$w_1$") a.set_ylabel("$w_2$") a.set_aspect("equal") a.set_xlim(-1.6, 1.8) a.set_ylim(-1.6, 1.8) a.set_title("$q = %d$: $\\mathbf{w}^\\star = (%.3f, %.3f)$" % (q, ws[0], ws[1]), fontsize=9) common.panel_label(a, lab, "lower right") out["q%d" % q] = ws.tolist() fig.tight_layout() common.savefig(fig, 3, "fig3-4", out_dir) return out def fig3_5_6(out_dir=None, L=100): """Figures 3.5 and 3.6 together (bias-variance); returns the numbers.""" lnl_show = [2.6, -0.31, -2.4] lnl = np.round(np.linspace(-3, 3, 61), 2) allv = sorted(set(list(lnl) + lnl_show)) R = bias_variance_experiment(allv, L=L) return R, lnl_show def fig3_5(out_dir=None): """Figure 3.5: 20 of the 100 fits (left) and their average with sin(2 pi x) (right), ln lambda = 2.6, -0.31, -2.4.""" R, show = fig3_5_6() fig, axes = common.panel_grid(3, 2, size=(3.0, 2.2), sharex=True) xg = R["xg"] for i, ll in enumerate(show): F = R["fits"][ll] for l in range(20): axes[i, 0].plot(xg, F[l], color=common.RED, lw=0.8, alpha=0.75) axes[i, 1].plot(xg, np.sin(2 * np.pi * xg), color=common.GREEN, lw=1.6) axes[i, 1].plot(xg, F.mean(0), color=common.RED, lw=1.6) for j in (0, 1): axes[i, j].set_ylim(-1.5, 1.5) axes[i, j].set_xlim(0, 1) axes[i, j].set_yticks([-1, 0, 1]) axes[i, 0].set_ylabel("$t$") k = list(R["lnlam"]).index(ll) axes[i, 0].set_title("$\\ln\\lambda = %g$ (bias$^2$ %.4f, var %.4f)" % (ll, R["bias2"][k], R["var"][k]), fontsize=8.5, loc="left") common.panel_label(axes[i, 0], "abc"[i], "lower left") common.panel_label(axes[i, 1], "def"[i], "lower left") axes[2, 0].set_xlabel("$x$") axes[2, 1].set_xlabel("$x$") fig.tight_layout() common.savefig(fig, 3, "fig3-5", out_dir) return {str(ll): (float(R["bias2"][list(R["lnlam"]).index(ll)]), float(R["var"][list(R["lnlam"]).index(ll)])) for ll in show} def fig3_6(out_dir=None): """Figure 3.6: squared bias, variance, their sum and the test error against ln lambda.""" R, _ = fig3_5_6() fig, ax = plt.subplots(figsize=(4.6, 3.3)) x = R["lnlam"] ax.plot(x, R["bias2"], color=common.BLUE, label="(bias)$^2$") ax.plot(x, R["var"], color=common.RED, label="variance") ax.plot(x, R["bias2"] + R["var"], color=common.MAGENTA, label="(bias)$^2$ + variance") ax.plot(x, R["test"], color="k", label="test error") ax.axhline(NOISE_SD ** 2, color="k", lw=0.8, ls=":") ax.text(x[0] + 0.05, NOISE_SD ** 2 + 0.004, "noise variance $\\sigma^2$ = %.2f" % NOISE_SD ** 2, fontsize=7.5) ax.set_xlabel(r"$\ln\lambda$") ax.set_ylim(0, 0.15) ax.set_xlim(x[0], x[-1]) ax.legend(loc="upper center", fontsize=8, frameon=False, ncol=2) fig.tight_layout() common.savefig(fig, 3, "fig3-6", out_dir) k = int(np.argmin(R["bias2"] + R["var"])) return dict(best_lnlam=float(x[k]), bias2=float(R["bias2"][k]), var=float(R["var"][k]), test=float(R["test"][k])) def _heat(ax, f, lim=1, n=140): g = np.linspace(-lim, lim, n) W0, W1 = np.meshgrid(g, g) Z = f(W0, W1) ax.imshow(Z, origin="lower", extent=(-lim, lim, -lim, lim), cmap="jet", aspect="auto", rasterized=True) def fig3_7(out_dir=None): """Figure 3.7: sequential Bayesian learning of a line (alpha = 2, beta = 25): likelihood, prior/posterior, data space.""" alpha, beta, N = 2.0, 25.0, 20 rng = JSRandom(SEED) x, t = line_data(rng, N) rows = [0, 1, 2, 20] fig, axes = common.panel_grid(4, 3, size=(2.5, 2.4)) xs = np.array([-1.0, 1.0]) Phi = np.column_stack([np.ones(N), x]) sample_rng = np.random.default_rng(7) out = {} for r, n in enumerate(rows): mN, SN = posterior(Phi[:n], t[:n], alpha, beta) if n else (np.zeros(2), np.eye(2) / alpha) if n: # likelihood of the latest point alone, as a function of w xl, tl = x[n - 1], t[n - 1] _heat(axes[r, 0], lambda a, b: np.exp(-beta / 2 * (tl - a - b * xl) ** 2)) axes[r, 0].plot(-0.3, 0.5, "+", color="w", ms=9, mew=2) else: axes[r, 0].axis("off") Sinv = np.linalg.inv(SN) _heat(axes[r, 1], lambda a, b: np.exp(-0.5 * (Sinv[0, 0] * (a - mN[0]) ** 2 + 2 * Sinv[0, 1] * (a - mN[0]) * (b - mN[1]) + Sinv[1, 1] * (b - mN[1]) ** 2))) if n: axes[r, 1].plot(-0.3, 0.5, "+", color="w", ms=9, mew=2) W = sample_rng.multivariate_normal(mN, SN, size=6) for w in W: axes[r, 2].plot(xs, w[0] + w[1] * xs, color="red", lw=1.6) if n: axes[r, 2].plot(x[:n], t[:n], "o", mfc="none", mec="blue", ms=5, mew=1.4) axes[r, 2].set_xlim(-1, 1) axes[r, 2].set_ylim(-1, 1) for c in (0, 1): axes[r, c].set_xlim(-1, 1) axes[r, c].set_ylim(-1, 1) axes[r, c].set_xlabel("$w_0$", labelpad=-6) axes[r, c].set_ylabel("$w_1$", labelpad=-8) axes[r, c].set_xticks([-1, 0, 1]) axes[r, c].set_yticks([-1, 0, 1]) axes[r, 2].set_xlabel("$x$", labelpad=-6) axes[r, 2].set_ylabel("$y$", labelpad=-8) axes[r, 2].set_xticks([-1, 0, 1]) axes[r, 2].set_yticks([-1, 0, 1]) out[n] = dict(mean=mN.tolist(), sd=np.sqrt(np.diag(SN)).tolist()) for c, h in enumerate(["likelihood", "prior/posterior", "data space"]): axes[0, c].set_title(h, fontsize=10) fig.tight_layout() common.savefig(fig, 3, "fig3-7", out_dir) return out def _sin_gauss_fit(N, seed=5, noise=0.2, alpha=2.0, beta=25.0, M=9, s=0.1): x, t = sine_data(JSRandom(seed), N, noise) Phi = gauss_design(x, M, s) mN, SN = posterior(Phi, t, alpha, beta) return x, t, mN, SN, beta, M, s def fig3_8(out_dir=None): """Figure 3.8: predictive distribution (mean and +-1 sd) for 9 Gaussian basis functions, N = 1, 2, 4, 25.""" fig, axes = common.panel_grid(2, 2, size=(3.1, 2.5)) xx = np.linspace(0, 1, 200) Px = None for a, N, lab in zip(axes.ravel(), [1, 2, 4, 25], "abcd"): x, t, mN, SN, beta, M, s = _sin_gauss_fit(N) Px = gauss_design(xx, M, s) mu, var = predictive(Px, mN, SN, beta) a.plot(xx, np.sin(2 * np.pi * xx), color=common.GREEN, lw=1.5) a.fill_between(xx, mu - np.sqrt(var), mu + np.sqrt(var), color=common.RED, alpha=0.25, lw=0) a.plot(xx, mu, color=common.RED, lw=1.6) a.plot(x, t, "o", mfc="none", mec=common.BLUE, ms=5, mew=1.3) a.set_xlim(0, 1) a.set_ylim(-1.5, 1.5) a.set_xlabel("$x$") a.set_ylabel("$t$") common.panel_label(a, lab, "upper right") a.set_title("$N = %d$" % N, fontsize=9.5, loc="left") fig.tight_layout() common.savefig(fig, 3, "fig3-8", out_dir) def fig3_9(out_dir=None): """Figure 3.9: five functions y(x, w) with w drawn from the posterior, N = 1, 2, 4, 25.""" fig, axes = common.panel_grid(2, 2, size=(3.1, 2.5)) xx = np.linspace(0, 1, 200) rng = np.random.default_rng(11) for a, N, lab in zip(axes.ravel(), [1, 2, 4, 25], "abcd"): x, t, mN, SN, beta, M, s = _sin_gauss_fit(N) Px = gauss_design(xx, M, s) for w in rng.multivariate_normal(mN, SN, size=5): a.plot(xx, Px @ w, color=common.RED, lw=1.3) a.plot(x, t, "o", mfc="none", mec=common.BLUE, ms=5, mew=1.3) a.set_xlim(0, 1) a.set_ylim(-1.5, 1.5) a.set_xlabel("$x$") a.set_ylabel("$t$") common.panel_label(a, lab, "upper right") a.set_title("$N = %d$" % N, fontsize=9.5, loc="left") fig.tight_layout() common.savefig(fig, 3, "fig3-9", out_dir) KERNEL_ALPHA, KERNEL_BETA = 5e-3, 11.1 def _kernel(kind): x = np.linspace(-1, 1, 200) mu = np.linspace(-1, 1, 11) if kind == "gauss": Phi = design_matrix(x, "gauss", mu=mu, s=0.2) elif kind == "sigmoid": Phi = design_matrix(x, "sigmoid", mu=mu, s=0.1) else: Phi = design_matrix(x, "poly", M=10) _, SN = posterior(Phi, np.zeros(len(x)), KERNEL_ALPHA, KERNEL_BETA) return x, equivalent_kernel(Phi, Phi, SN, KERNEL_BETA), Phi, SN def fig3_10(out_dir=None): """Figure 3.10: the equivalent kernel k(x, x') for the Gaussian basis functions, with three slices.""" x, K, _, _ = _kernel("gauss") fig = plt.figure(figsize=(7.2, 3.3)) gs = fig.add_gridspec(1, 2, width_ratios=[1.1, 1]) a = fig.add_subplot(gs[0]) b = fig.add_subplot(gs[1]) b.imshow(K, origin="lower", extent=(-1, 1, -1, 1), cmap="jet", rasterized=True) for xi in [0.6, -0.1, -0.6]: j = int(np.argmin(np.abs(x - xi))) a.plot(x, K[j] / K.max() * 0.6 + x[j], color="blue", lw=1.5) a.axhline(x[j], color=common.ORANGE, lw=1) b.axhline(x[j], color=common.ORANGE, lw=1) a.set_xlim(-1, 1) a.set_ylim(-1, 1) a.set_xlabel("$x'$") a.set_ylabel("$x$ (slices of $k(x, x')$ at $x$ = 0.6, -0.1, -0.6)", fontsize=8) b.set_xlabel("$x'$") b.set_ylabel("$x$") common.panel_label(a, "a", "upper left") common.panel_label(b, "b", "upper left") fig.tight_layout() common.savefig(fig, 3, "fig3-10", out_dir) return dict(k_max=float(K.max())) def fig3_11(out_dir=None): """Figure 3.11: equivalent kernels k(0, x') for polynomial (left) and sigmoidal (right) basis functions.""" fig, axes = common.panel_grid(1, 2, size=(3.4, 2.5)) for a, kind, lab in zip(axes[0], ["poly", "sigmoid"], "ab"): x, K, _, _ = _kernel(kind) j = int(np.argmin(np.abs(x))) a.plot(x, K[j], color=common.RED, lw=1.6) a.set_xlim(-1, 1) a.set_xlabel("$x'$") a.set_ylabel("$k(0, x')$") a.set_title({"poly": "polynomial", "sigmoid": "sigmoidal"}[kind], fontsize=9.5) common.panel_label(a, lab, "upper left") fig.tight_layout() common.savefig(fig, 3, "fig3-11", out_dir) def fig3_14(out_dir=None): """Figure 3.14: log evidence against the polynomial order M (alpha = 5e-3, beta = 1/0.3^2, N = 10).""" x, t = sine_data(JSRandom(SEED + 40), 10, 0.3, grid=True) ev = [log_evidence(design_matrix(x, "poly", M=M), t, 5e-3, 1 / 0.09) for M in range(9)] fig, axes = common.panel_grid(1, 2, size=(3.4, 2.7)) a, b = axes[0] a.plot(range(9), ev, "o-", color=common.RED, mfc="none", mew=1.6) a.set_xlabel("$M$") a.set_ylabel(r"$\ln p(\mathbf{t}\mid\alpha,\beta)$") a.set_xticks(range(0, 9, 2)) common.panel_label(a, "a", "lower right") xx = np.linspace(0, 1, 200) for M, c in [(0, common.BLUE), (3, common.RED), (8, common.GREEN)]: mN, _ = posterior(design_matrix(x, "poly", M=M), t, 5e-3, 1 / 0.09) b.plot(xx, design_matrix(xx, "poly", M=M) @ mN, color=c, label="$M = %d$" % M) b.plot(xx, np.sin(2 * np.pi * xx), "k:", lw=1) b.plot(x, t, "o", mfc="none", mec=common.BLUE, ms=5) b.set_ylim(-1.6, 1.6) b.set_xlabel("$x$") b.legend(fontsize=7.5, loc="upper right") common.panel_label(b, "b", "lower left") fig.tight_layout() common.savefig(fig, 3, "fig3-14", out_dir) return dict(best_M=int(np.argmax(ev)), ev=[round(float(v), 2) for v in ev]) def fig3_15(out_dir=None): """Figure 3.15: likelihood contours (red) and prior contours (green) in the eigenbasis of the Hessian.""" lam = np.array([1.0, 9.0]) alpha = 2.0 wml = np.array([1.5, 0.9]) # in the rotated coordinates (u1, u2) mN = lam * wml / (alpha + lam) fig, ax = plt.subplots(figsize=(3.8, 3.8)) g = np.linspace(-1.2, 2.6, 300) X, Y = np.meshgrid(g, g) E = 0.5 * (lam[0] * (X - wml[0]) ** 2 + lam[1] * (Y - wml[1]) ** 2) ax.contour(X, Y, E, levels=[0.3, 1.0, 2.5], colors=common.RED, linewidths=1.3) th = np.linspace(0, 2 * np.pi, 200) for r in (0.6, 1.2): ax.plot(r * np.cos(th), r * np.sin(th), color=common.GREEN, lw=1.3) ax.plot(*wml, "o", color=common.RED, ms=4) ax.plot(*mN, "o", color="k", ms=4) ax.annotate(r"$\mathbf{w}_{ML}$", wml, xytext=(6, 4), textcoords="offset points", fontsize=11) ax.annotate(r"$\mathbf{w}_{MAP}=\mathbf{m}_N$", mN, xytext=(6, -14), textcoords="offset points", fontsize=11) ax.set_xlabel("$u_1$ (small $\\lambda_1 = 1$: poorly determined)") ax.set_ylabel("$u_2$ ($\\lambda_2 = 9$: well determined)") ax.set_aspect("equal") ax.set_title("$m_{N,i} = \\lambda_i/(\\alpha+\\lambda_i)\\, w_{ML,i}$: shrink factors %.2f, %.2f ($\\alpha = %g$)" % (lam[0] / (alpha + lam[0]), lam[1] / (alpha + lam[1]), alpha), fontsize=8) fig.tight_layout() common.savefig(fig, 3, "fig3-15", out_dir) def _evidence_data(): x, t = sine_data(JSRandom(SEED + 70), 25, 0.3) return x, t, gauss_design(x, 9, 0.1) def fig3_16(out_dir=None): """Figure 3.16: gamma and 2 alpha E_W against ln alpha; log evidence and test error; (c) the alpha, beta iteration.""" x, t, Phi = _evidence_data() beta = 1 / 0.09 la = np.linspace(-5, 5, 201) rows = evidence_path(Phi, t, beta, la) gam = np.array([r[1] for r in rows]) twoaE = np.array([r[3] for r in rows]) ev = np.array([r[2] for r in rows]) xt, tt = sine_data(JSRandom(SEED + 71), 1000, 0.3) Pt = gauss_design(xt, 9, 0.1) test = np.array([np.mean((Pt @ r[0] - tt) ** 2) for r in rows]) fig, axes = common.panel_grid(1, 3, size=(3.2, 2.7)) a, b, c = axes[0] a.plot(la, gam, color=common.RED, label="$\\gamma$") a.plot(la, twoaE, color=common.BLUE, label="$2\\alpha E_W(\\mathbf{m}_N)$") a.set_ylim(0, 12) a.set_xlabel(r"$\ln\alpha$") a.legend(fontsize=8, loc="upper right") common.panel_label(a, "a", "upper left") b.plot(la, ev, color=common.RED, label=r"$\ln p(\mathbf{t}|\alpha,\beta)$") b2 = b.twinx() b2.plot(la, test, color=common.BLUE, label="test error") b2.set_ylim(0, test.max() * 1.05) b.set_xlabel(r"$\ln\alpha$") b.set_ylabel("log evidence", color=common.RED) b2.set_ylabel("test error", color=common.BLUE) b.set_ylim(ev.max() - 25, ev.max() + 2) common.panel_label(b, "b", "lower left") al, be, g, trace = evidence_iterate(Phi, t, 1.0, 1.0, fix_beta=False) it = np.arange(len(trace)) c.plot(it, [np.log(r[0]) for r in trace], "o-", color=common.RED, ms=3, label=r"$\ln\alpha$") c.plot(it, [np.log(r[1]) for r in trace], "s-", color=common.BLUE, ms=3, label=r"$\ln\beta$") c.plot(it, [r[2] for r in trace], "^-", color=common.GREEN, ms=3, label=r"$\gamma$") c.set_xlabel("iteration") c.legend(fontsize=8) common.panel_label(c, "c", "upper right") fig.tight_layout() common.savefig(fig, 3, "fig3-16", out_dir) cross = la[int(np.argmin(np.abs(gam - twoaE)))] return dict(ln_alpha_cross=float(cross), ln_alpha_evidence=float(la[int(np.argmax(ev))]), ln_alpha_test=float(la[int(np.argmin(test))]), alpha_beta_fixed_point=(al, be, g), iters=len(trace)) def fig3_17(out_dir=None): """Figure 3.17: the 10 posterior-mean weights against the effective number of parameters gamma (alpha from infinity to 0).""" x, t, Phi = _evidence_data() beta = 1 / 0.09 la = np.linspace(14, -12, 600) rows = evidence_path(Phi, t, beta, la) gam = np.array([r[1] for r in rows]) W = np.array([r[0] for r in rows]) fig, ax = plt.subplots(figsize=(4.8, 3.5)) order = np.argsort(-W[-1]) for i in range(10): ax.plot(gam, W[:, i], color=COLS[i % 6], lw=1.3) ax.text(10.15, W[-1, i], str(i), color=COLS[i % 6], fontsize=8, va="center") ax.set_xlabel(r"$\gamma$") ax.set_ylabel("$w_i$") ax.set_xlim(0, 10.6) fig.tight_layout() common.savefig(fig, 3, "fig3-17", out_dir) return dict(gamma_end=float(gam[-1])) FIGURES = { "3.1": fig3_1, "3.2": fig3_2, "3.3": fig3_3, "3.4": fig3_4, "3.5": fig3_5, "3.6": fig3_6, "3.7": fig3_7, "3.8": fig3_8, "3.9": fig3_9, "3.10": fig3_10, "3.11": fig3_11, "3.14": fig3_14, "3.15": fig3_15, "3.16": fig3_16, "3.17": fig3_17, } if __name__ == "__main__": common.run(FIGURES, 3)