"""Numerical checks of the derivations in Bishop (2006), Chapter 4 (Linear Models for Classification).

Each check prints PASS or FAIL with the numbers behind it; the last line counts them.
Run from this folder: python3 checks_ch04.py (a few seconds). The algorithms checked are the ones in
ch04_linear_classification.py, which also produce the figures. Gradients and Hessians are compared
with finite differences, closed forms with numerical optimisation or quadrature.
"""
import math

import numpy as np
from scipy import integrate, optimize
from scipy.special import expit, ndtr

import common
import ch04_linear_classification as m

RESULTS = []


def check(name, ok, detail=""):
    RESULTS.append(bool(ok))
    print("%s  %s%s" % ("PASS" if ok else "FAIL", name, ("  [" + detail + "]") if detail else ""))


def fd_grad(f, x, eps=1e-6):
    g = np.zeros_like(x)
    for i in range(len(x)):
        e = np.zeros_like(x)
        e[i] = eps
        g[i] = (f(x + e) - f(x - e)) / (2 * eps)
    return g


rng = np.random.default_rng(4)

# ---------------------------------------------------------------------------------------- 4.1
X2, y2 = m.make_two_class(outliers=0)
Phi2 = m.design(X2)
T2 = m.one_hot(y2)
W = m.least_squares_classifier(Phi2, T2)
check("least squares: pseudo-inverse solution (4.16) satisfies the normal equations Phi^T (Phi W - T) = 0",
      np.abs(Phi2.T @ (Phi2 @ W - T2)).max() < 1e-9, "max |residual| = %.1e" % np.abs(Phi2.T @ (Phi2 @ W - T2)).max())
check("least squares with a bias: the K outputs sum to one for every x (Exercise 4.2)",
      np.allclose((Phi2 @ W).sum(1), 1.0) and np.allclose((m.design(rng.normal(size=(5, 2)) * 5) @ W).sum(1), 1.0))

# Fisher = LDA
Xf, yf = m.make_fisher_data()
wF, mk, SW = m.fisher_two_class(Xf, yf)
pi, mu, Sig = m.gaussian_generative_fit(Xf, yf)
w_lda = np.linalg.solve(Sig, mu[1] - mu[0])
cos = wF @ w_lda / np.linalg.norm(w_lda)
check("Fisher direction (4.30) equals the direction of the Gaussian generative (LDA) classifier (4.66)", abs(abs(cos) - 1) < 1e-12, "cosine = %.15f" % cos)
_, _, SB = m.class_means_scatter(Xf, yf)
Jbest = m.fisher_criterion(wF, SW, SB)
Jrand = max(m.fisher_criterion(r, SW, SB) for r in rng.normal(size=(20000, 2)))
ang = np.linspace(0, np.pi, 100001)
Jgrid = np.array([m.fisher_criterion(np.array([np.cos(a), np.sin(a)]), SW, SB) for a in ang[::50]])
check("Fisher criterion J(w) (4.26) is maximised by w_F: no random or grid direction beats it", Jbest >= Jrand - 1e-9 and Jbest >= Jgrid.max() - 1e-9,
      "J(w_F) = %.4f, best of 20000 random directions %.4f, best on a grid %.4f" % (Jbest, Jrand, Jgrid.max()))
check("S_B = (m_2 - m_1)(m_2 - m_1)^T N_1 N_2 / N for two classes, so S_B w is parallel to m_2 - m_1 (text before (4.30))",
      np.allclose(SB, np.outer(mk[1] - mk[0], mk[1] - mk[0]) * (yf == 0).sum() * (yf == 1).sum() / len(yf)))

# multiclass Fisher
X3, y3 = m.make_three_class()
mk3, SW3, SB3 = m.class_means_scatter(X3, y3)
ST3 = (X3 - X3.mean(0)).T @ (X3 - X3.mean(0))
check("total scatter S_T = S_W + S_B (4.45)", np.allclose(ST3, SW3 + SB3))
J = lambda Wm: float(np.trace(np.linalg.solve(Wm.T @ SW3 @ Wm, Wm.T @ SB3 @ Wm)))
Wf, vals = m.fisher_multiclass(X3, y3, 1)
Jr = max(J(rng.normal(size=(2, 1))) for _ in range(5000))
check("multiclass Fisher: for one direction J(w) = Tr{(W S_W W^T)^-1 (W S_B W^T)} is maximised by the top eigenvector of S_W^-1 S_B, J = the top eigenvalue",
      abs(J(Wf) - vals[0]) < 1e-8 and J(Wf) >= Jr - 1e-9, "J = %.4f, eigenvalues %s, best of 5000 random directions %.4f" % (J(Wf), np.round(vals, 4), Jr))
check("with all D = 2 directions (an invertible W) J equals the sum of the eigenvalues whatever W is", abs(J(rng.normal(size=(2, 2))) - vals.sum()) < 1e-8, "sum = %.4f" % vals.sum())
Z = np.hstack([X3, np.random.default_rng(1).normal(size=(len(X3), 2))])
_, v4 = m.fisher_multiclass(Z, y3, 2)
check("rank of S_B is at most K - 1: three classes in 4-D give two positive eigenvalues and two (numerically) zero ones", v4[1] > 1e-3 and abs(v4[2]) < 1e-9 * v4[0] and abs(v4[3]) < 1e-9 * v4[0],
      "eigenvalues %s" % np.array2string(v4, precision=4))

# perceptron
def separable(rng, n, margin_dir=None):
    wt = np.array([0.3, 1.0, -0.7])
    Xp = rng.uniform(-2, 2, size=(2000, 2))
    Ph = m.design(Xp)
    s = Ph @ wt / np.linalg.norm(wt)
    keep = np.abs(s) > 0.3
    Ph, s = Ph[keep][:n], s[keep][:n]
    return Ph, np.sign(s), wt / np.linalg.norm(wt)


ok, worst = True, 0.0
for trial in range(25):
    Ph, tt, wstar = separable(np.random.default_rng(100 + trial), 60)
    w, hist = m.perceptron(Ph, tt)
    gamma = (tt * (Ph @ wstar)).min()
    R = np.linalg.norm(Ph, axis=1).max()
    bound = (R / gamma) ** 2
    ok &= bool((tt * (Ph @ w) > 0).all()) and len(hist) <= bound
    worst = max(worst, len(hist) / bound)
check("perceptron converges on separable data within the bound (R / gamma)^2 (perceptron convergence theorem, Section 4.1.7)", ok,
      "25 data sets, largest number of updates / bound = %.3f" % worst)
Xp = np.vstack([np.array(m.PERCEPTRON_RED), np.array(m.PERCEPTRON_BLUE)])
tp = np.r_[np.ones(5), -np.ones(5)]
wp, hp = m.perceptron(Xp, tp, w0=m.PERCEPTRON_W0)
check("Figure 4.7 data: two updates (w <- w + t phi, (4.55)), the first on (0.23, 0.72), the second on (0.45, 0.30), final w = (0.26, 0.12)",
      len(hp) == 2 and hp[0][1] == 0 and hp[1][1] == 1 and np.allclose(wp, [0.26, 0.12]) and (tp * (Xp @ wp) > 0).all(), "w = %s" % np.round(wp, 3))
# perceptron criterion decreases at each update
EP = lambda w: -sum(tp[i] * (Xp[i] @ w) for i in range(10) if tp[i] * (Xp[i] @ w) <= 0)
check("each update changes the perceptron criterion (4.54) by -t phi^T t phi < 0 for the visited point (the point's own error falls)",
      all(-(tp[n] * (Xp[n] @ wa)) < -(tp[n] * (Xp[n] @ wb)) for wb, n, wa in hp))

# ---------------------------------------------------------------------------------------- 4.2
Xg, yg = m.make_two_class()
pi, mu, Sig = m.gaussian_generative_fit(Xg, yg)


def neg_ll(theta):
    p1 = expit(theta[0])
    mu0, mu1 = theta[1:3], theta[3:5]
    L = np.array([[math.exp(theta[5]), 0], [theta[6], math.exp(theta[7])]])
    S = L @ L.T
    ll = 0.0
    for k, (pk, mk_) in enumerate([(1 - p1, mu0), (p1, mu1)]):
        ll += np.sum(np.log(pk) + common.gauss_logpdf(Xg[yg == k], mk_, S))
    return -ll


L0 = np.linalg.cholesky(Sig)
theta0 = np.r_[0.0, Xg[yg == 0].mean(0) + 0.3, Xg[yg == 1].mean(0) - 0.3, 0.0, 0.0, 0.0, 0.0]
res = optimize.minimize(neg_ll, theta0, method="BFGS", options=dict(gtol=1e-8, maxiter=2000))
L = np.array([[math.exp(res.x[5]), 0], [res.x[6], math.exp(res.x[7])]])
check("Gaussian generative model: closed-form ML (4.75)-(4.80) agrees with numerical maximisation of the joint log likelihood (4.71)",
      abs(expit(res.x[0]) - pi[1]) < 1e-4 and np.abs(res.x[1:3] - mu[0]).max() < 1e-3 and np.abs(res.x[3:5] - mu[1]).max() < 1e-3 and np.abs(L @ L.T - Sig).max() < 1e-3,
      "pi_1: %.5f vs %.5f; max |Sigma diff| = %.1e" % (pi[1], expit(res.x[0]), np.abs(L @ L.T - Sig).max()))
check("ML covariance is the weighted average of the class covariances, Sigma = S_W / N for equal shares (4.79-4.80)",
      np.allclose(Sig, sum((yg == k).mean() * np.cov(Xg[yg == k].T, bias=True) for k in range(2))))
W2, w02 = m.generative_posterior_params(pi, mu, Sig)
w_, b_ = W2[0] - W2[1], w02[0] - w02[1]
Xt = rng.normal(size=(50, 2)) * 3
post_direct = np.array([pi[k] * np.exp(common.gauss_logpdf(Xt, mu[k], Sig)) for k in range(2)])
post_direct = post_direct / post_direct.sum(0)
check("posterior p(C_1|x) is a logistic sigmoid of w^T x + w_0 with w = Sigma^-1 (mu_1 - mu_2) and w_0 from (4.66)-(4.67)",
      np.allclose(post_direct[0], expit(Xt @ w_ + b_), atol=1e-12), "max diff = %.1e" % np.abs(post_direct[0] - expit(Xt @ w_ + b_)).max())
pi3, mu3, Sig3 = m.gaussian_generative_fit(X3, y3)
W3, w03 = m.generative_posterior_params(pi3, mu3, Sig3)
post3 = np.array([pi3[k] * np.exp(common.gauss_logpdf(Xt, mu3[k], Sig3)) for k in range(3)])
post3 = (post3 / post3.sum(0)).T
check("K classes: posterior is softmax(a_k) with a_k = w_k^T x + w_k0, w_k = Sigma^-1 mu_k (4.62)-(4.65)",
      np.allclose(post3, m.softmax(Xt @ W3.T + w03), atol=1e-12))
a = np.linspace(-6, 6, 13)
check("sigmoid identities: sigma(-a) = 1 - sigma(a) (4.60), d sigma/da = sigma (1 - sigma) (4.88), logit sigma^-1(y) = ln(y/(1-y)) (4.61)",
      np.allclose(expit(-a), 1 - expit(a)) and np.allclose((expit(a + 1e-6) - expit(a - 1e-6)) / 2e-6, expit(a) * (1 - expit(a)), atol=1e-9)
      and np.allclose(np.log(expit(a) / (1 - expit(a))), a))

# ---------------------------------------------------------------------------------------- 4.3
Xl, yl = m.blobs(m.JSRandom(7), 25, 0.8)
Phl = m.design(Xl)
tl = yl.astype(float)
w0_ = rng.normal(size=3)
for alpha in (0.0, 0.7):
    g, H = m.logistic_grad_hess(w0_, Phl, tl, alpha)
    gfd = fd_grad(lambda w: -m.logistic_loglik(w, Phl, tl, alpha), w0_)
    Hfd = np.array([fd_grad(lambda w, i=i: m.logistic_grad_hess(w, Phl, tl, alpha)[0][i], w0_) for i in range(3)])
    check("logistic gradient (4.96) and Hessian (4.97) match finite differences, alpha = %.1f" % alpha,
          np.abs(g - gfd).max() < 1e-6 and np.abs(H - Hfd).max() < 1e-5, "grad err %.1e, Hessian err %.1e" % (np.abs(g - gfd).max(), np.abs(H - Hfd).max()))
w, info = m.irls_logistic(Phl, tl)
dll = np.diff(info["loglik"])
g, _ = m.logistic_grad_hess(w, Phl, tl)
check("IRLS: the log likelihood never decreases, the gradient vanishes at the end", (dll >= -1e-12).all() and np.abs(g).max() < 1e-8,
      "%d iterations, ln L from %.3f to %.4f, |grad| = %.1e" % (info["iters"], info["loglik"][0], info["loglik"][-1], np.abs(g).max()))
res = optimize.minimize(lambda w: -m.logistic_loglik(w, Phl, tl), np.zeros(3), jac=lambda w: m.logistic_grad_hess(w, Phl, tl)[0], method="BFGS", options=dict(gtol=1e-10))
check("IRLS solution equals the BFGS solution (the error function is convex, Exercise 4.15)", np.abs(w - res.x).max() < 1e-5, "max |diff| = %.1e" % np.abs(w - res.x).max())
# weighted least-squares form (4.99)-(4.101)
wt = rng.normal(size=3)
y = expit(Phl @ wt)
R = y * (1 - y)
z = Phl @ wt - (y - tl) / R
w_new = np.linalg.solve(Phl.T @ (Phl * R[:, None]), Phl.T @ (R * z))
g, H = m.logistic_grad_hess(wt, Phl, tl)
check("one Newton-Raphson step equals the weighted least-squares solution w = (Phi^T R Phi)^-1 Phi^T R z (4.99)-(4.101)",
      np.allclose(w_new, wt - np.linalg.solve(H, g), atol=1e-10))
Xs, ys = m.blobs(m.JSRandom(3), 20, 3.0, 0.4)
Pss = m.design(Xs)
_, infos = m.irls_logistic(Pss, ys.astype(float), max_iter=30)
nrm = np.linalg.norm(infos["path"], axis=1)
check("linearly separable data, no prior: |w| grows at every IRLS step and the log likelihood tends to 0 (Section 4.3.2)",
      (np.diff(nrm) > 0).all() and infos["loglik"][-1] > -1e-6, "|w| after 5, 15, 30 steps: %.1f, %.1f, %.1f; ln L = %.1e" % (nrm[5], nrm[15], nrm[30], infos["loglik"][-1]))
w_reg, _ = m.irls_logistic(Pss, ys.astype(float), alpha=0.1)
check("with a Gaussian prior (alpha > 0) the same data give a finite w (MAP)", np.linalg.norm(w_reg) < 20, "|w_MAP| = %.2f" % np.linalg.norm(w_reg))

# multiclass softmax
A = rng.normal(size=(4, 3))
Y = m.softmax(A, axis=1)
Jac = np.array([[(Y[0, k] * ((k == j) - Y[0, j])) for j in range(3)] for k in range(3)])
fdJ = np.array([[(m.softmax(A + np.eye(3)[j][None] * 1e-6 * (np.arange(4) == 0)[:, None], axis=1)[0, k]
                  - m.softmax(A - np.eye(3)[j][None] * 1e-6 * (np.arange(4) == 0)[:, None], axis=1)[0, k]) / 2e-6 for j in range(3)] for k in range(3)])
check("softmax derivative d y_k / d a_j = y_k (I_kj - y_j) (4.106)", np.abs(Jac - fdJ).max() < 1e-9, "max err %.1e" % np.abs(Jac - fdJ).max())
Phm, Tm = m.design(X3), m.one_hot(y3)
Wm = rng.normal(size=(3, 3)) * 0.3
gm, Hm = m.multiclass_grad_hess(Wm, Phm, Tm, 0.2)
f = lambda v: -m.multiclass_loglik(v.reshape(3, 3).T, Phm, Tm, 0.2)
gfd = fd_grad(f, Wm.T.ravel())
Hfd = np.array([fd_grad(lambda v, i=i: m.multiclass_grad_hess(v.reshape(3, 3).T, Phm, Tm, 0.2)[0][i], Wm.T.ravel()) for i in range(9)])
check("multiclass gradient (4.109) and Hessian (4.110, without the erratum's minus sign) match finite differences",
      np.abs(gm - gfd).max() < 1e-5 and np.abs(Hm - Hfd).max() < 1e-4, "grad err %.1e, Hessian err %.1e" % (np.abs(gm - gfd).max(), np.abs(Hm - Hfd).max()))
Wl = m.multiclass_irls(Phm, Tm, alpha=1e-4)
gl, _ = m.multiclass_grad_hess(Wl, Phm, Tm, 1e-4)
check("multiclass IRLS reaches a zero gradient and classifies the Figure 4.5 data better than least squares",
      np.abs(gl).max() < 1e-6 and ((Phm @ Wl).argmax(1) == y3).mean() > ((Phm @ m.least_squares_classifier(Phm, Tm)).argmax(1) == y3).mean() + 0.1,
      "accuracy %.3f (logistic) vs %.3f (least squares)" % (((Phm @ Wl).argmax(1) == y3).mean(), ((Phm @ m.least_squares_classifier(Phm, Tm)).argmax(1) == y3).mean()))

# probit
lam = math.sqrt(math.pi / 8)
check("probit: Phi(a) = (1 + erf(a / sqrt 2)) / 2 (4.115) and the slopes of sigma(a) and Phi(lambda a) at 0 agree for lambda^2 = pi/8 (Section 4.3.5)",
      np.allclose(ndtr(a), 0.5 * (1 + np.vectorize(math.erf)(a / math.sqrt(2)))) and abs(lam / math.sqrt(2 * math.pi) - 0.25) < 1e-15,
      "slope of Phi(lambda a) at 0 = %.12f" % (lam / math.sqrt(2 * math.pi)))
wp0 = rng.normal(size=3) * 0.5
gp = m.probit_grad(wp0, Phl, tl)
gfd = fd_grad(lambda w: -m.probit_loglik(w, Phl, tl), wp0)
check("probit gradient of the negative log likelihood matches finite differences", np.abs(gp - gfd).max() < 1e-6, "max err %.1e" % np.abs(gp - gfd).max())
wpr, ip = m.probit_irls(Phl, tl)
resp = optimize.minimize(lambda w: -m.probit_loglik(w, Phl, tl), np.zeros(3), jac=lambda w: m.probit_grad(w, Phl, tl), method="BFGS", options=dict(gtol=1e-10))
check("probit IRLS (Fisher scoring) increases the log likelihood and equals the BFGS optimum",
      (np.diff(ip["loglik"]) >= -1e-10).all() and np.abs(wpr - resp.x).max() < 1e-5, "%d iterations, max |diff| = %.1e; |w_probit| / |w_logistic| = %.3f (about 1/%.3f expected)"
      % (ip["iters"], np.abs(wpr - resp.x).max(), np.linalg.norm(wpr) / np.linalg.norm(w), 1 / lam))

# ---------------------------------------------------------------------------------------- 4.4
ln_f, dln_f, d2ln_f = m.fig414_functions(20.0, 4.0)
z0, Acurv, lnZl = m.laplace_1d(ln_f, dln_f, d2ln_f, 0.0)
eps = 1e-4
A_fd = -(ln_f(z0 + eps) - 2 * ln_f(z0) + ln_f(z0 - eps)) / eps ** 2
Zq, _ = integrate.quad(lambda z: math.exp(ln_f(z)), -10, 10, points=[z0], limit=300)
check("Laplace, p(z) ~ exp(-z^2/2) sigma(20 z + 4) (Figure 4.14): the mode, the curvature A (4.127) by finite differences, and ln Z (4.126) against quadrature",
      abs(dln_f(z0)) < 1e-10 and abs(A_fd - Acurv) < 1e-4 and abs(lnZl - math.log(Zq)) < 0.1,
      "z0 = %.4f, A = %.4f (finite difference %.4f), ln Z: Laplace %.4f, quadrature %.4f" % (z0, Acurv, A_fd, lnZl, math.log(Zq)))
lg = lambda z: -0.5 * z * z * 3.0
zz, AA, lz = m.laplace_1d(lg, lambda z: -3.0 * z, lambda z: -3.0, 1.0)
check("Laplace approximation is exact for a Gaussian: ln Z = 1/2 ln(2 pi / 3)", abs(lz - 0.5 * math.log(2 * math.pi / 3)) < 1e-12)

# evidence in 1-D logistic regression: quadrature versus Laplace (4.137), and BIC (4.139)
def one_d_evidence(N, alpha, seed):
    r = np.random.default_rng(seed)
    x = r.normal(size=N)
    t = (r.uniform(size=N) < expit(1.5 * x)).astype(float)
    Ph = x[:, None]
    ll = lambda w: m.logistic_loglik(np.array([w]), Ph, t)
    prior = lambda w: math.exp(0.5 * math.log(alpha / (2 * math.pi)) - 0.5 * alpha * w * w)
    wmap, _ = m.irls_logistic(Ph, t, alpha=alpha)
    Z, _ = integrate.quad(lambda w: math.exp(ll(w) - ll(wmap[0])) * prior(w), -15, 15, points=[wmap[0]], limit=400)
    ln_ex = math.log(Z) + ll(wmap[0])
    _, S, lap = m.laplace_logistic(Ph, t, alpha)
    return ln_ex, lap, Ph, t


errs = []
for N in (5, 20, 100, 500):
    ex, lap, Ph, t = one_d_evidence(N, 1.0, 11 + N)
    errs.append(abs(ex - lap))
check("Laplace evidence (4.137) against numerical integration of the evidence in 1-D (Bayesian logistic regression, alpha = 1): error below 0.15 and shrinking with N",
      max(errs) < 0.15 and errs[-1] < errs[0], "|ln p(D) quadrature - Laplace| for N = 5, 20, 100, 500: %s" % ", ".join("%.4f" % e for e in errs))
ex, lap, Ph, t = one_d_evidence(500, 1.0, 511)
bic = m.bic_logistic(Ph, t)
check("for large N the Laplace evidence is close to the BIC (4.139): they differ by O(1) while ln p(D) is O(N)", abs(lap - bic) < 3.0, "N = 500: Laplace %.3f, BIC %.3f" % (lap, bic))

# ---------------------------------------------------------------------------------------- 4.5
mu_a, s2 = 1.3, 2.0
exact, _ = integrate.quad(lambda a_: expit(a_) * math.exp(-0.5 * (a_ - mu_a) ** 2 / s2) / math.sqrt(2 * math.pi * s2), -30, 30, limit=300)
quad_probit, _ = integrate.quad(lambda a_: ndtr(lam * a_) * math.exp(-0.5 * (a_ - mu_a) ** 2 / s2) / math.sqrt(2 * math.pi * s2), -30, 30, limit=300)
closed = ndtr(lam * mu_a / math.sqrt(1 + lam ** 2 * s2))
check("convolution of Phi(lambda a) with N(a | mu, sigma^2) equals Phi(lambda mu / sqrt(1 + lambda^2 sigma^2)) (4.152)", abs(quad_probit - closed) < 1e-10, "%.12f vs %.12f" % (quad_probit, closed))
approx = float(expit(m.kappa(s2) * mu_a))
check("sigma-convolution: sigma(kappa(sigma^2) mu) (4.153)-(4.155) approximates the exact integral", abs(approx - exact) < 0.01, "exact %.5f, approximation %.5f" % (exact, approx))
mx = 0.0
for mu_ in np.linspace(-6, 6, 25):
    for s_ in (0.1, 1.0, 4.0, 16.0):
        e_, _ = integrate.quad(lambda a_: expit(a_) * math.exp(-0.5 * (a_ - mu_) ** 2 / s_) / math.sqrt(2 * math.pi * s_), -60, 60, limit=300)
        mx = max(mx, abs(e_ - float(expit(m.kappa(s_) * mu_))))
check("the approximation error stays below 0.02 over mu in [-6, 6] and sigma^2 in {0.1, 1, 4, 16}", mx < 0.02, "max error %.4f" % mx)
# predictive against Monte Carlo from the Laplace posterior
Xb, yb = m.blobs(m.JSRandom(11), 8, 0.6)
Pb = m.design(Xb)
wm, S, _ = m.laplace_logistic(Pb, yb.astype(float), 0.1)
check("Laplace posterior: S_N^-1 = alpha I + Phi^T R Phi (4.143, erratum S_N -> S_N^-1) and the mode has zero gradient of the log posterior",
      np.allclose(np.linalg.inv(S), 0.1 * np.eye(3) + Pb.T @ (Pb * (expit(Pb @ wm) * (1 - expit(Pb @ wm)))[:, None])) and np.abs(m.logistic_grad_hess(wm, Pb, yb.astype(float), 0.1)[0]).max() < 1e-9)
test = np.array([[1, 0.0, 0.0], [1, 3.0, -2.0], [1, -4.0, 4.0], [1, 6.0, 6.0]])
samp = np.random.default_rng(5).multivariate_normal(wm, S, size=400000)
mc = expit(samp @ test.T).mean(0)
pr = m.bayes_predictive(test, wm, S)
check("Bayesian predictive (4.155) agrees with Monte Carlo averages of sigma(w^T phi) under q(w) to within 0.02", np.abs(mc - pr).max() < 0.02,
      "max |MC - (4.155)| = %.4f; MAP probabilities %s versus predictive %s" % (np.abs(mc - pr).max(), np.round(expit(test @ wm), 3), np.round(pr, 3)))
check("the predictive probability lies closer to 1/2 than the MAP probability away from the data (contours splay out), and the p = 1/2 boundary is unchanged (kappa > 0 leaves the sign of mu_a)",
      (np.abs(pr - 0.5) <= np.abs(expit(test @ wm) - 0.5) + 1e-12).all() and ((pr > 0.5) == (test @ wm > 0)).all())

print("\n%d checks: %d PASS, %d FAIL" % (len(RESULTS), sum(RESULTS), len(RESULTS) - sum(RESULTS)))