"""Numerical checks of the derivations in Bishop (2006), Chapter 12 (continuous latent variables). Run `python3 checks_ch12.py`: each check prints PASS or FAIL with the numbers behind it, and the exit status is 1 if any check fails. The algorithms are those of ch12_continuous_latent.py. """ import sys import numpy as np from scipy.optimize import minimize import common import ch12_continuous_latent as c 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 rand_cov_data(rng, N, D, scales): """N points in D dimensions with standard deviations `scales` along random orthonormal axes.""" Q, _ = np.linalg.qr(rng.standard_normal((D, D))) return (rng.standard_normal((N, D)) * scales) @ Q.T + rng.standard_normal(D) def subspace_gap(A, B): """Largest difference between the orthogonal projections onto span(A) and span(B).""" Pa = A @ np.linalg.pinv(A) Pb = B @ np.linalg.pinv(B) return np.abs(Pa - Pb).max() def main(): rng = np.random.default_rng(2006) # ------------------------------------------------------------------ Section 12.1 X = rand_cov_data(rng, 200, 5, np.array([3.0, 2.0, 1.0, 0.5, 0.2])) xbar, lam, U = c.pca(X) S = np.cov(X.T, bias=True) total = np.trace(S) worst = 0.0 for _ in range(20): Q, _ = np.linalg.qr(rng.standard_normal((5, 2))) # a random 2-D subspace Xt = xbar + (X - xbar) @ Q @ Q.T J = ((X - Xt) ** 2).sum(1).mean() var = np.trace(Q.T @ S @ Q) worst = max(worst, abs(var + J - total)) check("projected variance + reconstruction error = total variance Tr S, random 2-D subspaces", worst < 1e-10, "max |var + J - Tr S| = %.1e" % worst) J_pca = [c.distortion(X, xbar, U, M) for M in range(6)] J_formula = [lam[M:].sum() for M in range(6)] check("distortion J of the PCA subspace = sum of the discarded eigenvalues (12.18)", np.allclose(J_pca, J_formula, atol=1e-10), "J(M=2) = %.6f" % J_pca[2]) Js = [] for _ in range(200): Q, _ = np.linalg.qr(rng.standard_normal((5, 2))) Xt = xbar + (X - xbar) @ Q @ Q.T Js.append(((X - Xt) ** 2).sum(1).mean()) check("no random 2-D subspace has a smaller error than the principal one (Exercise 12.2)", min(Js) >= J_pca[2] - 1e-12, "PCA %.4f, best of 200 random %.4f" % (J_pca[2], min(Js))) u = rng.standard_normal((1000, 5)) u /= np.linalg.norm(u, axis=1, keepdims=True) q = np.einsum("nd,de,ne->n", u, S, u) check("u^T S u <= lambda_1 for unit u, with equality at u_1 (12.5), (12.6)", q.max() <= lam[0] + 1e-12 and abs(U[:, 0] @ S @ U[:, 0] - lam[0]) < 1e-12, "max over 1000 random u: %.4f, lambda_1 = %.4f" % (q.max(), lam[0])) # PCA for N < D (Section 12.1.4) Xs = rng.standard_normal((15, 40)) @ rng.standard_normal((40, 40)) _, lam_big, U_big = c.pca(Xs) _, lam_small, U_small = c.pca_small_n(Xs, 10) norms = np.linalg.norm(U_small, axis=0) cos = np.abs(np.sum(U_small * U_big[:, :10], 0)) check("N < D trick (12.28)-(12.30): same eigenvalues, unit eigenvectors (Ex 12.3), same directions", np.allclose(lam_small, lam_big[:10]) and np.allclose(norms, 1) and np.allclose(cos, 1), "max eigenvalue gap %.1e, min |cos| %.12f" % (np.abs(lam_small - lam_big[:10]).max(), cos.min())) check("at least D - N + 1 eigenvalues of S are zero when N < D", np.sum(lam_big < 1e-9 * lam_big[0]) >= 40 - 15 + 1, "%d zero eigenvalues" % np.sum(lam_big < 1e-9 * lam_big[0])) # whitening (12.24), (12.25) for name, Xw in (("Old Faithful", common.old_faithful()), ("5-D test data", X)): Y = c.whiten(Xw) check("whitened %s: zero mean and identity covariance (12.25)" % name, np.allclose(Y.mean(0), 0, atol=1e-12) and np.allclose(np.cov(Y.T, bias=True), np.eye(Y.shape[1]), atol=1e-10)) # ------------------------------------------------------------------ Section 12.2.1: PPCA X = rand_cov_data(rng, 300, 6, np.array([2.0, 1.5, 0.4, 0.3, 0.3, 0.2])) N, D, M = 300, 6, 2 mu, W, s2 = c.ppca_ml(X, M) _, lam, U = c.pca(X) check("sigma^2_ML = average of the discarded eigenvalues (12.46)", abs(s2 - lam[M:].mean()) < 1e-14) # Woodbury (12.40) and the determinant lemma C = W @ W.T + s2 * np.eye(D) Mm = W.T @ W + s2 * np.eye(M) Cinv_wb = (np.eye(D) - W @ np.linalg.solve(Mm, W.T)) / s2 check("C^{-1} by the Woodbury identity (12.40), (12.41) equals the direct inverse", np.allclose(Cinv_wb, np.linalg.inv(C), atol=1e-12), "max diff %.1e" % np.abs(Cinv_wb - np.linalg.inv(C)).max()) ld = (D - M) * np.log(s2) + np.linalg.slogdet(Mm)[1] check("ln|C| = (D - M) ln sigma^2 + ln|M| (matrix determinant lemma)", abs(ld - np.linalg.slogdet(C)[1]) < 1e-12) ll_direct = common.gauss_logpdf(X, mu, C).sum() check("ppca_loglik (12.43) via (12.40) equals the direct Gaussian log density", abs(c.ppca_loglik(X, mu, W, s2) - ll_direct) < 1e-8, "%.6f" % ll_direct) S = np.cov(X.T, bias=True) ll_trace = -N / 2 * (D * np.log(2 * np.pi) + np.linalg.slogdet(C)[1] + np.trace(np.linalg.solve(C, S))) check("log likelihood at mu = xbar equals the trace form (12.44)", abs(ll_trace - ll_direct) < 1e-8) # the ML solution beats random W, rotated W gives the same value (12.39) ll_ml = c.ppca_loglik(X, mu, W, s2) others = [c.ppca_loglik(X, mu, rng.standard_normal((D, M)), s2) for _ in range(50)] others += [c.ppca_loglik(X, mu, W * f, s2 * g) for f in (0.9, 1.1) for g in (0.8, 1.25)] check("log likelihood at W_ML (12.45), sigma^2_ML (12.46) exceeds random W and perturbed values", ll_ml > max(others), "%.4f > %.4f" % (ll_ml, max(others))) R, _ = np.linalg.qr(rng.standard_normal((M, M))) check("W R gives the same likelihood for any orthogonal R (12.39)", abs(c.ppca_loglik(X, mu, W @ R, s2) - ll_ml) < 1e-9) # stationary point: numerical gradient of ln p with respect to W and sigma^2 vanishes def f(p): return c.ppca_loglik(X, mu, p[:-1].reshape(D, M), p[-1]) p0 = np.r_[W.ravel(), s2] h = 1e-6 grad = np.array([(f(p0 + h * e) - f(p0 - h * e)) / (2 * h) for e in np.eye(len(p0))]) check("(12.45), (12.46) is a stationary point: finite-difference gradient ~ 0", np.abs(grad).max() < 1e-4, "max |grad| = %.1e" % np.abs(grad).max()) # other eigenvector subsets give stationary points (saddles) with lower likelihood (Tipping and # Bishop 1999); on data whose eigenvalues are close enough for every subset to be valid Xe = rand_cov_data(rng, 300, 6, np.sqrt([3.0, 2.0, 1.5, 1.2, 1.0, 0.9])) mue, We, s2e = c.ppca_ml(Xe, 2) _, lame, Ue = c.pca(Xe) lls, grads = [], [] for sub in ([0, 2], [1, 2], [0, 3], [2, 3]): rest = [i for i in range(6) if i not in sub] s2s = lame[rest].mean() if np.all(lame[sub] > s2s): Ws = Ue[:, sub] * np.sqrt(lame[sub] - s2s) lls.append(c.ppca_loglik(Xe, mue, Ws, s2s)) fe = lambda p: c.ppca_loglik(Xe, mue, p[:-1].reshape(6, 2), p[-1]) pe = np.r_[Ws.ravel(), s2s] grads.append(max(abs((fe(pe + 1e-6 * e) - fe(pe - 1e-6 * e)) / 2e-6) for e in np.eye(len(pe)))) check("other eigenvector subsets are stationary points with a lower likelihood (saddles)", len(lls) >= 2 and max(lls) < c.ppca_loglik(Xe, mue, We, s2e) and max(grads) < 1e-4, "%d subsets: best %.3f < ML %.3f, max |grad| %.1e" % (len(lls), max(lls), c.ppca_loglik(Xe, mue, We, s2e), max(grads))) # variance along retained and discarded directions check("model variance u_i^T C u_i = lambda_i along retained eigenvectors, sigma^2 orthogonal to them", np.allclose([U[:, i] @ C @ U[:, i] for i in range(M)], lam[:M]) and np.allclose([U[:, i] @ C @ U[:, i] for i in range(M, D)], s2)) # M = D gives C = S (12.47) muD, WD, s2D = c.ppca_ml(X, D - 1) check("with M = D - 1 the ML covariance equals the sample covariance (12.47, Ex 12.14)", np.allclose(WD @ WD.T + s2D * np.eye(D), S, atol=1e-12)) # posterior (12.42) against the joint Gaussian of (z, x) x0 = X[:5] Ez, Sz = c.ppca_posterior(x0, mu, W, s2) Cinv = np.linalg.inv(C) Ez_joint = (W.T @ Cinv @ (x0 - mu).T).T # E[z|x] = W^T C^{-1} (x - mu) Sz_joint = np.eye(M) - W.T @ Cinv @ W # cov = I - W^T C^{-1} W check("posterior (12.42) equals the conditional of the joint Gaussian of (z, x)", np.allclose(Ez, Ez_joint, atol=1e-12) and np.allclose(Sz, Sz_joint, atol=1e-12), "max diff %.1e" % max(np.abs(Ez - Ez_joint).max(), np.abs(Sz - Sz_joint).max())) # sigma^2 -> 0: posterior mean is the orthogonal projection (12.50) Ez0, _ = c.ppca_posterior(x0, mu, W, 1e-12) proj = np.linalg.solve(W.T @ W, W.T @ (x0 - mu).T).T check("as sigma^2 -> 0 the posterior mean becomes (W^T W)^{-1} W^T (x - xbar) (12.50)", np.allclose(Ez0, proj, atol=1e-8)) # shrinkage towards the origin for sigma^2 > 0 (Exercise 12.12) check("for sigma^2 > 0 the latent projection is shrunk towards the origin (Ex 12.12)", np.all(np.linalg.norm(Ez @ np.linalg.cholesky(W.T @ W), axis=1) < np.linalg.norm(proj @ np.linalg.cholesky(W.T @ W), axis=1))) # ------------------------------------------------------------------ Section 12.2.2: EM _, Wem, s2em, ll = c.ppca_em(X, M, n_iter=3000, sigma2=1.0, rng=np.random.default_rng(1)) check("PPCA EM (12.54)-(12.57): the log likelihood never decreases", np.all(np.diff(ll) > -1e-9), "%.3f -> %.6f over %d steps" % (ll[0], ll[-1], len(ll) - 1)) check("PPCA EM converges to the closed form: same likelihood, sigma^2 and W W^T (W up to R)", abs(ll[-1] - ll_ml) < 1e-6 and abs(s2em - s2) < 1e-8 and np.allclose(Wem @ Wem.T, W @ W.T, atol=1e-6), "ll %.8f vs %.8f; sigma^2 %.8f vs %.8f" % (ll[-1], ll_ml, s2em, s2)) Rhat = np.linalg.lstsq(W, Wem, rcond=None)[0] check("the EM solution is W_ML R with R orthogonal (rotational non-identifiability)", np.allclose(Rhat @ Rhat.T, np.eye(M), atol=1e-6)) # the M step (12.56), (12.57) maximises the expected complete-data log likelihood (12.53) Xc = X - X.mean(0) W1, s21 = rng.standard_normal((D, M)), 0.7 Minv = np.linalg.inv(W1.T @ W1 + s21 * np.eye(M)) Ezz_n = Xc @ W1 @ Minv sumEzz = N * s21 * Minv + Ezz_n.T @ Ezz_n def Q(Wq, s2q): # (12.53) with mu = xbar, including the M/2 ln(2 pi) term of the errata return -(N * D / 2 * np.log(2 * np.pi * s2q) + 0.5 * np.trace(sumEzz) + ((Xc ** 2).sum() - 2 * np.sum(Ezz_n * (Xc @ Wq)) + np.trace(sumEzz @ Wq.T @ Wq)) / (2 * s2q) + N * M / 2 * np.log(2 * np.pi)) Wn, s2n = c.ppca_em_step(Xc, W1, s21) pert = max(Q(Wn + 1e-3 * rng.standard_normal((D, M)), s2n * np.exp(1e-3 * rng.standard_normal())) for _ in range(200)) check("the M step maximises the expected complete-data log likelihood (12.53)", Q(Wn, s2n) >= pert, "Q(new) = %.6f >= best perturbation %.6f" % (Q(Wn, s2n), pert)) # sigma^2 -> 0 EM (12.58), (12.59) finds the principal subspace Wpca = c.em_pca(X, rng.standard_normal((D, M)), n_iter=500) check("EM for PCA in the limit sigma^2 -> 0 (12.58), (12.59) converges to the principal subspace", subspace_gap(Wpca, U[:, :M]) < 1e-8, "projection gap %.1e" % subspace_gap(Wpca, U[:, :M])) Xc10 = c.EM_PCA_POINTS - c.EM_PCA_POINTS.mean(0) _, Wn = c.em_pca_step(Xc10, np.array([[2.0], [-1.0]])) Z_old = Xc10 @ np.array([2.0, -1.0]) / 5.0 resid = Xc10 - np.outer(Z_old, Wn[:, 0]) check("after an M step of EM for PCA the residuals are orthogonal to the old W (Fig 12.12c)", np.abs(resid @ np.array([2.0, -1.0])).max() < 1e-12) # missing data: exact EM on a small problem against direct maximisation Xm = rand_cov_data(np.random.default_rng(7), 40, 5, np.array([2.0, 1.0, 0.3, 0.3, 0.3])) obs = c.missing_mask(40, 5, 0.3, np.random.default_rng(8)) mu_m, W_m, s2_m, llm = c.ppca_em_missing(Xm, obs, 2, n_iter=5000, rng=np.random.default_rng(9), tol=1e-14) check("EM with missing data: the observed-data log likelihood never decreases", np.all(np.diff(llm) > -1e-9), "%.3f -> %.6f in %d steps" % (llm[0], llm[-1], len(llm) - 1)) def nll(p): return -c.ppca_loglik_missing(Xm, obs, p[:5], p[5:15].reshape(5, 2), np.exp(p[15])) best = -np.inf for seed in range(3): r = np.random.default_rng(100 + seed) p0 = np.r_[np.nanmean(np.where(obs, Xm, np.nan), 0), r.standard_normal(10), 0.0] res = minimize(nll, p0, method="BFGS", options=dict(maxiter=5000, gtol=1e-8)) best = max(best, -res.fun) check("EM with missing data reaches the maximum found by direct quasi-Newton optimisation", abs(llm[-1] - best) < 1e-5, "EM %.7f, BFGS %.7f" % (llm[-1], best)) mu_c, W_c, s2_c, llc = c.ppca_em_missing(X, np.ones_like(X, bool), M, n_iter=3000, rng=np.random.default_rng(3), tol=1e-13) check("with nothing missing, the missing-data EM gives the closed-form ML solution", abs(llc[-1] - ll_ml) < 1e-6 and np.allclose(mu_c, mu, atol=1e-6)) # ------------------------------------------------------------------ Section 12.2.3: Bayesian PCA Xb = c.bpca_data(rng=np.random.default_rng(1214)) _, Wb, s2b, alpha, hist = c.bayesian_pca(Xb, 9, n_iter=2000, W=np.random.default_rng(5).standard_normal((10, 9))) kept = np.sum(alpha < 1e6) check("Bayesian PCA (12.62), (12.63) on the Fig 12.14 data keeps exactly 3 of 9 columns", kept == 3, "alpha of kept columns: %s" % np.array2string(np.sort(alpha)[:kept], precision=2)) check("at convergence alpha_i = D / ||w_i||^2 for the kept columns (12.62)", np.allclose(alpha[alpha < 1e6], 10 / (Wb[:, alpha < 1e6] ** 2).sum(0), rtol=1e-6)) Xc = Xb - Xb.mean(0) Minv = np.linalg.inv(Wb.T @ Wb + s2b * np.eye(9)) Ez = Xc @ Wb @ Minv sumEzz = len(Xb) * s2b * Minv + Ez.T @ Ez Wfix = (Xc.T @ Ez) @ np.linalg.inv(sumEzz + s2b * np.diag(alpha)) # ARD breaks the rotational symmetry among the kept columns only weakly, so EM creeps slowly check("the Bayesian PCA solution is (numerically) a fixed point of the modified M step (12.63)", np.abs(Wfix - Wb).max() < 1e-4 * np.abs(Wb).max(), "max change %.1e" % np.abs(Wfix - Wb).max()) # ------------------------------------------------------------------ Section 12.2.4: factor analysis Xf = rand_cov_data(rng, 400, 6, np.array([2.0, 1.2, 0.5, 0.4, 0.3, 0.2])) Xf = Xf * np.array([1.0, 3.0, 0.5, 1.0, 2.0, 1.0]) W0f = np.random.default_rng(11).standard_normal((6, 2)) psi0 = np.diag(np.cov(Xf.T, bias=True)).copy() muf, Wf, psi, llf = c.factor_analysis_em(Xf, 2, n_iter=3000, W=W0f, psi=psi0) check("factor analysis EM (12.66)-(12.70): the log likelihood never decreases", np.all(np.diff(llf) > -1e-9), "%.3f -> %.6f" % (llf[0], llf[-1])) G = np.linalg.inv(np.eye(2) + Wf.T @ (Wf / psi[:, None])) Cf = Wf @ Wf.T + np.diag(psi) check("E step (12.66) with G (12.68) equals the joint-Gaussian posterior W^T C^{-1} (x - mu)", np.allclose(G @ Wf.T @ np.diag(1 / psi), Wf.T @ np.linalg.inv(Cf), atol=1e-10)) # equivariance under rescaling (Section 12.2.4, Exercise 12.25) Dsc = np.diag([10.0, 1, 1, 1, 1, 1]) # started from the correspondingly rescaled point, every EM iterate is exactly the rescaled one _, Wf2, psi2, llf2 = c.factor_analysis_em(Xf @ Dsc, 2, n_iter=3000, W=Dsc @ W0f, psi=psi0 * np.diag(Dsc) ** 2) Cf2 = Wf2 @ Wf2.T + np.diag(psi2) check("factor analysis: rescaling a variable rescales Psi and the row of W (C -> D C D, Ex 12.25)", np.allclose(Wf2, Dsc @ Wf, atol=1e-9) and np.allclose(psi2, psi * np.diag(Dsc) ** 2, rtol=1e-9) and np.allclose(llf2, llf - len(Xf) * np.log(10.0), atol=1e-6), "max rel. diff of C %.1e" % (np.abs(Cf2 - Dsc @ Cf @ Dsc).max() / np.abs(Cf2).max())) _, Wp1, s2p1 = c.ppca_ml(Xf, 2) _, Wp2, s2p2 = c.ppca_ml(Xf @ Dsc, 2) Cp1 = Wp1 @ Wp1.T + s2p1 * np.eye(6) Cp2 = Wp2 @ Wp2.T + s2p2 * np.eye(6) rel_fa = np.abs(Cf2 - Dsc @ Cf @ Dsc).max() / np.abs(Cf2).max() rel_pp = np.abs(Cp2 - Dsc @ Cp1 @ Dsc).max() / np.abs(Cp2).max() check("PPCA is not equivariant under rescaling: C changes by much more than D C D", rel_pp > 1e-2 > 1e3 * rel_fa, "relative change: PPCA %.1e, factor analysis %.1e" % (rel_pp, rel_fa)) R2, _ = np.linalg.qr(rng.standard_normal((6, 6))) _, Wp3, s2p3 = c.ppca_ml(Xf @ R2.T, 2) check("PPCA is equivariant under rotations of data space: W -> R W (Section 12.2.4)", np.allclose(Wp3 @ Wp3.T, R2 @ Wp1 @ Wp1.T @ R2.T, atol=1e-9) and abs(s2p3 - s2p1) < 1e-12) # ------------------------------------------------------------------ Section 12.3: kernel PCA Xk = rng.standard_normal((30, 2)) * [1.0, 0.5] + [0.3, -0.2] def phi(A): # explicit feature map of k(x, x') = (1 + x^T x')^2 in 2-D x1, x2 = A[:, 0], A[:, 1] r2 = np.sqrt(2.0) return np.column_stack([np.ones(len(A)), r2 * x1, r2 * x2, x1 ** 2, x2 ** 2, r2 * x1 * x2]) def kpoly(A, B): return (1.0 + A @ B.T) ** 2 F = phi(Xk) check("explicit feature map reproduces the polynomial kernel", np.allclose(F @ F.T, kpoly(Xk, Xk))) Fc = F - F.mean(0) Kt = c.centre_gram(kpoly(Xk, Xk)) check("centred Gram matrix (12.85) equals the Gram matrix of explicitly centred features (12.83)", np.allclose(Kt, Fc @ Fc.T, atol=1e-10), "max diff %.1e" % np.abs(Kt - Fc @ Fc.T).max()) model = c.kernel_pca(Xk, kpoly, 4) lam_f = np.sort(np.linalg.eigvalsh(Fc.T @ Fc / len(Xk)))[::-1][:4] check("eigenvalues of K~ are N lambda_i, lambda_i the feature-space covariance eigenvalues (12.80)", np.allclose(model["eig"] / len(Xk), lam_f, atol=1e-10)) A = model["A"] check("normalisation (12.81): lambda_i N a_i^T a_i = 1 so that v_i^T v_i = 1", np.allclose(model["eig"] * (A ** 2).sum(0), 1) and np.allclose(np.linalg.norm(Fc.T @ A, axis=0), 1)) Xnew = rng.standard_normal((7, 2)) yk = c.kernel_pca_project(model, Xnew) _, _, Uf = np.linalg.svd(Fc, full_matrices=False) yf = (phi(Xnew) - F.mean(0)) @ Uf[:4].T check("kernel PCA projections (12.82) equal linear PCA projections in feature space (up to sign)", np.allclose(np.abs(yk), np.abs(yf), atol=1e-8), "max diff %.1e" % np.abs(np.abs(yk) - np.abs(yf)).max()) yk_nc = c.kernel_pca_project(dict(model, centre=False), Xnew) shift = yk_nc - yk check("without centring the test kernel, (12.87) differs only by a constant per component", np.allclose(shift, shift[0], atol=1e-10)) model_lin = c.kernel_pca(Xk, lambda A, B: A @ B.T, 2) _, lam_lin, U_lin = c.pca(Xk) check("a linear kernel recovers standard PCA (Exercise 12.27)", np.allclose(model_lin["eig"] / len(Xk), lam_lin) and np.allclose(np.abs(c.kernel_pca_project(model_lin, Xnew)), np.abs((Xnew - Xk.mean(0)) @ U_lin), atol=1e-10)) Xt = c.kpca_toy_data() mt = c.kernel_pca(Xt, c.gaussian_kernel, 8) check("Fig 12.17 data: two leading eigenvalues near 21.7 (book: 21.72, 21.65), then a gap", 18 < mt["eig"][1] < mt["eig"][0] < 25 and mt["eig"][2] < 6, "eigenvalues %s" % np.array2string(mt["eig"], precision=2)) # ------------------------------------------------------------------ Section 12.4.3: GTM Xo, _ = common.oil_flow("trn") alpha_gtm = 0.01 g = c.gtm_fit(Xo, n_iter=30, alpha=alpha_gtm) # the monotone quantity of the penalised EM is ln p(X | W, beta) + ln p(W | alpha) pen = [] for it in range(31): gi = c.gtm_fit(Xo, n_iter=it, alpha=alpha_gtm) if it in (0, 1, 2, 5, 10, 20, 30) else None if gi is not None: pen.append(gi["loglik"][-1] - 0.5 * alpha_gtm * np.sum(gi["W"] ** 2)) check("GTM EM: the penalised log likelihood never decreases (iterations 0, 1, 2, 5, 10, 20, 30)", np.all(np.diff(pen) > -1e-6), "%.1f -> %.1f" % (pen[0], pen[-1])) check("GTM responsibilities sum to one over the latent grid", np.allclose(g["R"].sum(0), 1)) # ------------------------------------------------------------------ book numbers (Figures 12.3, 12.4) _, lam3, _ = c.threes_pca() check("MNIST test threes, pixels 0-255: lambda_1..4 = 3.4, 2.8, 2.4, 1.6 x 10^5 as in Fig 12.3", np.allclose(np.round(lam3[:4] / 1e5, 1), [3.4, 2.8, 2.4, 1.6]), "%s" % np.array2string(lam3[:4] / 1e5, precision=3)) print("\n%d checks: %d PASS, %d FAIL" % (len(RESULTS), sum(RESULTS), len(RESULTS) - sum(RESULTS))) return 0 if all(RESULTS) else 1 if __name__ == "__main__": sys.exit(main())