"""Bishop (2006), Chapter 6: Kernel Methods. Reproduces the numerical figures of Chapter 6 of C. M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006), pp. 291-323, written from scratch with numpy, scipy and matplotlib. Every algorithm is a small function meant to be read next to the lecture notes (ch06-kernel-methods.html): * Section 6.1, dual representations: ridge_primal and ridge_dual (the dual solution a = (K + lambda I)^(-1) t of (6.8) and the prediction (6.9)). * Section 6.2, constructing kernels: kernels from basis functions (6.10) (basis_kernel), the polynomial kernel (6.12), the Gaussian kernel (6.23), the sigmoidal kernel (6.37), the closure rules (6.13)-(6.22) (see checks_ch06.py) and the exponential kernel (6.56). * Section 6.3, radial basis functions: normalised basis functions (6.41) and the Nadaraya-Watson model (6.45)-(6.48) with its Gaussian-mixture conditional density (nadaraya_watson). * Section 6.4, Gaussian processes: kernel classes with analytic derivatives (GaussKernel, ExpKernel; the covariance (6.63) and the ARD kernel (6.71)-(6.72)), gp_predict ((6.66)-(6.67)), gp_log_marginal ((6.69) and its gradient (6.70)), gp_fit (maximum marginal likelihood by L-BFGS), and GP classification with the Laplace approximation: laplace_mode (Newton-Raphson (6.83), fixed point (6.84)), laplace_predict ((6.87)-(6.88) and (4.153)) and laplace_evidence ((6.90) with the gradient (6.91)-(6.94)). Figures (the code builds every numerical figure of the chapter): 6.1 to 6.12. Choices made where the book does not say (the functions are chosen to look like the book's, not copied): * 6.1: polynomial basis x^j, j = 1..10; 21 Gaussian basis functions of width 0.1 spaced 0.1 apart (every second one is drawn in the top panel; the kernel uses all 21, which gives the peak height 1.77 of the book); 11 logistic sigmoids sigma((x - mu)/0.1) spaced 0.2 apart (peak 5.5 as in the errata-corrected figure). Kernel plotted against x with x' = -0.5, 0, 0. * 6.2: three Gaussian basis functions with centres -0.5, 0.1, 0.65 and widths 0.1, 0.13, 0.25 (read off the figure). * 6.3: the sinusoidal data set of Appendix A with N = 10 equally spaced inputs, noise sd 0.3 (seeded), Gaussian kernel sd 0.1. * 6.4, 6.5: prior samples on [-1, 1]; the same standard-normal draws are used in every panel (like the book, where the panels are visibly related). Exponential kernel with theta = 4. * 6.6: Gaussian kernel (theta0, theta1) = (1, 4), noise sd 0.5, 12 random inputs. * 6.7: correlation 0.5 between the two function values, t1 = 0.95. * 6.8: the first 7 of N = 10 equally spaced points of the sinusoidal data (the three right-most omitted, as in the book), kernel (6.63) with theta0 = 1, theta1 = 12, beta^(-1) = 0.3^2. * 6.9: ARD prior samples with the stated eta, drawn over [-4, 4]^2 (the book's wiggles imply a larger input range than the tick labels show). * 6.10: the book's recipe (x1 Gaussian, t = sin(2 pi x1) + noise, x2 = x1 + noise, x3 independent Gaussian) with our own scales; the marginal likelihood is maximised by L-BFGS-B (100 iterations) instead of scaled conjugate gradients, so the curves differ in detail. * 6.11: a sample of a(x) from a GP with Gaussian kernel (theta0 = 16, theta1 = 16) and sigma(a). * 6.12: the book's data set (Appendix A, Figure A.7) is not available, so a synthetic two-class set of the same type is made (blue: one Gaussian; red: mixture of two Gaussians; equal priors, 100 points each); the kernel amplitude and length scale are learnt by maximising the Laplace evidence (6.90) with the analytic gradient. Run `python3 ch06_kernel_methods.py --list` for the figure list. """ import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from scipy.stats import multivariate_normal import common from common import RED, BLUE, GREEN, MAGENTA, CYAN, ORANGE, PURPLE, sigmoid CH = 6 JITTER = 1e-6 # ======================================================================================== kernels def sqdist(X, Y): """Matrix of squared Euclidean distances between the rows of X and Y.""" X, Y = np.atleast_2d(X), np.atleast_2d(Y) return np.maximum((X ** 2).sum(1)[:, None] + (Y ** 2).sum(1)[None, :] - 2 * X @ Y.T, 0.0) def basis_kernel(Phi_x, Phi_y): """k(x, x') = phi(x)^T phi(x') (6.1), (6.10), for feature matrices with one row per input.""" return Phi_x @ Phi_y.T def k_poly(X, Y, M=2, c=0.0): """Polynomial kernel (x^T x' + c)^M, (6.12) (c = 0 is the homogeneous case with M = 2).""" return (np.atleast_2d(X) @ np.atleast_2d(Y).T + c) ** M def k_gauss(X, Y, sigma=1.0): """Gaussian kernel exp(-||x - x'||^2 / 2 sigma^2), (6.23).""" return np.exp(-sqdist(X, Y) / (2 * sigma ** 2)) def k_sigmoid(X, Y, a=1.0, b=0.0): """Sigmoidal kernel tanh(a x^T x' + b), (6.37) (not positive semidefinite in general).""" return np.tanh(a * np.atleast_2d(X) @ np.atleast_2d(Y).T + b) def k_exp(X, Y, theta=1.0): """Exponential (Ornstein-Uhlenbeck) kernel exp(-theta |x - x'|), (6.56).""" return np.exp(-theta * np.sqrt(sqdist(X, Y))) class GaussKernel: """Kernel (6.63)/(6.71)-(6.72): k = th0 exp(-1/2 sum_i eta_i (x_i - x'_i)^2) + th2 + th3 x^T x'. The free parameters are ln th0 and ln eta_i (one ln eta per input dimension, or a single shared one when tie=True, which is (6.63) with eta = theta_1); th2 and th3 are fixed constants. The derivatives dK/d(ln parameter) are analytic and used by gp_log_marginal and laplace_evidence. Note (6.63) writes theta_1/2 ||x - x'||^2, so eta = theta_1. """ def __init__(self, th0=1.0, eta=1.0, D=1, tie=False, th2=0.0, th3=0.0): eta = np.broadcast_to(np.asarray(eta, float), (1 if tie else D,)).copy() self.lp = np.concatenate([[np.log(th0)], np.log(eta)]) self.D, self.tie, self.th2, self.th3 = D, tie, th2, th3 def eta(self, lp=None): lp = self.lp if lp is None else lp e = np.exp(lp[1:]) return np.full(self.D, e[0]) if self.tie else e def K(self, X, Y=None, lp=None): lp = self.lp if lp is None else lp X = np.atleast_2d(X) Y = X if Y is None else np.atleast_2d(Y) e = self.eta(lp) d2 = sqdist(X * np.sqrt(e), Y * np.sqrt(e)) return np.exp(lp[0]) * np.exp(-0.5 * d2) + self.th2 + self.th3 * X @ Y.T def dK(self, X, lp=None): """List of dK/d(ln parameter) at the training inputs X, in the order of lp.""" lp = self.lp if lp is None else lp X = np.atleast_2d(X) e = self.eta(lp) E = np.exp(lp[0]) * np.exp(-0.5 * sqdist(X * np.sqrt(e), X * np.sqrt(e))) out = [E] if self.tie: out.append(E * (-0.5 * e[0] * sqdist(X, X))) else: for i in range(self.D): out.append(E * (-0.5 * e[i] * (X[:, i][:, None] - X[:, i][None, :]) ** 2)) return out def with_params(self, lp): k = GaussKernel.__new__(GaussKernel) k.__dict__.update(self.__dict__) k.lp = np.array(lp, float) return k class ExpKernel: """Exponential kernel th0 exp(-theta ||x - x'||), (6.56); parameters ln th0 and ln theta.""" def __init__(self, th0=1.0, theta=1.0): self.lp = np.log([th0, theta]) def K(self, X, Y=None, lp=None): lp = self.lp if lp is None else lp X = np.atleast_2d(X) Y = X if Y is None else np.atleast_2d(Y) return np.exp(lp[0]) * np.exp(-np.exp(lp[1]) * np.sqrt(sqdist(X, Y))) def dK(self, X, lp=None): lp = self.lp if lp is None else lp X = np.atleast_2d(X) d = np.sqrt(sqdist(X, X)) E = np.exp(lp[0]) * np.exp(-np.exp(lp[1]) * d) return [E, E * (-np.exp(lp[1]) * d)] def with_params(self, lp): k = ExpKernel.__new__(ExpKernel) k.__dict__.update(self.__dict__) k.lp = np.array(lp, float) return k # ======================================================================================== Section 6.1 def ridge_primal(Phi, t, lam): """Regularised least squares in the primal: w = (Phi^T Phi + lam I)^(-1) Phi^T t, from (6.2).""" return np.linalg.solve(Phi.T @ Phi + lam * np.eye(Phi.shape[1]), Phi.T @ t) def ridge_dual(K, t, lam): """Dual coefficients a = (K + lam I_N)^(-1) t, (6.8); predictions are k(x)^T a, (6.9).""" return np.linalg.solve(K + lam * np.eye(len(t)), t) # ======================================================================================== Section 6.3 def nadaraya_watson(x, X, t, sigma): """Nadaraya-Watson regression (6.45)-(6.48) with an isotropic Gaussian f(x, t) of sd sigma. Returns the weights k(x, x_n) (6.46) (each row sums to one), the conditional mean (6.45) and the conditional variance of the Gaussian mixture (6.48): sigma^2 + sum_n k(x, x_n) t_n^2 - (sum_n k(x, x_n) t_n)^2. """ x = np.atleast_1d(x).astype(float) W = np.exp(-(x[:, None] - X[None, :]) ** 2 / (2 * sigma ** 2)) W /= W.sum(1, keepdims=True) m = W @ t v = sigma ** 2 + W @ t ** 2 - m ** 2 return W, m, v # ======================================================================================== Section 6.4 def gp_predict(kern, X, t, Xs, beta): """Predictive mean and covariance of GP regression: (6.66)-(6.67) with C_N = K + beta^(-1) I. Returns m(x*) = k^T C^(-1) t and the noise-free variance c - k^T C^(-1) k, c = k(x*, x*); add 1/beta for the predictive variance of a new target. Also returns the full covariance of the function values at Xs (needed for posterior samples). """ C = kern.K(X) + np.eye(len(X)) / beta ks = kern.K(X, Xs) L = np.linalg.cholesky(C) alpha = np.linalg.solve(L.T, np.linalg.solve(L, t)) V = np.linalg.solve(L, ks) cov = kern.K(Xs) - V.T @ V return ks.T @ alpha, np.diag(cov).copy(), cov def gp_log_marginal(kern, lp, lnbeta, X, t, grad=True): """ln p(t | theta) = -1/2 ln|C| - 1/2 t^T C^(-1) t - N/2 ln 2 pi, (6.69), and its gradient (6.70). The parameters are lp (the kernel's ln parameters) followed by ln beta; the gradient is with respect to these logarithms: dC/d ln beta = -beta^(-1) I. """ N = len(t) beta = np.exp(lnbeta) C = kern.K(X, lp=lp) + np.eye(N) / beta L = np.linalg.cholesky(C) alpha = np.linalg.solve(L.T, np.linalg.solve(L, t)) val = -np.log(np.diag(L)).sum() - 0.5 * t @ alpha - 0.5 * N * np.log(2 * np.pi) if not grad: return val Ci = np.linalg.solve(L.T, np.linalg.solve(L, np.eye(N))) dCs = kern.dK(X, lp) + [-np.eye(N) / beta] g = [-0.5 * np.sum(Ci * dC) + 0.5 * alpha @ dC @ alpha for dC in dCs] # (6.70) return val, np.array(g) def gp_fit(kern, X, t, beta0=10.0, maxiter=100, fit_noise=True): """Maximise ln p(t | theta) (6.69) over the ln parameters and ln beta with L-BFGS-B. Returns the fitted kernel, beta and the trace of the parameters after each iteration (rows: iterations, columns: kernel ln parameters then ln beta). """ nk = len(kern.lp) x0 = np.concatenate([kern.lp, [np.log(beta0)]]) trace = [x0.copy()] def f(p): pp = p.copy() if not fit_noise: pp[-1] = x0[-1] try: v, g = gp_log_marginal(kern, pp[:nk], pp[nk], X, t) except np.linalg.LinAlgError: return 1e10, np.zeros_like(p) if not fit_noise: g[-1] = 0 return -v, -g res = minimize(f, x0, jac=True, method="L-BFGS-B", bounds=[(-12, 12)] * len(x0), callback=lambda p: trace.append(p.copy()), options=dict(maxiter=maxiter, ftol=1e-14, gtol=1e-9)) p = res.x.copy() if not fit_noise: p[-1] = x0[-1] return kern.with_params(p[:nk]), float(np.exp(p[nk])), np.array(trace) def laplace_mode(C, t, iters=50, tol=1e-12, damped=False): """Mode a*_N of p(a_N | t_N) by the Newton-Raphson iteration (6.83). a_new = C (I + W C)^(-1) (t - sigma + W a), with W = diag(sigma (1 - sigma)) (6.82). Returns a*, and the history of the objective Psi(a) (6.80, without constants) and of the gradient norm; with damped=True a step is halved until Psi does not decrease. """ N = len(t) a = np.zeros(N) hist = [psi(a, C, t)] gnorm = [np.linalg.norm(grad_psi(a, C, t))] for _ in range(iters): s = sigmoid(a) W = s * (1 - s) anew = C @ np.linalg.solve(np.eye(N) + W[:, None] * C, t - s + W * a) if damped: step = 1.0 while psi(a + step * (anew - a), C, t) < hist[-1] - 1e-14 and step > 1e-8: step /= 2 anew = a + step * (anew - a) a = anew hist.append(psi(a, C, t)) gnorm.append(np.linalg.norm(grad_psi(a, C, t))) if gnorm[-1] < tol: break return a, np.array(hist), np.array(gnorm) def psi(a, C, t): """Psi(a) = ln p(a) + ln p(t | a) without the constants not depending on a: -1/2 a^T C^(-1) a + t^T a - sum ln(1 + e^a) (6.80); the term -1/2 ln|C| is a constant too.""" return -0.5 * a @ np.linalg.solve(C, a) + t @ a - np.sum(np.logaddexp(0, a)) def grad_psi(a, C, t): """Gradient of Psi: t - sigma(a) - C^(-1) a, (6.81).""" return t - sigmoid(a) - np.linalg.solve(C, a) def laplace_predict(kern, X, t, Xs, nu=1e-4, lp=None, a=None, probit=True): """GP classification with the Laplace approximation (Section 6.4.6). Returns the predictive mean E[a*|t] = k^T (t - sigma) (6.87), the variance var[a*|t] = c - k^T (W^(-1) + C)^(-1) k (6.88), and p(t* = 1 | t) from the probit approximation (4.153), sigma(kappa(var) mean) with kappa = (1 + pi var / 8)^(-1/2). The covariance is C = K + nu I (6.75). """ N = len(t) C = kern.K(X, lp=lp) + nu * np.eye(N) if a is None: a, _, _ = laplace_mode(C, t) s = sigmoid(a) W = s * (1 - s) ks = kern.K(X, Xs, lp=lp) mean = ks.T @ (t - s) sw = np.sqrt(W) B = np.eye(N) + sw[:, None] * C * sw[None, :] V = np.linalg.solve(np.linalg.cholesky(B), sw[:, None] * ks) cs = np.diag(kern.K(Xs, lp=lp)) + nu var = cs - (V ** 2).sum(0) # (W^-1 + C)^-1 = W^1/2 B^-1 W^1/2 p = sigmoid(mean / np.sqrt(1 + np.pi * var / 8)) if probit else sigmoid(mean) return mean, var, p def laplace_evidence(kern, lp, X, t, nu=1e-4, grad=True): """Laplace approximation to ln p(t | theta), (6.90), and its gradient (6.91)-(6.94). ln p(t|theta) = Psi(a*) - 1/2 ln|W + C^(-1)| + N/2 ln 2 pi, with Psi including -1/2 ln|C| - N/2 ln 2 pi (the book's (6.80); the erratum removes a stray +const). Gradient with respect to the kernel's ln parameters: explicit part (6.91), implicit part through a*_N (6.92) with da*/dtheta from (6.94). Note on (6.93)-(6.94): differentiating a* = C (t - sigma(a*)) gives (I + C W) da* = dC (t - sigma), so da*/dtheta = (I + C W)^(-1) dC (t - sigma). The printed (6.94) has (I + W C)^(-1), which is not the same matrix (a transposition slip; it is not in the published errata). The finite-difference check in checks_ch06.py confirms the corrected form. """ N = len(t) C = kern.K(X, lp=lp) + nu * np.eye(N) a, _, _ = laplace_mode(C, t) s = sigmoid(a) W = s * (1 - s) sw = np.sqrt(W) B = np.eye(N) + sw[:, None] * C * sw[None, :] LB = np.linalg.cholesky(B) LC = np.linalg.cholesky(C) Ci_a = np.linalg.solve(LC.T, np.linalg.solve(LC, a)) lnC = 2 * np.log(np.diag(LC)).sum() psi_star = -0.5 * a @ Ci_a - 0.5 * N * np.log(2 * np.pi) - 0.5 * lnC + t @ a - np.sum(np.logaddexp(0, a)) # |W + C^-1| = |C^-1| |I + C W| = |I + W^1/2 C W^1/2| / |C| ln_det_H = 2 * np.log(np.diag(LB)).sum() - lnC val = psi_star - 0.5 * ln_det_H + 0.5 * N * np.log(2 * np.pi) if not grad: return val # (I + C W)^-1 W = W^1/2 B^-1 W^1/2 ; (I + W C)^-1 = I - W^1/2 B^-1 W^1/2 C ... use explicit solves Bi = np.linalg.solve(B, np.eye(N)) R = sw[:, None] * Bi * sw[None, :] # (I + W C)^-1 W = W^1/2 B^-1 W^1/2 (symmetric) ICW = np.eye(N) + C * W[None, :] # I + C W ICW_inv_C = np.linalg.solve(ICW, C) # (I + C W)^-1 C g = [] for dC in kern.dK(X, lp): explicit = 0.5 * Ci_a @ dC @ Ci_a - 0.5 * np.trace(R @ dC) # (6.91) da = np.linalg.solve(ICW, dC @ (t - s)) # (6.94), see the note below implicit = -0.5 * np.sum(np.diag(ICW_inv_C) * s * (1 - s) * (1 - 2 * s) * da) # (6.92) g.append(explicit + implicit) return val, np.array(g) def gpc_fit(kern, X, t, nu=1e-4, maxiter=60): """Learn the kernel parameters of a GP classifier by maximising the Laplace evidence (6.90).""" nk = len(kern.lp) def f(p): try: v, g = laplace_evidence(kern, p, X, t, nu) except np.linalg.LinAlgError: return 1e10, np.zeros(nk) return -v, -g res = minimize(f, kern.lp.copy(), jac=True, method="L-BFGS-B", bounds=[(-8, 8)] * nk, options=dict(maxiter=maxiter)) return kern.with_params(res.x), -res.fun # ======================================================================================== data sets def sinusoid(rng, N, noise=0.3, x=None): """The synthetic regression set of Appendix A: t = sin(2 pi x) + N(0, noise^2).""" x = np.sort(rng.uniform(0, 1, N)) if x is None else x return x, np.sin(2 * np.pi * x) + noise * rng.standard_normal(len(x)) # Two-class synthetic data in the style of Appendix A (Figure A.7): the blue class comes from one # Gaussian, the red class from a mixture of two; equal priors. CLASS_BLUE = dict(w=[1.0], mu=[[-0.4, 0.1]], cov=[np.diag([0.55, 1.1])]) CLASS_RED = dict(w=[0.5, 0.5], mu=[[1.0, 1.0], [1.1, -0.9]], cov=[np.diag([0.35, 0.45]), np.diag([0.45, 0.45])]) def class_density(P, X): return sum(w * multivariate_normal(m, c).pdf(X) for w, m, c in zip(P["w"], P["mu"], P["cov"])) def two_class_data(rng, n_per_class=100): """Draw n points per class from CLASS_BLUE (t = 0) and CLASS_RED (t = 1); returns X, t.""" def draw(P): k = rng.choice(len(P["w"]), size=n_per_class, p=P["w"]) return np.array([rng.multivariate_normal(P["mu"][j], P["cov"][j]) for j in k]) X = np.vstack([draw(CLASS_BLUE), draw(CLASS_RED)]) return X, np.r_[np.zeros(n_per_class), np.ones(n_per_class)] def true_posterior_red(X): """p(red | x) with equal priors (the colour scale of Figure A.7).""" pr, pb = class_density(CLASS_RED, X), class_density(CLASS_BLUE, X) return pr / (pr + pb) def prior_samples(K, z, jitter=JITTER): """Draws L z from N(0, K), with the standard normal draws z (columns) shared between panels.""" L = np.linalg.cholesky(K + jitter * np.eye(len(K))) return L @ z # ======================================================================================== figures def _box(ax, xl=None, yl=None): if xl is not None: ax.set_xlim(*xl) if yl is not None: ax.set_ylim(*yl) def fig6_1(out_dir=None): """Figure 6.1: kernels built from basis functions: polynomial, Gaussian, sigmoid (6.10).""" x = np.linspace(-1, 1, 400) cols = [RED, GREEN, BLUE, MAGENTA, CYAN, "#b5a800", "#777777"] sig = 0.1 bases = { "poly": (lambda u: np.stack([u ** j for j in range(1, 11)], 1), range(10), -0.5), "gauss": (lambda u: np.stack([np.exp(-(u - m) ** 2 / (2 * sig ** 2)) for m in np.linspace(-1, 1, 21)], 1), range(0, 21, 2), 0.0), "sigm": (lambda u: np.stack([sigmoid((u - m) / sig) for m in np.linspace(-1, 1, 11)], 1), range(11), 0.0), } fig, axes = common.panel_grid(2, 3, size=(3.1, 2.5)) for c, (name, (phi, show, x0)) in enumerate(bases.items()): P = phi(x) ax = axes[0, c] for r, j in enumerate(show): ax.plot(x, P[:, j], color=cols[r % len(cols)], lw=1.2) ax.set_xlim(-1, 1) k = basis_kernel(P, phi(np.array([x0])))[:, 0] ax = axes[1, c] ax.plot(x, k, color=BLUE, lw=2) ax.plot([x0], [ax.get_ylim()[0]], "x", color=RED, ms=7, mew=2, clip_on=False) ax.set_xlim(-1, 1) if name == "poly": ax.axhline(0, color="k", lw=0.6, ls=":") print(" fig 6.1 %s: max of k = %.3f, min = %.3f" % (name, k.max(), k.min())) fig.tight_layout() return common.savefig(fig, CH, "fig6-1", out_dir) def fig6_2(out_dir=None): """Figure 6.2: Gaussian basis functions and the normalised basis functions (6.41).""" x = np.linspace(-1, 1, 500) mu, s = np.array([-0.5, 0.1, 0.65]), np.array([0.1, 0.13, 0.25]) P = np.exp(-(x[:, None] - mu) ** 2 / (2 * s ** 2)) fig, axes = common.panel_grid(1, 2, size=(4.2, 3.0)) for ax, Q in zip(axes[0], (P, P / P.sum(1, keepdims=True))): for j, c in enumerate((BLUE, GREEN, RED)): ax.plot(x, Q[:, j], color=c) ax.set_xlim(-1, 1) ax.set_ylim(0, 1.02) fig.tight_layout() return common.savefig(fig, CH, "fig6-2", out_dir) def fig6_3(out_dir=None): """Figure 6.3: Nadaraya-Watson regression with isotropic Gaussian kernels (6.45)-(6.48).""" rng = np.random.default_rng(3) X, t = sinusoid(rng, 10, 0.3, x=np.linspace(0, 1, 10)) sig = 0.1 xs = np.linspace(0, 1, 400) W, m, v = nadaraya_watson(xs, X, t, sig) fig, ax = plt.subplots(figsize=(5.0, 3.6)) ax.fill_between(xs, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v), color=RED, alpha=0.15, lw=0) ax.plot(xs, np.sin(2 * np.pi * xs), color=GREEN, lw=2) ax.plot(xs, m, color=RED, lw=2) ax.plot(X, t, "o", color=BLUE, ms=4, zorder=5) th = np.linspace(0, 2 * np.pi, 80) for xn, tn in zip(X, t): ax.plot(xn + sig * np.cos(th), tn + sig * np.sin(th), color=BLUE, lw=1) ax.set_xlim(-0.08, 1.08) ax.set_ylim(-1.6, 1.6) print(" fig 6.3: weights sum to 1 (min %.6f, max %.6f)" % (W.sum(1).min(), W.sum(1).max())) return common.savefig(fig, CH, "fig6-3", out_dir) def fig6_4(out_dir=None): """Figure 6.4: samples from a GP with a Gaussian kernel (left) and an exponential kernel (right).""" x = np.linspace(-1, 1, 300)[:, None] z = np.random.default_rng(1).standard_normal((len(x), 5)) cols = [RED, BLUE, "#7b1fa2", CYAN, "#0a7a0a"] fig, axes = common.panel_grid(1, 2, size=(4.2, 3.0)) for ax, K in zip(axes[0], (1.0 * np.exp(-0.5 * 4 * sqdist(x, x)), k_exp(x, x, 4.0))): Y = prior_samples(K, z) for j in range(5): ax.plot(x[:, 0], Y[:, j], color=cols[j], lw=1.6) ax.set_xlim(-1, 1) ax.set_ylim(-3, 3) fig.tight_layout() return common.savefig(fig, CH, "fig6-4", out_dir) def fig6_5(out_dir=None): """Figure 6.5: GP prior samples for the covariance (6.63), (theta0..theta3) as in the panel titles.""" x = np.linspace(-1, 1, 300)[:, None] z = np.random.default_rng(1).standard_normal((len(x), 5)) cols = [RED, BLUE, "#7b1fa2", CYAN, "#0a7a0a"] thetas = [(1, 4, 0, 0), (9, 4, 0, 0), (1, 64, 0, 0), (1, 0.25, 0, 0), (1, 4, 10, 0), (1, 4, 0, 5)] fig, axes = common.panel_grid(2, 3, size=(3.3, 2.7)) ylims = [3, 9, 3, 3, 9, 4] for ax, th, yl in zip(axes.ravel(), thetas, ylims): K = th[0] * np.exp(-0.5 * th[1] * sqdist(x, x)) + th[2] + th[3] * x @ x.T # (6.63) Y = prior_samples(K, z) for j in range(5): ax.plot(x[:, 0], Y[:, j], color=cols[j], lw=1.5) ax.set_xlim(-1, 1) ax.set_ylim(-yl, yl) ax.set_title("(%.2f, %.2f, %.2f, %.2f)" % th, fontsize=9) fig.tight_layout() return common.savefig(fig, CH, "fig6-5", out_dir) def fig6_6(out_dir=None): """Figure 6.6: a function drawn from the GP prior, its values y_n at inputs x_n, and the noisy t_n (6.57)-(6.62).""" rng = np.random.default_rng(6) xg = np.linspace(-1, 1, 300) xn = np.sort(rng.uniform(-0.95, 0.95, 12)) xa = np.concatenate([xg, xn]) K = np.exp(-0.5 * 4 * sqdist(xa[:, None], xa[:, None])) y = prior_samples(K, rng.standard_normal((len(xa), 1)))[:, 0] tn = y[300:] + 0.5 * rng.standard_normal(12) fig, ax = plt.subplots(figsize=(4.8, 3.3)) ax.plot(xg, y[:300], color=BLUE, lw=1.8) for a, b, c in zip(xn, y[300:], tn): ax.plot([a, a], [b, c], color="#b5a800", lw=1.3) ax.plot(xn, y[300:], "o", mfc="none", mec=RED, ms=6, mew=1.4) ax.plot(xn, tn, "o", mfc="none", mec=GREEN, ms=6, mew=1.4) ax.set_xlim(-1, 1) ax.set_ylim(-3, 3) ax.set_xlabel("$x$") ax.set_ylabel("$t$", rotation=0) return common.savefig(fig, CH, "fig6-6", out_dir) def fig6_7(out_dir=None): """Figure 6.7: one training point and one test point: the joint p(t1, t2) and the conditional p(t2 | t1) (6.66)-(6.67).""" rho, t1 = 0.5, 0.95 C = np.array([[1, rho], [rho, 1]]) g = np.linspace(-2, 2, 300) T1, T2 = np.meshgrid(g, g) pdf = multivariate_normal(np.zeros(2), C).pdf(np.dstack([T1, T2])) peak = 1 / (2 * np.pi * np.sqrt(np.linalg.det(C))) fig, ax = plt.subplots(figsize=(4.4, 4.0)) ax.contour(T1, T2, pdf, levels=np.sort(peak * np.exp(-0.5 * np.array([0.5, 1.5, 3.0, 5.0]))), colors=RED, linewidths=1.1) m, var = rho * t1, 1 - rho ** 2 cond = np.exp(-0.5 * (g - m) ** 2 / var) / np.sqrt(2 * np.pi * var) ax.plot(0.9 * cond, g, color=GREEN, lw=2) ax.axvline(t1, color=BLUE, lw=1) ax.axhline(m, color=GREEN, lw=0.9, xmax=(t1 + 2) / 4) ax.plot(t1, 0, "o", mfc="none", mec=BLUE, ms=6, mew=1.5, zorder=5) ax.plot(0, m, "o", mfc="none", mec=GREEN, ms=6, mew=1.5, zorder=5) ax.plot([-2, 2], [-2, 2], "k:", lw=0.8) ax.axhline(0, color="k", lw=0.7) ax.axvline(0, color="k", lw=0.7) ax.text(0.06, m + 0.1, "$m(x_2)$", fontsize=9, ha="left") ax.set_xlabel("$t_1$") ax.set_ylabel("$t_2$", rotation=0) ax.set_xlim(-2, 2) ax.set_ylim(-2, 2) ax.set_aspect("equal") print(" fig 6.7: conditional mean %.3f = rho t1, variance %.3f = 1 - rho^2" % (m, var)) return common.savefig(fig, CH, "fig6-7", out_dir) def fig6_8(out_dir=None): """Figure 6.8: GP regression on the sinusoidal data with the three right-most points omitted.""" rng = np.random.default_rng(8) x, t = sinusoid(rng, 10, 0.3, x=np.linspace(0, 1, 10)) x, t = x[:7], t[:7] kern = GaussKernel(1.0, 12.0, D=1, tie=True) xs = np.linspace(0, 1, 300)[:, None] m, v, _ = gp_predict(kern, x[:, None], t, xs, 1 / 0.3 ** 2) sd = np.sqrt(v + 0.3 ** 2) fig, ax = plt.subplots(figsize=(4.8, 3.3)) ax.fill_between(xs[:, 0], m - 2 * sd, m + 2 * sd, color=RED, alpha=0.15, lw=0) ax.plot(xs[:, 0], np.sin(2 * np.pi * xs[:, 0]), color=GREEN, lw=2) ax.plot(xs[:, 0], m, color=RED, lw=2) ax.plot(x, t, "o", mfc="none", mec=BLUE, ms=6, mew=1.5) ax.set_xlim(-0.05, 1.05) ax.set_ylim(-1.45, 1.4) print(" fig 6.8: predictive sd at x = 0.3 and 1.0: %.2f, %.2f" % (sd[90], sd[-1])) return common.savefig(fig, CH, "fig6-8", out_dir) def fig6_9(out_dir=None): """Figure 6.9: samples from the ARD prior (6.71): eta = (1, 1) (left) and (1, 0.01) (right).""" n = 36 g = np.linspace(-4, 4, n) G1, G2 = np.meshgrid(g, g) P = np.c_[G1.ravel(), G2.ravel()] z = np.random.default_rng(9).standard_normal((len(P), 1)) fig = plt.figure(figsize=(8.4, 4.0)) for c, eta in enumerate(((1.0, 1.0), (1.0, 0.01))): K = GaussKernel(1.0, eta, D=2).K(P) Y = prior_samples(K, z, 1e-6)[:, 0].reshape(n, n) ax = fig.add_subplot(1, 2, c + 1, projection="3d") ax.plot_surface(G1, G2, Y, cmap="copper", linewidth=0, antialiased=True, rasterized=True) ax.set_xlabel("$x_1$", labelpad=-4) ax.set_ylabel("$x_2$", labelpad=-4) ax.set_zlim(-3, 3) ax.tick_params(labelsize=7, pad=-2) ax.set_title("$\\eta_1 = 1,\\ \\eta_2 = %s$" % ("1" if c == 0 else "0.01"), fontsize=9) ax.view_init(28, -55) fig.tight_layout() return common.savefig(fig, CH, "fig6-9", out_dir) def ard_data(rng, N=100): """The ARD data set of Section 6.4.4: x1 Gaussian, t = sin(2 pi x1) + noise, x2 = x1 + noise, x3 independent.""" x1 = 0.35 * rng.standard_normal(N) t = np.sin(2 * np.pi * x1) + 0.1 * rng.standard_normal(N) x2 = x1 + 0.05 * rng.standard_normal(N) x3 = 0.35 * rng.standard_normal(N) return np.c_[x1, x2, x3], t def fig6_10(out_dir=None): """Figure 6.10: ARD: the hyperparameters eta_1, eta_2, eta_3 against the iteration number of the marginal-likelihood optimisation.""" X, t = ard_data(np.random.default_rng(10)) kern = GaussKernel(1.0, 1.0, D=3) kf, beta, tr = gp_fit(kern, X, t, beta0=10.0, maxiter=100) eta = np.exp(tr[:, 1:4]) fig, ax = plt.subplots(figsize=(4.8, 3.3)) for j, c in enumerate((RED, GREEN, BLUE)): ax.semilogy(np.arange(len(eta)), eta[:, j], color=c, lw=1.6) ax.set_xlim(0, 100) ax.set_xlabel("iteration") ax.set_ylabel("$\\eta_i$") ax.legend(["$\\eta_1$", "$\\eta_2$", "$\\eta_3$"], frameon=False, loc="center right") print(" fig 6.10: final eta = %s after %d iterations, beta = %.1f" % (np.round(eta[-1], 4), len(eta) - 1, beta)) return common.savefig(fig, CH, "fig6-10", out_dir) def fig6_11(out_dir=None): """Figure 6.11: a sample a(x) from a GP prior (left) and the sigmoid-transformed sample sigma(a(x)) (right).""" x = np.linspace(-1, 1, 300)[:, None] K = 16 * np.exp(-0.5 * 16 * sqdist(x, x)) a = prior_samples(K, np.random.default_rng(11).standard_normal((300, 1)))[:, 0] fig, axes = common.panel_grid(1, 2, size=(4.2, 3.0)) axes[0, 0].plot(x[:, 0], a, color=RED, lw=1.8) axes[0, 1].plot(x[:, 0], sigmoid(a), color=RED, lw=1.8) for ax in axes[0]: ax.set_xlim(-1, 1) axes[0, 1].set_ylim(0, 1) fig.tight_layout() return common.savefig(fig, CH, "fig6-11", out_dir) def fig6_12(out_dir=None): """Figure 6.12: GP classification (Laplace): data and decision boundaries (left); predicted probability (right).""" rng = np.random.default_rng(12) X, t = two_class_data(rng, 100) kern, ev = gpc_fit(GaussKernel(4.0, 1.0, D=2, tie=True), X, t) g1, g2 = np.linspace(-2.6, 2.6, 80), np.linspace(-2.8, 3.2, 90) G1, G2 = np.meshgrid(g1, g2) P = np.c_[G1.ravel(), G2.ravel()] _, _, p = laplace_predict(kern, X, t, P) p = p.reshape(G1.shape) ptrue = true_posterior_red(P).reshape(G1.shape) err_gp = np.mean((laplace_predict(kern, X, t, X)[2] > 0.5) != (t == 1)) # test error against the true boundary on a large fresh sample Xt, tt = two_class_data(np.random.default_rng(112), 2000) pt = laplace_predict(kern, X, t, Xt)[2] print(" fig 6.12: theta0 = %.2f, eta = %.2f, ln evidence %.2f; training error %.3f; test error GP %.3f, Bayes %.3f" % (np.exp(kern.lp[0]), np.exp(kern.lp[1]), ev, err_gp, np.mean((pt > 0.5) != (tt == 1)), np.mean((true_posterior_red(Xt) > 0.5) != (tt == 1)))) fig, axes = common.panel_grid(1, 2, size=(4.2, 3.9)) ax = axes[0, 0] ax.plot(*X[t == 0].T, "o", mfc="none", mec=BLUE, ms=5) ax.plot(*X[t == 1].T, "x", color=RED, ms=5) ax.contour(G1, G2, ptrue, levels=[0.5], colors=GREEN, linewidths=1.4) ax.contour(G1, G2, p, levels=[0.5], colors="k", linewidths=1.2) ax = axes[0, 1] ax.pcolormesh(G1, G2, p, cmap="coolwarm", vmin=0, vmax=1, shading="auto", rasterized=True) ax.contour(G1, G2, p, levels=[0.5], colors="k", linewidths=1.2) for ax in axes[0]: ax.set_xlim(-2.6, 2.6) ax.set_ylim(-2.8, 3.2) fig.tight_layout() return common.savefig(fig, CH, "fig6-12", out_dir) FIGURES = {"6.1": fig6_1, "6.2": fig6_2, "6.3": fig6_3, "6.4": fig6_4, "6.5": fig6_5, "6.6": fig6_6, "6.7": fig6_7, "6.8": fig6_8, "6.9": fig6_9, "6.10": fig6_10, "6.11": fig6_11, "6.12": fig6_12} if __name__ == "__main__": common.run(FIGURES, 6)