{ "cells": [ { "cell_type": "markdown", "id": "60cc263e", "metadata": {}, "source": [ "# CS6140 Machine Learning — Fall 2026\n", "# HW4 Starter — Generative Methods\n", "\n", "This notebook implements `HW4_26F.html`. For Problems 1, 3, and 5 you'll first solve it with a\n", "well-known library (a sanity-check baseline), then implement the same thing from scratch and\n", "evaluate it exactly the same way, so your two numbers should land right next to each other once\n", "everything is filled in. Problem 2 has no single library call that reproduces all three likelihood\n", "variants, so it's validated by scoring all three the same way instead; Problem 4 has no standard\n", "library counterpart for a Binomial mixture, so it's validated by recovering known, self-generated\n", "ground-truth parameters (same idea as Problem 3).\n", "\n", "**Problems 1-4 (required):**\n", "1. Gaussian Discriminant Analysis (GDA): LDA vs. QDA, k-fold CV — **Spambase** [40 points]\n", "2. Naive Bayes: Bernoulli / Gaussian / non-parametric (histogram), 10-fold CV, ROC/threshold\n", " discussion — **Spambase** [50 points]\n", "3. EM for a Gaussian mixture, own E/M steps, recovered vs. true parameters — **2gaussian.txt**,\n", " **3gaussian.txt** [50 points]\n", "4. EM for a mixture of Binomials, own E/M steps, self-generated coin-flip data — **generated\n", " data** [50 points]\n", "\n", "**Problem 5 (required, Parts A/B only; Part C optional/no-credit):**\n", "5. Bayesian linear regression: Ridge as the MAP estimate under a Gaussian weight prior, and the\n", " full posterior/predictive distribution — **Housing** (bridges back to HW1's closed-form\n", " Ridge) [60 points]. Part C (choosing λ via the evidence, and a worked-through cautionary tale\n", " about comparing that to cross-validation) is optional, no credit -- see\n", " `../lecture_notes/bayesian_ridge_regression.pdf` if you want the full depth on it.\n", "\n", "**How to work through this notebook.** Everything runs except the lines marked\n", "`### TODO_STUDENT ###` — those have their code replaced with a `...` placeholder and a one-line\n", "description of what to implement there. Read the surrounding (working) code and the markdown\n", "above each section for context, then replace each `...` with real code. A cell with a `...` in it\n", "will raise a `SyntaxError` or produce obviously wrong output until you fill it in — that's\n", "expected, not a bug. Work top to bottom: later cells depend on classes/functions built in earlier\n", "ones (e.g. Problem 3's `gaussian_mixture_em` and Problem 4's `binomial_mixture_em` are independent\n", "of each other, but Problem 5's Part C -- if you attempt it -- reuses Part B's\n", "`bayesian_ridge_posterior`). Once a \"library baseline\" cell's TODOs are filled in and run, the\n", "printed number is a target: your matching \"from scratch\" implementation just below it should land\n", "on (or very near) the same value.\n", "\n", "Datasets are read from the shared course `data/` folder (paths relative to\n", "`3_generative_models/HW3/`, i.e. two levels up, same depth as the solution's `code_HW3/`) and\n", "`2gaussian.txt`/`3gaussian.txt`, which live right here next to this notebook.\n", "\n", "**Requirements.**\n", "1. Fill in every `TODO_STUDENT` blank according to the algorithm steps covered in lecture and the\n", " linked materials, and make sure the notebook runs top to bottom without errors.\n", "2. Understand the *whole* notebook — including the provided/given code, not just the blanks you\n", " filled in — well enough to explain any part of it during office hours.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "074640ac", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import multivariate_normal\n", "from sklearn.model_selection import StratifiedKFold, KFold, train_test_split\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.metrics import accuracy_score, roc_curve, auc\n", "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis\n", "from sklearn.naive_bayes import GaussianNB, BernoulliNB\n", "from sklearn.mixture import GaussianMixture\n", "from sklearn.linear_model import Ridge, BayesianRidge\n", "\n", "DATA = \"../../data\"\n", "np.random.seed(42)\n" ] }, { "cell_type": "markdown", "id": "ce8617bc", "metadata": {}, "source": [ "## Data loading" ] }, { "cell_type": "code", "execution_count": null, "id": "585e7007", "metadata": {}, "outputs": [], "source": [ "# ---- Spambase (binary classification): one file, used raw (unscaled) throughout this HW.\n", "# GDA/Naive Bayes are per-feature or per-class density estimators, so unlike SVM/kernel methods\n", "# they don't need feature scaling for correctness -- only Problem 1's covariance inversion\n", "# benefits from it, which is handled locally there.\n", "spam_raw = np.loadtxt(f\"{DATA}/spambase/spambase.data\", delimiter=\",\")\n", "X_spam, y_spam = spam_raw[:, :-1], spam_raw[:, -1].astype(int) # y in {0,1}, 1 = spam\n", "\n", "print(\"Spambase:\", X_spam.shape, \" spam fraction:\", y_spam.mean().round(3))\n" ] }, { "cell_type": "markdown", "id": "0b33a6d4", "metadata": {}, "source": [ "## PROBLEM 1 — Gaussian Discriminant Analysis (GDA)   [40 points]\n", "\n", "Fit class-conditional Gaussians to Spambase under two covariance assumptions -- a single shared\n", "covariance (**LDA**) and a separate covariance per class (**QDA**) -- and compare them with\n", "k-fold cross validation.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4bf9b223", "metadata": {}, "outputs": [], "source": [ "def gaussian_log_density(X, mean, cov):\n", " \"\"\"log N(x; mean, cov) for every row of X, via slogdet (stable even for Spambase's 57x57\n", " covariance, where a raw determinant would under/overflow long before it's ever logged).\"\"\"\n", " d = mean.shape[0]\n", " sign, logdet = np.linalg.slogdet(cov)\n", " inv_cov = np.linalg.inv(cov)\n", " diff = X - mean\n", " quad = np.einsum... ### TODO_STUDENT ### the Mahalanobis quadratic form (x-mean)^T @ inv_cov @ (x-mean), for every row at once\n", " return -0.5 * (d * np.log(2 * np.pi) + logdet + quad)\n", "\n", "\n", "def gda_fit(X, y, shared_cov):\n", " \"\"\"MLE class means + covariance(s). shared_cov=True pools into one covariance (LDA);\n", " shared_cov=False estimates one covariance per class (QDA). A small ridge is added to each\n", " covariance purely for invertibility -- Spambase's features are highly correlated.\"\"\"\n", " classes = np.unique(y)\n", " d = X.shape[1]\n", " means, covs, priors, ridge = {}, {}, {}, 1e-6 * np.eye(d)\n", " for c in classes:\n", " Xc = X[y == c]\n", " priors[c] = len(Xc) / len(X)\n", " means[c] = ... ### TODO_STUDENT ### class-conditional mean (MLE)\n", " centered = Xc - means[c]\n", " covs[c] = ... ### TODO_STUDENT ### class-conditional covariance (MLE, biased by n not n-1)\n", " if shared_cov:\n", " pooled = ... ### TODO_STUDENT ### pool the per-class covariances, weighted by class size, into one shared covariance\n", " covs = {c: pooled for c in classes}\n", " return {\"classes\": classes, \"means\": means, \"covs\": {c: covs[c] + ridge for c in classes},\n", " \"priors\": priors}\n", "\n", "\n", "def gda_predict(X, model):\n", " log_post = np.stack(\n", " [np.log(model[\"priors\"][c]) + gaussian_log_density(X, model[\"means\"][c], model[\"covs\"][c])\n", " for c in model[\"classes\"]], axis=1)\n", " return model[\"classes\"][np.argmax(log_post, axis=1)]\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c4edd4b3", "metadata": {}, "outputs": [], "source": [ "def k_fold_gda(X, y, shared_cov, k=10, seed=42):\n", " skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=seed)\n", " train_accs, test_accs = [], []\n", " for train_idx, test_idx in skf.split(X, y):\n", " model = gda_fit(X[train_idx], y[train_idx], shared_cov)\n", " train_accs.append(accuracy_score(y[train_idx], gda_predict(X[train_idx], model)))\n", " test_accs.append(accuracy_score(y[test_idx], gda_predict(X[test_idx], model)))\n", " return np.mean(train_accs), np.mean(test_accs)\n", "\n", "\n", "def k_fold_library(X, y, make_model, k=10, seed=42):\n", " \"\"\"Same folds, sklearn's LDA/QDA as a library baseline.\"\"\"\n", " skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=seed)\n", " train_accs, test_accs = [], []\n", " for train_idx, test_idx in skf.split(X, y):\n", " clf = make_model... ### TODO_STUDENT ### fit sklearn's LDA/QDA as the library baseline for this fold\n", " train_accs.append(accuracy_score(y[train_idx], clf.predict(X[train_idx])))\n", " test_accs.append(accuracy_score(y[test_idx], clf.predict(X[test_idx])))\n", " return np.mean(train_accs), np.mean(test_accs)\n", "\n", "\n", "results = []\n", "for name, shared in [(\"LDA (shared cov)\", True), (\"QDA (separate cov)\", False)]:\n", " tr, te = k_fold_gda(X_spam, y_spam, shared_cov=shared, k=10)\n", " results.append({\"model\": f\"{name} — scratch\", \"train_acc\": tr, \"test_acc\": te})\n", "\n", "# QDA needs its own small ridge (reg_param) here too -- Spambase's per-class covariances are\n", "# rank-deficient enough that sklearn >=1.9 raises a hard LinAlgError with the unregularized\n", "# default (older sklearn only warned and limped along); reg_param=1e-3 matches our own scratch\n", "# ridge in spirit and keeps this baseline runnable across sklearn versions\n", "for name, make_model in [(\"LDA (shared cov)\", lambda: LinearDiscriminantAnalysis()),\n", " (\"QDA (separate cov)\", lambda: QuadraticDiscriminantAnalysis(reg_param=1e-3))]:\n", " tr, te = k_fold_library(X_spam, y_spam, make_model, k=10)\n", " results.append({\"model\": f\"{name} — sklearn\", \"train_acc\": tr, \"test_acc\": te})\n", "\n", "gda_results = pd.DataFrame(results).set_index(\"model\").round(4)\n", "gda_results\n" ] }, { "cell_type": "markdown", "id": "40975777", "metadata": {}, "source": [ "**Does the data look normally distributed?** Compare LDA's and QDA's train/test accuracy in the\n", "table above and discuss.\n" ] }, { "cell_type": "markdown", "id": "4f40503e", "metadata": {}, "source": [ "## PROBLEM 2 — Naive Bayes   [50 points]\n", "\n", "Three Naive Bayes variants — Bernoulli, Gaussian, and a non-parametric histogram — sharing one\n", "scoring convention: every `*_fit` function returns a small model dict/object, and every model can\n", "produce a **log-odds score** `log P(spam|x) - log P(non-spam|x)`, which both classifies (sign) and\n", "feeds the ROC/threshold analysis below.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f9d01f06", "metadata": {}, "outputs": [], "source": [ "def bernoulli_nb_fit(X, y, alpha=1.0):\n", " mu = ... ### TODO_STUDENT ### per-feature Bernoulli threshold: the training-set mean of each feature\n", " X_bin = (X > mu).astype(int)\n", " p1, priors = {}, {}\n", " for c in (0, 1):\n", " Xc = X_bin[y == c]\n", " priors[c] = len(Xc) / len(X)\n", " p1[c] = ... ### TODO_STUDENT ### additive-smoothed P(feature above threshold | class)\n", " return {\"mu\": mu, \"priors\": priors, \"p1\": p1}\n", "\n", "\n", "def bernoulli_nb_log_odds(X, model):\n", " X_bin = (X > model[\"mu\"]).astype(int)\n", " def class_loglik(c):\n", " p1 = model[\"p1\"][c]\n", " return ... ### TODO_STUDENT ### sum of per-feature Bernoulli log-likelihoods, plus the log prior\n", " return class_loglik(1) - class_loglik(0)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "41c1e25d", "metadata": {}, "outputs": [], "source": [ "def gaussian_nb_fit(X, y, var_smoothing=1e-9):\n", " eps = var_smoothing * X.var(axis=0).max() # same variance floor convention as sklearn's GaussianNB\n", " mean, var, priors = {}, {}, {}\n", " for c in (0, 1):\n", " Xc = X[y == c]\n", " priors[c] = len(Xc) / len(X)\n", " mean[c] = ... ### TODO_STUDENT ### class-conditional mean per feature (MLE)\n", " var[c] = ... ### TODO_STUDENT ### class-conditional variance per feature (MLE, plus a numerical-stability floor)\n", " return {\"priors\": priors, \"mean\": mean, \"var\": var}\n", "\n", "\n", "def gaussian_nb_log_odds(X, model):\n", " def class_loglik(c):\n", " mean, var = model[\"mean\"][c], model[\"var\"][c]\n", " log_p = ... ### TODO_STUDENT ### per-feature Gaussian log-density\n", " return np.log(model[\"priors\"][c]) + log_p.sum(axis=1)\n", " return class_loglik(1) - class_loglik(0)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5018c954", "metadata": {}, "outputs": [], "source": [ "class HistogramNaiveBayes:\n", " \"\"\"Non-parametric per-feature likelihood: a small histogram whose bin edges are placed at\n", " the feature's overall range and its two class-conditional means, so the bins can actually\n", " separate the classes even with very few bins.\"\"\"\n", "\n", " def __init__(self, alpha=1.0):\n", " self.alpha = alpha\n", "\n", " def fit(self, X, y):\n", " self.priors = {c: (np.sum(y == c) + self.alpha) / (len(y) + 2 * self.alpha) for c in (0, 1)}\n", " self.edges, self.bin_probs = [], {0: [], 1: []}\n", " for j in range(X.shape[1]):\n", " col = X[:, j]\n", " mean0, mean1 = col[y == 0].mean(), col[y == 1].mean()\n", " lo, hi = sorted((mean0, mean1))\n", " edges = ... ### TODO_STUDENT ### data-driven bin edges: feature range plus its two class-conditional means\n", " if len(edges) < 2:\n", " edges = [col.min(), col.min() + 1e-9]\n", " self.edges.append(edges)\n", " for c in (0, 1):\n", " counts, _ = np.histogram(col[y == c], bins=edges)\n", " self.bin_probs[c].append(...) ### TODO_STUDENT ### additive-smoothed P(feature in bin | class)\n", " return self\n", "\n", " def _class_loglik(self, X, c):\n", " log_lik = np.full(X.shape[0], np.log(self.priors[c]))\n", " for j in range(X.shape[1]):\n", " bin_idx = np.searchsorted(self.edges[j], X[:, j]) - 1\n", " bin_idx = np.clip(bin_idx, 0, len(self.bin_probs[c][j]) - 1)\n", " log_lik += ... ### TODO_STUDENT ### look up and accumulate each feature's log bin-probability\n", " return log_lik\n", "\n", " def log_odds(self, X):\n", " return self._class_loglik(X, 1) - self._class_loglik(X, 0)\n", "\n", "\n", "hist_nb_fit = lambda X, y: HistogramNaiveBayes().fit(X, y)\n", "hist_nb_log_odds = lambda X, model: model.log_odds(X)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2d740a5e", "metadata": {}, "outputs": [], "source": [ "NB_VARIANTS = {\n", " \"Bernoulli\": (bernoulli_nb_fit, bernoulli_nb_log_odds),\n", " \"Gaussian\": (gaussian_nb_fit, gaussian_nb_log_odds),\n", " \"Histogram\": (hist_nb_fit, hist_nb_log_odds),\n", "}\n", "\n", "def evaluate_nb(fit_fn, score_fn, X, y, k=10, seed=42):\n", " skf = StratifiedKFold(n_splits=k, shuffle=True, random_state=seed)\n", " accs = []\n", " for train_idx, test_idx in skf.split(X, y):\n", " model = fit_fn(X[train_idx], y[train_idx])\n", " preds = (score_fn(X[test_idx], model) > 0).astype(int)\n", " accs.append(accuracy_score(y[test_idx], preds))\n", " return np.array(accs)\n", "\n", "nb_cv_results = {}\n", "for name, (fit_fn, score_fn) in NB_VARIANTS.items():\n", " accs = evaluate_nb(fit_fn, score_fn, X_spam, y_spam, k=10)\n", " nb_cv_results[name] = accs\n", " print(f\"{name:10s} 10-fold accuracy: {accs.mean():.4f} +/- {accs.std():.4f}\")\n", "\n", "# sanity check: sklearn's GaussianNB against our from-scratch Gaussian variant (no train/test\n", "# split needed for a single-shot comparison, since both are fit and scored identically below)\n", "Xtr, Xte, ytr, yte = train_test_split(X_spam, y_spam, test_size=0.2, random_state=0, stratify=y_spam)\n", "sk_gnb = GaussianNB... ### TODO_STUDENT ### fit sklearn's GaussianNB as a library baseline for the Gaussian variant\n", "scratch_gnb = gaussian_nb_fit(Xtr, ytr)\n", "print(f\"Gaussian NB sanity check on one 80/20 split — sklearn: \"\n", " f\"{accuracy_score(yte, sk_gnb.predict(Xte)):.4f}, scratch: \"\n", " f\"{accuracy_score(yte, (gaussian_nb_log_odds(Xte, scratch_gnb) > 0).astype(int)):.4f}\")\n", "\n", "# same idea for the Bernoulli variant: sklearn's BernoulliNB only takes a single scalar\n", "# `binarize` threshold, not our per-feature mean vector, so we binarize with our own mu first\n", "# and pass binarize=None (telling sklearn the data is already binary) for a fair comparison\n", "mu_tr = Xtr.mean(axis=0)\n", "Xtr_bin, Xte_bin = (Xtr > mu_tr).astype(int), (Xte > mu_tr).astype(int)\n", "sk_bnb = BernoulliNB... ### TODO_STUDENT ### fit sklearn's BernoulliNB (pre-binarized, alpha=1.0) as a library baseline for the Bernoulli variant\n", "scratch_bnb = bernoulli_nb_fit(Xtr, ytr)\n", "print(f\"Bernoulli NB sanity check on the same 80/20 split — sklearn: \"\n", " f\"{accuracy_score(yte, sk_bnb.predict(Xte_bin)):.4f}, scratch: \"\n", " f\"{accuracy_score(yte, (bernoulli_nb_log_odds(Xte, scratch_bnb) > 0).astype(int)):.4f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "452161c0", "metadata": {}, "outputs": [], "source": [ "# ROC curves + AUC on one held-out split, using each model's log-odds as the decision score\n", "Xtr, Xte, ytr, yte = train_test_split(X_spam, y_spam, test_size=0.2, random_state=42, stratify=y_spam)\n", "\n", "plt.figure(figsize=(5, 5))\n", "roc_data = {}\n", "for name, (fit_fn, score_fn) in NB_VARIANTS.items():\n", " model = fit_fn(Xtr, ytr)\n", " scores = score_fn(Xte, model)\n", " fpr, tpr, thresh = roc_curve(yte, scores)\n", " roc_data[name] = (fpr, tpr, thresh)\n", " plt.plot(fpr, tpr, label=f\"{name} (AUC={auc(fpr, tpr):.3f})\")\n", "plt.plot([0, 1], [0, 1], \"k--\", alpha=0.3)\n", "plt.xlabel(\"False positive rate (legit email flagged as spam)\")\n", "plt.ylabel(\"True positive rate (spam caught)\")\n", "plt.title(\"Naive Bayes ROC — held-out 20%\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5cb7e0d4", "metadata": {}, "outputs": [], "source": [ "# What false-positive rate does each achievable log-odds threshold tau buy us?\n", "best_name = max(nb_cv_results, key=lambda n: nb_cv_results[n].mean())\n", "fpr, tpr, thresh = roc_data[best_name]\n", "table = pd.DataFrame({\"tau\": thresh, \"FPR\": fpr, \"TPR\": tpr})\n", "table.iloc[np.linspace(0, len(table) - 1, 8).astype(int)]\n" ] }, { "cell_type": "markdown", "id": "068a59b8", "metadata": {}, "source": [ "**Choosing a deployment threshold τ.** Using the table above, report and justify what value of\n", "τ you would deploy in a real e-mail spam filter.\n" ] }, { "cell_type": "markdown", "id": "6a0c0194", "metadata": {}, "source": [ "## PROBLEM 3 — EM on generated data   [50 points]\n", "\n", "Own E/M steps for a Gaussian mixture, run from random initial parameters, then checked against\n", "the known generating parameters (EM's component *order* is arbitrary, so components are matched\n", "to the closest true mean before comparing).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c544eef9", "metadata": {}, "outputs": [], "source": [ "def gaussian_mixture_em(X, K, max_iter=300, tol=1e-6, seed=0):\n", " \"\"\"From-scratch EM for a K-component Gaussian mixture. Returns (means, covs, weights,\n", " log-likelihood history).\"\"\"\n", " rng = np.random.RandomState(seed)\n", " n, d = X.shape\n", " means = ... ### TODO_STUDENT ### initialize each mean at a random data point\n", " covs = ... ### TODO_STUDENT ### initialize each covariance at the overall data covariance\n", " weights = np.full(K, 1.0 / K)\n", " loglik_hist = []\n", "\n", " for it in range(max_iter):\n", " # E-step: responsibilities r[i, k] = P(z_i = k | x_i)\n", " dens = np.stack(...) ### TODO_STUDENT ### the (unnormalized) joint density weight_k * N(x_i; mean_k, cov_k) for every point and component\n", " total = dens.sum(axis=1, keepdims=True)\n", " r = ... ### TODO_STUDENT ### normalize each row into a responsibility distribution over the K components\n", " loglik_hist.append(np.sum(np.log(total)))\n", " if it > 0 and abs(loglik_hist[-1] - loglik_hist[-2]) < tol:\n", " break\n", "\n", " # M-step\n", " Nk = r.sum(axis=0)\n", " weights = ... ### TODO_STUDENT ### mixture weights = average responsibility per component\n", " means = ... ### TODO_STUDENT ### each mean = responsibility-weighted average of the data\n", " for k in range(K):\n", " diff = X - means[k]\n", " covs[k] = ... ### TODO_STUDENT ### each covariance = responsibility-weighted scatter matrix\n", "\n", " return means, covs, weights, loglik_hist\n" ] }, { "cell_type": "code", "execution_count": null, "id": "20444338", "metadata": {}, "outputs": [], "source": [ "def match_components(est_means, true_means):\n", " \"\"\"EM's component labels are arbitrary -- greedily match each estimated mean to its\n", " nearest not-yet-used true mean so recovered and true parameters line up for comparison.\"\"\"\n", " remaining = list(range(len(true_means)))\n", " order = []\n", " for m in est_means:\n", " dists = [np.linalg.norm(m - true_means[j]) for j in remaining]\n", " best = remaining.pop(int(np.argmin(dists)))\n", " order.append(best)\n", " return order\n", "\n", "\n", "def compare_to_truth(means, covs, weights, true_means, true_covs, true_weights):\n", " order = match_components(means, true_means)\n", " rows = []\n", " for est_idx, true_idx in enumerate(order):\n", " rows.append({\n", " \"component\": true_idx,\n", " \"true_mean\": np.round(true_means[true_idx], 2).tolist(),\n", " \"est_mean\": np.round(means[est_idx], 2).tolist(),\n", " \"mean_abs_diff\": round(float(np.abs(means[est_idx] - true_means[true_idx]).mean()), 3),\n", " \"true_weight\": round(true_weights[true_idx], 3),\n", " \"est_weight\": round(float(weights[est_idx]), 3),\n", " \"cov_frob_diff\": round(float(np.linalg.norm(covs[est_idx] - true_covs[true_idx])), 3),\n", " })\n", " return pd.DataFrame(rows).set_index(\"component\").sort_index()\n", "\n", "\n", "def plot_gmm_contours(X, means, covs, weights, title):\n", " x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1\n", " y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1\n", " xx, yy = np.meshgrid(np.linspace(x_min, x_max, 150), np.linspace(y_min, y_max, 150))\n", " grid = np.column_stack([xx.ravel(), yy.ravel()])\n", " density = sum(w * multivariate_normal(m, c).pdf(grid) for w, m, c in zip(weights, means, covs))\n", " plt.figure(figsize=(5, 5))\n", " plt.scatter(X[:, 0], X[:, 1], s=3, alpha=0.25)\n", " plt.contour(xx, yy, density.reshape(xx.shape), levels=10)\n", " plt.title(title); plt.xlabel(\"x\"); plt.ylabel(\"y\")\n", " plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c3a90d37", "metadata": {}, "outputs": [], "source": [ "# ---- 2gaussian.txt: true generating parameters (from HW4_26F.html) ----\n", "X_2g = np.loadtxt(\"2gaussian.txt\")\n", "true_means_2g = [np.array([3, 3]), np.array([7, 4])]\n", "true_covs_2g = [np.array([[1, 0], [0, 3]]), np.array([[1, 0.5], [0.5, 1]])]\n", "true_n_2g = [2000, 4000]\n", "true_weights_2g = [n / sum(true_n_2g) for n in true_n_2g]\n", "\n", "means_2g, covs_2g, weights_2g, ll_2g = gaussian_mixture_em(X_2g, K=2, seed=0)\n", "compare_to_truth(means_2g, covs_2g, weights_2g, true_means_2g, true_covs_2g, true_weights_2g)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b703f332", "metadata": {}, "outputs": [], "source": [ "plot_gmm_contours(X_2g, means_2g, covs_2g, weights_2g, \"2gaussian.txt — recovered mixture\")\n", "plt.figure(); plt.plot(ll_2g); plt.xlabel(\"iteration\"); plt.ylabel(\"log-likelihood\")\n", "plt.title(\"2gaussian.txt — EM convergence\"); plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7347a360", "metadata": {}, "outputs": [], "source": [ "# ---- 3gaussian.txt: true generating parameters (from HW4_26F.html) ----\n", "X_3g = np.loadtxt(\"3gaussian.txt\")\n", "true_means_3g = [np.array([3, 3]), np.array([7, 4]), np.array([5, 7])]\n", "true_covs_3g = [np.array([[1, 0], [0, 3]]), np.array([[1, 0.5], [0.5, 1]]), np.array([[1, 0.2], [0.2, 1]])]\n", "true_n_3g = [2000, 3000, 5000]\n", "true_weights_3g = [n / sum(true_n_3g) for n in true_n_3g]\n", "\n", "means_3g, covs_3g, weights_3g, ll_3g = gaussian_mixture_em(X_3g, K=3, seed=0)\n", "compare_to_truth(means_3g, covs_3g, weights_3g, true_means_3g, true_covs_3g, true_weights_3g)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4b767c42", "metadata": {}, "outputs": [], "source": [ "plot_gmm_contours(X_3g, means_3g, covs_3g, weights_3g, \"3gaussian.txt — recovered mixture\")\n", "plt.figure(); plt.plot(ll_3g); plt.xlabel(\"iteration\"); plt.ylabel(\"log-likelihood\")\n", "plt.title(\"3gaussian.txt — EM convergence\"); plt.show()\n", "\n", "# library sanity check: sklearn's GaussianMixture, same K, should reach the same means/log-likelihood\n", "for name, X, K in [(\"2gaussian\", X_2g, 2), (\"3gaussian\", X_3g, 3)]:\n", " skgmm = GaussianMixture... ### TODO_STUDENT ### fit sklearn's GaussianMixture as a library baseline, same K\n", " print(f\"{name}: sklearn means =\\n{np.round(skgmm.means_, 2)}\")\n" ] }, { "cell_type": "markdown", "id": "d0981adf", "metadata": {}, "source": [ "Compare the true-vs-recovered tables and contour plots above, and sklearn's `GaussianMixture`\n", "means, to your own results, and discuss.\n" ] }, { "cell_type": "markdown", "id": "81c1c819", "metadata": {}, "source": [ "## PROBLEM 4 — EM with a Mixture of Binomials   [50 points]\n", "\n", "No off-the-shelf library implements a Binomial-mixture EM, so this problem is validated the same\n", "way as Problem 3: generate data from known (p, r, π), then check that EM recovers them.\n", "\n", "**A) Generate the data.**\n" ] }, { "cell_type": "code", "execution_count": null, "id": "aef2ca40", "metadata": {}, "outputs": [], "source": [ "def generate_coin_flip_data(p, r, pi, K=10, M=200, seed=0):\n", " \"\"\"Simulate M sequences of K flips from a 2-coin mixture: coin 1 (head-probability p) is\n", " chosen with probability pi, coin 2 (head-probability r) otherwise. Returns the M x K matrix\n", " of 0/1 outcomes (the coin identity itself is not returned -- EM has to recover it).\"\"\"\n", " rng = np.random.RandomState(seed)\n", " coin1 = ... ### TODO_STUDENT ### decide which coin generated each sequence, with probability pi for coin 1\n", " head_prob = np.where(coin1, p, r)\n", " flips = ... ### TODO_STUDENT ### simulate K flips of the chosen coin, for every sequence at once\n", " return flips\n", "\n", "flips_demo = generate_coin_flip_data(p=0.75, r=0.4, pi=0.8, K=10, M=5, seed=0)\n", "print(\"example generated sequences (rows = sequences, columns = the K=10 flips):\")\n", "print(flips_demo)\n" ] }, { "cell_type": "markdown", "id": "f2164f80", "metadata": {}, "source": [ "**B) Recover the parameters.**" ] }, { "cell_type": "code", "execution_count": null, "id": "658cd5a0", "metadata": {}, "outputs": [], "source": [ "def binomial_mixture_em(flips, n_components=2, max_iter=300, tol=1e-8, seed=0):\n", " \"\"\"From-scratch EM for a mixture of Binomial(K, p_k) components -- one scalar\n", " head-probability per coin, shared across all K flips in a sequence (not a separate parameter\n", " per flip position: within one sequence every flip comes from the same coin).\"\"\"\n", " rng = np.random.RandomState(seed)\n", " M, K = flips.shape\n", " heads = flips.sum(axis=1) # sufficient statistic per sequence\n", " p = rng.uniform(0.2, 0.8, n_components)\n", " weights = np.full(n_components, 1.0 / n_components)\n", " loglik_hist = []\n", "\n", " for it in range(max_iter):\n", " # E-step. The shared C(K, heads_i) binomial-coefficient factor cancels in the\n", " # responsibility ratio, so it's dropped here -- and is never added back, including in\n", " # loglik_hist below. That's fine for parameter estimation (it's a per-sequence additive\n", " # constant, identical across components and iterations, so it affects neither\n", " # responsibilities nor convergence) but means loglik_hist is the mixture log-density up\n", " # to that missing constant, not the true observed-data log-likelihood -- only the\n", " # convergence *plot*'s absolute y-axis scale is affected.\n", " log_lik = ... ### TODO_STUDENT ### per-sequence, per-coin Binomial log-likelihood (up to the shared C(K,heads) factor)\n", " m = log_lik.max(axis=1, keepdims=True)\n", " unnorm = weights[None, :] * np.exp(log_lik - m)\n", " total = unnorm.sum(axis=1, keepdims=True)\n", " r = ... ### TODO_STUDENT ### normalize into responsibilities over the n_components coins\n", " loglik_hist.append(float(np.sum(np.log(total) + m)))\n", " if it > 0 and abs(loglik_hist[-1] - loglik_hist[-2]) < tol:\n", " break\n", "\n", " # M-step\n", " Nk = r.sum(axis=0)\n", " weights = ... ### TODO_STUDENT ### update mixture weights\n", " p = ... ### TODO_STUDENT ### update each coin's head-probability: responsibility-weighted fraction of heads\n", "\n", " return p, weights, loglik_hist\n" ] }, { "cell_type": "code", "execution_count": null, "id": "af115527", "metadata": {}, "outputs": [], "source": [ "# Try several (p, r, pi) settings; match recovered coins to the true ones by sorting on head\n", "# probability (with only two components, sorting both true and estimated pairs ascending is an\n", "# easy, unambiguous match -- and keeps the \"weights\" columns aligned with the \"(p, r)\" columns).\n", "settings = [(0.75, 0.4, 0.8), (0.6, 0.3, 0.5), (0.9, 0.1, 0.3)]\n", "rows = []\n", "for p_true, r_true, pi_true in settings:\n", " flips = generate_coin_flip_data(p_true, r_true, pi_true, K=10, M=300, seed=1)\n", " p_est, w_est, ll = binomial_mixture_em(flips, n_components=2, seed=0)\n", "\n", " true_p, true_w = np.array([p_true, r_true]), np.array([pi_true, 1 - pi_true])\n", " true_order, est_order = np.argsort(true_p), np.argsort(p_est)\n", " rows.append({\n", " \"true (p, r)\": tuple(true_p[true_order].round(3)),\n", " \"est (p, r)\": tuple(p_est[est_order].round(3)),\n", " \"true weights\": tuple(true_w[true_order].round(3)),\n", " \"est weights\": tuple(w_est[est_order].round(3)),\n", " \"iterations\": len(ll),\n", " })\n", "\n", "binomial_em_results = pd.DataFrame(rows)\n", "binomial_em_results\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ad71a52e", "metadata": {}, "outputs": [], "source": [ "# convergence plot for the first setting above\n", "flips = generate_coin_flip_data(*settings[0], K=10, M=300, seed=1)\n", "_, _, ll = binomial_mixture_em(flips, n_components=2, seed=0)\n", "plt.figure(); plt.plot(ll); plt.xlabel(\"iteration\"); plt.ylabel(\"log-likelihood\")\n", "plt.title(f\"Binomial mixture EM convergence — p={settings[0][0]}, r={settings[0][1]}, pi={settings[0][2]}\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "2f3526e8", "metadata": {}, "source": [ "Compare the recovered (p, r) and π to the true generating values above, and discuss.\n", "\n", "**[optional extension, no credit] C) T coins instead of two.** `binomial_mixture_em` already\n", "takes `n_components` as a parameter, so generalizing to T coins is a one-line change -- a quick\n", "demonstration with T=3:\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f1a3f825", "metadata": {}, "outputs": [], "source": [ "p_true_T = [0.85, 0.5, 0.15]\n", "w_true_T = [0.3, 0.4, 0.3]\n", "flips_T = np.zeros((300, 10), dtype=int)\n", "rng = np.random.RandomState(2)\n", "coin_idx = rng.choice(3, size=300, p=w_true_T)\n", "for i, k in enumerate(coin_idx):\n", " flips_T[i] = rng.binomial(1, p_true_T[k], size=10)\n", "\n", "p_est_T, w_est_T, _ = binomial_mixture_em(flips_T, n_components=3, seed=0)\n", "order_T = np.argsort(p_est_T)\n", "print(\"true p (sorted): \", np.round(sorted(p_true_T), 3))\n", "print(\"est p (sorted): \", np.round(p_est_T[order_T], 3))\n" ] }, { "cell_type": "markdown", "id": "5beb5bc1", "metadata": {}, "source": [ "## PROBLEM 5 — Bayesian Linear Regression: Ridge as MAP   [60 points]\n", "\n", "Bridges back to HW1's closed-form Ridge regression. Model: y = wᵀx + ε, ε ~ N(0, σ²), with a\n", "Gaussian prior w ~ N(0, τ²I). Throughout, `X_train_housing` / `X_test_housing` are already\n", "standardized (unit variance, zero training mean) and `y` is centered by its training mean — the\n", "centering half of that is the same trick sklearn's own `Ridge(fit_intercept=True)` uses internally\n", "to avoid needing to special-case an unpenalized intercept term (the standardizing half is an extra\n", "step we take on top, so one isotropic prior variance τ² is sensible across features of very\n", "different scales). See `../lecture_notes/bayesian_ridge_regression.pdf` for the full derivation.\n", "\n", "Parts A and B below are required. Part C (further down) is optional, no credit.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eca2737c", "metadata": {}, "outputs": [], "source": [ "# ---- Housing (regression), pre-split train/test files, standardized like HW1 ----\n", "train_housing = np.loadtxt(f\"{DATA}/housing_train.txt\")\n", "test_housing = np.loadtxt(f\"{DATA}/housing_test.txt\")\n", "X_train_housing, y_train_housing = train_housing[:, :-1], train_housing[:, -1]\n", "X_test_housing, y_test_housing = test_housing[:, :-1], test_housing[:, -1]\n", "\n", "housing_scaler = StandardScaler().fit(X_train_housing)\n", "X_train_housing = housing_scaler.transform(X_train_housing)\n", "X_test_housing = housing_scaler.transform(X_test_housing)\n", "\n", "y_train_mean = y_train_housing.mean()\n", "yc_train_housing = y_train_housing - y_train_mean # centered target -- removes the need for a bias term\n", "\n", "print(\"Housing train:\", X_train_housing.shape, \" test:\", X_test_housing.shape)\n" ] }, { "cell_type": "markdown", "id": "5b639774", "metadata": {}, "source": [ "**A) MAP = Ridge.** Posterior mean (the MAP estimate) vs. sklearn's `Ridge`, at a fixed λ." ] }, { "cell_type": "code", "execution_count": null, "id": "e88592c3", "metadata": {}, "outputs": [], "source": [ "def bayesian_ridge_posterior(X, y_centered, sigma2, lam):\n", " \"\"\"Closed-form Gaussian posterior N(mean_w, cov_w) over the (centered-data) regression\n", " weights: precision A = (1/sigma2) X^T X + (1/tau2) I, with lam = sigma2/tau2.\"\"\"\n", " A = ... ### TODO_STUDENT ### posterior precision (up to the 1/sigma2 factor): X^T X + lambda*I\n", " cov_w = sigma2 * np.linalg... ### TODO_STUDENT ### posterior covariance Sigma_w = sigma2 * (X^T X + lambda I)^-1\n", " mean_w = np.linalg... ### TODO_STUDENT ### posterior mean mu_w -- identical to the Ridge solution at this lambda\n", " return mean_w, cov_w\n", "\n", "\n", "lam_check = 5.0\n", "sk_ridge_check = Ridge... ### TODO_STUDENT ### fit sklearn's Ridge as a library check for Part A's \"MAP = Ridge\" claim\n", "mean_w_check, _ = bayesian_ridge_posterior(X_train_housing, yc_train_housing, sigma2=1.0, lam=lam_check)\n", "\n", "print(f\"lambda = {lam_check}\")\n", "print(\"max |sklearn Ridge coef - posterior mean (MAP)| =\",\n", " np.max(np.abs(sk_ridge_check.coef_ - mean_w_check)))\n", "print(\"sklearn intercept:\", round(sk_ridge_check.intercept_, 4),\n", " \" vs. training-set y-mean:\", round(y_train_mean, 4))\n" ] }, { "cell_type": "markdown", "id": "92e3064c", "metadata": {}, "source": [ "**B) The full posterior, not just its mode.** Predictive distribution at each test point:\n", "mean from the posterior mean, variance = leftover noise σ² plus the model's own parameter\n", "uncertainty xᵀΣ_w x.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ae4e7386", "metadata": {}, "outputs": [], "source": [ "def predictive_distribution(X, mean_w, cov_w, sigma2):\n", " pred_mean = X @ mean_w\n", " pred_var = ... ### TODO_STUDENT ### predictive variance = leftover noise sigma2 + parameter uncertainty x^T Sigma_w x\n", " return pred_mean, pred_var\n", "\n", "\n", "# sigma^2: a plug-in estimate from the (unregularized) least-squares residual variance on\n", "# the training data, held fixed while lambda varies below\n", "w_ols, *_ = np.linalg.lstsq(X_train_housing, yc_train_housing, rcond=None)\n", "sigma2_hat = np.var(yc_train_housing - X_train_housing @ w_ols)\n", "print(f\"sigma^2 (from OLS residual variance): {sigma2_hat:.3f}\")\n", "\n", "lam_demo = 5.0 # revisited properly via the evidence in the optional Part C below\n", "mean_w, cov_w = bayesian_ridge_posterior(X_train_housing, yc_train_housing, sigma2_hat, lam_demo)\n", "pred_mean, pred_var = predictive_distribution(X_test_housing, mean_w, cov_w, sigma2_hat)\n", "pred_mean = pred_mean + y_train_mean # back into original Housing-price units\n", "pred_std = np.sqrt(pred_var)\n", "\n", "lower, upper = pred_mean - 1.96 * pred_std, pred_mean + 1.96 * pred_std\n", "coverage = np.mean((y_test_housing >= lower) & (y_test_housing <= upper))\n", "print(f\"95% predictive-interval empirical coverage on the Housing test set: {coverage:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "08321f6d", "metadata": {}, "source": [ "**[optional, no credit] C) Choosing λ without touching a validation set.** The evidence\n", "log p(y|X) integrates w out entirely: y ~ N(0, C) with C = σ²I + τ²XXᵀ (a sum of two independent\n", "zero-mean Gaussians, Xw ~ N(0, τ²XXᵀ) and ε ~ N(0, σ²I)). Sweep λ, and compare the\n", "evidence-maximizing λ to the λ that k-fold cross validation on the *training* data would pick --\n", "not a single held-out test split, which turns out to be a surprisingly unreliable comparison\n", "target (see the lecture note for the full worked-through example).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7f7a8361", "metadata": {}, "outputs": [], "source": [ "def log_evidence(X, y_centered, sigma2, lam):\n", " n = X.shape[0]\n", " tau2 = sigma2 / lam\n", " C = ... ### TODO_STUDENT ### marginal covariance of y after integrating out w: sigma2*I + tau2*X X^T\n", " sign, logdet = np.linalg.slogdet(C)\n", " quad = y_centered @ np.linalg... ### TODO_STUDENT ### the quadratic form y^T C^-1 y (solved, not explicitly inverted)\n", " return -0.5 * (quad + logdet + n * np.log(2 * np.pi))\n", "\n", "\n", "def kfold_cv_mse(X, y_centered, lam, k=10, seed=0):\n", " \"\"\"Mean held-out MSE for one lambda, k-fold cross-validated *within the training set*\n", " (never touches the real test set) -- the standard, statistically sound way to estimate\n", " generalization error, unlike a single small held-out split.\"\"\"\n", " kf = KFold(n_splits=k, shuffle=True, random_state=seed)\n", " fold_mses = []\n", " for tr_idx, val_idx in kf.split(X):\n", " mean_w_fold, _ = bayesian_ridge_posterior(X[tr_idx], y_centered[tr_idx], sigma2_hat, lam)\n", " fold_mses.append(np.mean((y_centered[val_idx] - X[val_idx] @ mean_w_fold) ** 2))\n", " return np.mean(fold_mses)\n", "\n", "\n", "lambdas = np.logspace(-2, 3, 40) # same range as HW1's Ridge regularization-path sweep\n", "evidences = [log_evidence(X_train_housing, yc_train_housing, sigma2_hat, lam) for lam in lambdas]\n", "lam_evidence = lambdas[np.argmax(evidences)]\n", "\n", "cv_mses = [kfold_cv_mse(X_train_housing, yc_train_housing, lam) for lam in lambdas]\n", "lam_cv = lambdas[np.argmin(cv_mses)]\n", "\n", "test_mses = []\n", "for lam in lambdas:\n", " mean_w_l, _ = bayesian_ridge_posterior(X_train_housing, yc_train_housing, sigma2_hat, lam)\n", " pred_l = X_test_housing @ mean_w_l + y_train_mean\n", " test_mses.append(np.mean((y_test_housing - pred_l) ** 2))\n", "lam_test_best = lambdas[np.argmin(test_mses)]\n", "\n", "# library check: sklearn's BayesianRidge runs the real thing -- joint type-II ML / MacKay's\n", "# \"evidence procedure\" over BOTH sigma2 and tau2 (we only swept lambda at a fixed, plug-in\n", "# sigma2) -- if your from-scratch evidence sweep is correct, its argmax should land close to\n", "# sklearn's independently-optimized answer\n", "sk_bayes_ridge = BayesianRidge... ### TODO_STUDENT ### fit sklearn's BayesianRidge as a library check for the evidence-based lambda\n", "lam_sklearn = sk_bayes_ridge.lambda_ / sk_bayes_ridge.alpha_ # sklearn's (weight precision)/(noise precision) = our sigma2/tau2 = lambda\n", "\n", "print(f\"lambda maximizing the evidence (from-scratch, plug-in sigma2, no held-out data used): {lam_evidence:.3g}\")\n", "print(f\"lambda implied by sklearn's BayesianRidge (joint evidence maximization): {lam_sklearn:.3g}\")\n", "print(f\"lambda minimizing 10-fold CV MSE (training data only, the fair comparison): {lam_cv:.3g}\")\n", "print(f\"lambda minimizing the single 74-point test-split MSE (kept for contrast -- see below): {lam_test_best:.3g}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "bb0e4f63", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "axes[0].plot(lambdas, evidences, \"o-\")\n", "axes[0].axvline(lam_evidence, color=\"r\", linestyle=\"--\", label=f\"scratch argmax = {lam_evidence:.2g}\")\n", "axes[0].axvline(lam_sklearn, color=\"g\", linestyle=\":\", label=f\"sklearn BayesianRidge = {lam_sklearn:.2g}\")\n", "axes[0].set_xscale(\"log\"); axes[0].set_xlabel(\"lambda\"); axes[0].set_ylabel(\"log evidence\")\n", "axes[0].set_title(\"Evidence vs. lambda (never touches held-out data)\"); axes[0].legend()\n", "\n", "axes[1].plot(lambdas, cv_mses, \"o-\", color=\"seagreen\", label=\"10-fold CV (training data)\")\n", "axes[1].axvline(lam_cv, color=\"seagreen\", linestyle=\"--\", label=f\"CV argmin = {lam_cv:.2g}\")\n", "axes[1].plot(lambdas, test_mses, \"o-\", color=\"darkorange\", alpha=0.6, label=\"single 74-pt test split\")\n", "axes[1].axvline(lam_test_best, color=\"darkorange\", linestyle=\":\", label=f\"test argmin = {lam_test_best:.2g}\")\n", "axes[1].set_xscale(\"log\"); axes[1].set_xlabel(\"lambda\"); axes[1].set_ylabel(\"MSE\")\n", "axes[1].set_title(\"CV vs. single-split MSE\"); axes[1].legend(fontsize=8)\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "f498d402", "metadata": {}, "source": [ "Discuss: do the evidence-based and cross-validated λ choices land in the same region? How about\n", "the single test-split argmin? What do you gain (and what do you give up) by picking λ from the\n", "evidence instead of cross-validating -- note it never touches a validation fold, or the test set,\n", "at all.\n" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }