{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"Lecture 07: Linear Regression (2) – CS 189, Fall 2026
\n",
"\n",
"\n",
"In this lecture we take the probabilistic view: least squares as maximum likelihood, what the choice of noise model buys, regularization restated as a prior belief, and the bias-variance decomposition.\n",
"\n",
"*Reference: Bishop & Bishop,* Deep Learning: Foundations and Concepts, *§4.1.2–4.1.3, §4.2, §4.3.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:50.257617Z",
"iopub.status.busy": "2026-09-15T03:39:50.257404Z",
"iopub.status.idle": "2026-09-15T03:39:51.528496Z",
"shell.execute_reply": "2026-09-15T03:39:51.527187Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import scipy.stats as st\n",
"from scipy.optimize import minimize\n",
"from numpy.polynomial import legendre as L\n",
"from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.metrics import mean_squared_error\n",
"from plotly.subplots import make_subplots\n",
"import plotly.graph_objects as go\n",
"\n",
"import warnings\n",
"from sklearn.exceptions import ConvergenceWarning\n",
"warnings.filterwarnings(\"ignore\", category=ConvergenceWarning)\n",
"warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
"\n",
"import plotly.io as pio\n",
"pio.renderers.default = \"notebook_connected\"\n",
"\n",
"np.random.seed(42)\n",
"\n",
"C_DATA, C_FIT, C_ALT, C_RESID, C_SPAN = \"#003262\", \"#FDB515\", \"#C4820E\", \"#D55E00\", \"#00553A\"\n",
"plt.rcParams.update({\"figure.figsize\": (8, 5), \"axes.grid\": True, \"grid.alpha\": 0.3,\n",
" \"font.size\": 13, \"axes.titlesize\": 15, \"axes.labelsize\": 14})\n",
"\n",
"def build_Phi(x, D, bias=True):\n",
" \"Polynomial design matrix WITH the bias column, written out explicitly.\"\n",
" x = np.asarray(x).ravel()\n",
" cols = [np.ones_like(x)] if bias else []\n",
" cols += [x ** j for j in range(1, D)]\n",
" return np.column_stack(cols)\n",
"\n",
"# Shared course dataset, identical to Lecture 06.\n",
"_rng_course = np.random.default_rng(189)\n",
"n = 100\n",
"x_s = np.sort(_rng_course.random(n) * 2 - 1)\n",
"t_s = np.sin(5 * x_s) + 0.1 * _rng_course.standard_normal(n)\n",
"x_dense = np.linspace(-1, 1, 400)\n",
"\n",
"lin = LinearRegression().fit(x_s[:, None], t_s)\n",
"mse_lin = np.mean((t_s - lin.predict(x_s[:, None])) ** 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 0. **Where we left off**\n",
"\n",
"Lecture 06 fitted models and judged them by held-out error. Two questions were left open: where the squared error came from in the first place, and why the penalty $\\lambda\\lVert\\mathbf{w}\\rVert^2$ should be believed rather than merely observed to work."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. **Model Complexity**\n",
"\n",
"Now turn the dial on model complexity, with a held-out test set this time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"i_tr, i_te = train_test_split(np.arange(n), test_size=0.55, random_state=0)\n",
"\n",
"print(f\"{len(i_tr)} training points, {len(i_te)} test points\\n\")\n",
"print(f\"{'degree':>7} {'train MSE':>12} {'TEST MSE':>12} {'max |w_j|':>12}\")\n",
"for d in [1, 3, 5, 9, 16]:\n",
" P_tr, P_te = build_Phi(x_s[i_tr], d + 1), build_Phi(x_s[i_te], d + 1)\n",
" w = np.linalg.lstsq(P_tr, t_s[i_tr], rcond=None)[0]\n",
" print(f\"{d:>7} {np.mean((t_s[i_tr] - P_tr @ w)**2):>12.5f} \"\n",
" f\"{np.mean((t_s[i_te] - P_te @ w)**2):>12.5f} {np.abs(w).max():>12.1f}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds = list(range(1, 17))\n",
"TR, TE, MW = [], [], []\n",
"for d in ds:\n",
" P_tr, P_te = build_Phi(x_s[i_tr], d + 1), build_Phi(x_s[i_te], d + 1)\n",
" w = np.linalg.lstsq(P_tr, t_s[i_tr], rcond=None)[0]\n",
" TR.append(np.mean((t_s[i_tr] - P_tr @ w) ** 2))\n",
" TE.append(np.mean((t_s[i_te] - P_te @ w) ** 2))\n",
" MW.append(np.abs(w).max())\n",
"k_best = int(np.argmin(TE))\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(15, 5))\n",
"\n",
"axs[0].scatter(x_s, t_s, s=35, alpha=0.45, color=C_DATA, zorder=3, label=\"data\")\n",
"for d, col in zip([1, 5, 16], [C_ALT, C_SPAN, C_RESID]):\n",
" w = np.linalg.lstsq(build_Phi(x_s[i_tr], d + 1), t_s[i_tr], rcond=None)[0]\n",
" axs[0].plot(x_dense, build_Phi(x_dense, d + 1) @ w, lw=3, color=col, label=f\"degree {d}\")\n",
"axs[0].set(xlabel=\"$x$\", ylabel=\"$t$\", ylim=(-2, 2), title=\"Fitted functions\")\n",
"axs[0].legend(fontsize=13)\n",
"\n",
"axs[1].semilogy(ds, TR, \"o-\", lw=3.5, ms=9, color=C_DATA, label=\"training MSE\")\n",
"axs[1].semilogy(ds, TE, \"s-\", lw=3.5, ms=9, color=C_RESID, label=\"test MSE\")\n",
"axs[1].axvline(ds[k_best], color=C_ALT, ls=\"--\", lw=3,\n",
" label=f\"lowest test error: degree {ds[k_best]}\")\n",
"axs[1].set(xlabel=\"polynomial degree\", ylabel=\"MSE\", xticks=ds[::2],\n",
" title=\"Train and test error\")\n",
"axs[1].legend(fontsize=13)\n",
"\n",
"plt.tight_layout(); plt.show()\n",
"\n",
"print(f\"lowest test MSE {TE[k_best]:.5f} at degree {ds[k_best]}\")\n",
"print(f\"degree 16: train {TR[-1]:.5f}, test {TE[-1]:.5f}, max |w_j| {MW[-1]:.1f}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, axs = plt.subplots(1, 2, figsize=(15, 5))\n",
"\n",
"axs[0].scatter(x_s, t_s, s=35, alpha=0.45, color=C_DATA, zorder=3, label=\"data\")\n",
"for d, col in zip([1, 5, 16], [C_ALT, C_SPAN, C_RESID]):\n",
" w = np.linalg.lstsq(build_Phi(x_s[i_tr], d + 1), t_s[i_tr], rcond=None)[0]\n",
" axs[0].plot(x_dense, build_Phi(x_dense, d + 1) @ w, lw=3, color=col, label=f\"degree {d}\")\n",
"axs[0].set(xlabel=\"$x$\", ylabel=\"$t$\", ylim=(-2, 2),\n",
" title=\"Higher degree, lower training error\")\n",
"axs[0].legend(fontsize=13)\n",
"\n",
"axs[1].semilogy(ds, TR, \"o-\", lw=3.5, ms=9, color=C_DATA, label=\"training MSE\")\n",
"axs[1].semilogy(ds, MW, \"s--\", lw=3.5, ms=9, color=C_RESID, label=r\"max $|w_j|$\")\n",
"axs[1].set(xlabel=\"polynomial degree\", xticks=ds[::2],\n",
" title=\"...and exploding coefficients\")\n",
"axs[1].legend(fontsize=13)\n",
"\n",
"plt.tight_layout(); plt.show()\n",
"\n",
"print(\"max |w_j| by degree:\")\n",
"print(\" \" + \" \".join(f\"{d}:{m:,.0f}\" for d, m in zip(ds, MW)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. **Regularized Least Squares**\n",
"\n",
"Training error falls, test error turns around, and the coefficients explode. Regularization attacks the coefficients directly:\n",
"\n",
"$$\n",
"E(\\mathbf{w}) = \\underbrace{\\tfrac{1}{2}\\lVert \\mathbf{t} - \\Phi\\mathbf{w}\\rVert^2}_{E_D(\\mathbf{w})}\n",
"\\;+\\;\\lambda\\, \\underbrace{E_W(\\mathbf{w})}_{\\text{penalty}}\n",
"$$\n",
"\n",
"| | penalty | closed form? | effect on $\\mathbf{w}$ |\n",
"|---|---|---|---|\n",
"| **Ridge (L2)** | $\\tfrac12\\lVert\\mathbf{w}\\rVert_2^2$ | yes: $(\\Phi^\\top\\Phi + \\lambda I)^{-1}\\Phi^\\top\\mathbf{t}$ | shrinks all coefficients smoothly |\n",
"| **Lasso (L1)** | $\\lVert\\mathbf{w}\\rVert_1$ | no | drives some coefficients **exactly** to zero |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2.1 Coefficient paths"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"DEG = 10\n",
"i_tr2, i_te2 = train_test_split(np.arange(n), test_size=0.85, random_state=189)\n",
"# Standardize the polynomial columns: penalties are not scale-invariant, and x^10 is\n",
"# numerically tiny next to x^1. Skipping this is the most common bug in ridge/lasso demos.\n",
"Ptr_raw = build_Phi(x_s[i_tr2], DEG + 1, bias=False)\n",
"mu, sd = Ptr_raw.mean(0), Ptr_raw.std(0)\n",
"Ptr = (Ptr_raw - mu) / sd\n",
"Pte = (build_Phi(x_s[i_te2], DEG + 1, bias=False) - mu) / sd\n",
"ttr, tte = t_s[i_tr2], t_s[i_te2]\n",
"\n",
"lambdas = np.logspace(-6, 2, 40)\n",
"\n",
"def path(Model, penalty):\n",
" coefs, tr, te = [], [], []\n",
" for lam in lambdas:\n",
" m = Model(alpha=lam, fit_intercept=True, max_iter=500_000, tol=1e-8).fit(Ptr, ttr)\n",
" coefs.append(m.coef_.ravel())\n",
" tr.append(mean_squared_error(ttr, m.predict(Ptr)))\n",
" te.append(mean_squared_error(tte, m.predict(Pte)))\n",
" return np.array(coefs), np.array(tr), np.array(te)\n",
"\n",
"ridge_c, ridge_tr, ridge_te = path(Ridge, None)\n",
"lasso_c, lasso_tr, lasso_te = path(Lasso, None)\n",
"\n",
"print(f\"n_train = {len(i_tr2)} points, degree {DEG} -> {DEG + 1} parameters.\")\n",
"print(f\"unregularized (lambda -> 0) test MSE: {ridge_te[0]:.4f}\")\n",
"print(f\"ridge: best test MSE {ridge_te.min():.4f} at lambda = {lambdas[ridge_te.argmin()]:.4g}\"\n",
" f\" ({ridge_te[0]/ridge_te.min():.0f}x better)\")\n",
"print(f\"lasso: best test MSE {lasso_te.min():.4f} at lambda = {lambdas[lasso_te.argmin()]:.4g}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, axs = plt.subplots(2, 2, figsize=(13, 9))\n",
"for row, (coefs, tr, te, name) in enumerate(\n",
" [(ridge_c, ridge_tr, ridge_te, \"Ridge (L2)\"),\n",
" (lasso_c, lasso_tr, lasso_te, \"Lasso (L1)\")]):\n",
" ax = axs[row, 0]\n",
" for j in range(coefs.shape[1]):\n",
" ax.plot(lambdas, coefs[:, j], lw=2, label=f\"degree {j+1}\")\n",
" ax.axhline(0, color=\"k\", lw=1)\n",
" ax.set(xscale=\"log\", xlabel=r\"$\\lambda$\", ylabel=\"coefficient\",\n",
" title=f\"{name}: coefficient paths\")\n",
" ax.legend(fontsize=8, ncol=2)\n",
"\n",
" ax = axs[row, 1]\n",
" ax.plot(lambdas, tr, lw=3, color=C_DATA, label=\"train MSE\")\n",
" ax.plot(lambdas, te, lw=3, color=C_RESID, label=\"test MSE\")\n",
" ax.axvline(lambdas[te.argmin()], color=C_ALT, ls=\"--\", lw=2.5,\n",
" label=rf\"best $\\lambda$ = {lambdas[te.argmin()]:.2g}\")\n",
" ax.set(xscale=\"log\", yscale=\"log\", xlabel=r\"$\\lambda$\", ylabel=\"MSE\",\n",
" title=f\"{name}: train vs test\")\n",
" ax.legend(fontsize=10)\n",
"plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2.2 Lasso sets coefficients to zero; ridge does not\n",
"\n",
"Index convention: with the bias column excluded, column $j$ holds $x^{j+1}$, so the first row is **degree 1**, not degree 0."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show = [1e-6, 7e-4, 1e-2, 1e-1, 1e0]\n",
"for name, C in ((\"LASSO\", lasso_c), (\"RIDGE\", ridge_c)):\n",
" rows = {}\n",
" for lam in show:\n",
" j = int(np.argmin(np.abs(lambdas - lam)))\n",
" rows[f\"lam={lambdas[j]:.1e}\"] = C[j]\n",
" tbl = pd.DataFrame(rows, index=[f\"degree {d}\" for d in range(1, DEG + 1)])\n",
" print(f\"{name} coefficients\\n\"); print(tbl.round(3).to_string())\n",
" print(\"exact zeros per lambda:\", {c: int((tbl[c] == 0).sum()) for c in tbl.columns}, \"\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2.3 The constraint picture\n",
"\n",
"The penalized problem is equivalent to minimizing $E_D(\\mathbf{w})$ subject to $E_W(\\mathbf{w}) \\le c$ for some $c(\\lambda)$. The contours of $E_D$ grow until they first touch the constraint region. The L2 region is a **circle**, smooth everywhere; the L1 region is a **diamond** whose corners lie on the axes, and a corner is a point where one coordinate is exactly zero."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"W0 = np.array([3.2, 0.9])\n",
"th = np.deg2rad(25)\n",
"Rot = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])\n",
"A_MAT = Rot @ np.diag([6.0, 1.0]) @ Rot.T * 2.5\n",
"\n",
"def ridge_solution(lam, A=A_MAT, w0=W0):\n",
" w = np.linalg.solve(A + lam * np.eye(2), A @ w0)\n",
" return w, 0.5 * (w - w0) @ A @ (w - w0), 0.5 * (w @ w)\n",
"\n",
"def lasso_solution(lam, A=A_MAT, w0=W0, iters=4000):\n",
" \"Coordinate descent with soft-thresholding: reaches exact zeros, unlike Nelder-Mead.\"\n",
" w, b = np.zeros(2), A @ w0\n",
" for _ in range(iters):\n",
" for j in range(2):\n",
" rho = b[j] - A[j] @ w + A[j, j] * w[j]\n",
" w[j] = np.sign(rho) * max(abs(rho) - lam, 0.0) / A[j, j]\n",
" return w, 0.5 * (w - w0) @ A @ (w - w0), np.sum(np.abs(w))\n",
"\n",
"for lam in [0.0, 2.0, 6.0, 12.0]:\n",
" wr, wl = ridge_solution(lam)[0], lasso_solution(lam)[0]\n",
" print(f\"lambda={lam:5.1f} ridge w = [{wr[0]:6.3f} {wr[1]:6.3f}] \"\n",
" f\"lasso w = [{wl[0]:6.3f} {wl[1]:6.3f}]\"\n",
" + (\" <- w2 is exactly 0\" if wl[1] == 0.0 else \"\"))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def regularization_figure(solver, kind):\n",
" gx = np.linspace(W0[0] - 10, W0[0] + 10, 401)\n",
" gy = np.linspace(W0[1] - 10, W0[1] + 10, 401)\n",
" GX, GY = np.meshgrid(gx, gy)\n",
" U = np.stack([GX - W0[0], GY - W0[1]], axis=-1)\n",
" Z = 0.5 * np.einsum(\"...i,...i\", U, U @ A_MAT)\n",
" zmax = float(np.percentile(Z, 95))\n",
"\n",
" lams_curve = np.linspace(0.0, 15.0, 400)\n",
" ED, EW = np.array([[solver(l)[1], solver(l)[2]] for l in lams_curve]).T\n",
" Etot = ED + lams_curve * EW\n",
" ymax = float(Etot.max()) * 1.05\n",
"\n",
" def shape(lam):\n",
" w = solver(lam)[0]\n",
" if kind == \"ridge\":\n",
" c = np.linalg.norm(w)\n",
" a = np.linspace(0, 2 * np.pi, 400)\n",
" return c * np.cos(a), c * np.sin(a)\n",
" c = np.linalg.norm(w, ord=1)\n",
" return c * np.array([1, 0, -1, 0, 1]), c * np.array([0, 1, 0, -1, 0])\n",
"\n",
" norm_lbl = \"||w||₂ = c(λ)\" if kind == \"ridge\" else \"||w||₁ = c(λ)\"\n",
" fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45],\n",
" subplot_titles=(f\"{kind.capitalize()}: data contours and constraint region\",\n",
" \"Error decomposition vs λ\"))\n",
" fig.add_trace(go.Contour(x=gx, y=gy, z=np.clip(Z, 0, zmax), zmin=0, zmax=zmax,\n",
" colorscale=\"Blues\", reversescale=True, showscale=False, opacity=0.96,\n",
" contours=dict(start=0.01 * zmax, end=0.99 * zmax,\n",
" size=0.98 * zmax / 20, showlines=False)), row=1, col=1)\n",
" fig.add_trace(go.Scatter(x=[W0[0]], y=[W0[1]], mode=\"markers\",\n",
" marker=dict(symbol=\"star\", size=16, color=\"crimson\"),\n",
" name=\"unregularized optimum\"), row=1, col=1)\n",
" for y, nm in [(ED, \"E_D\"), (lams_curve * EW, \"λ·E_W\"), (Etot, \"E\")]:\n",
" fig.add_trace(go.Scatter(x=lams_curve, y=y, mode=\"lines\", line=dict(width=3), name=nm),\n",
" row=1, col=2)\n",
"\n",
" cx, cy = shape(0.0)\n",
" dyn = [go.Scatter(x=cx, y=cy, mode=\"lines\",\n",
" line=dict(width=5, color=\"darkmagenta\"), name=norm_lbl),\n",
" go.Scatter(x=[solver(0.0)[0][0]], y=[solver(0.0)[0][1]], mode=\"markers\",\n",
" marker=dict(size=14, symbol=\"x\", color=\"teal\"), name=\"ŵ(λ)\"),\n",
" go.Scatter(x=[0, 0], y=[0, ymax], mode=\"lines\",\n",
" line=dict(width=2, dash=\"dot\", color=\"teal\"), showlegend=False)]\n",
" for tr, col in zip(dyn, [1, 1, 2]):\n",
" fig.add_trace(tr, row=1, col=col)\n",
" dyn_ix = list(range(len(fig.data) - 3, len(fig.data)))\n",
"\n",
" lams = np.linspace(0.0, 15.0, 16)\n",
" fig.frames = [go.Frame(name=f\"{l:.2f}\", traces=dyn_ix, data=[\n",
" go.Scatter(x=shape(l)[0], y=shape(l)[1]),\n",
" go.Scatter(x=[solver(l)[0][0]], y=[solver(l)[0][1]]),\n",
" go.Scatter(x=[l, l], y=[0, ymax])]) for l in lams]\n",
"\n",
" fig.update_layout(\n",
" template=\"plotly_white\", height=560,\n",
" sliders=[dict(active=0, pad=dict(l=100, t=55), steps=[\n",
" {\"label\": f\"λ = {l:.1f}\", \"method\": \"animate\",\n",
" \"args\": [[f\"{l:.2f}\"], {\"mode\": \"immediate\",\n",
" \"frame\": {\"duration\": 0, \"redraw\": True},\n",
" \"transition\": {\"duration\": 0}}]} for l in lams])],\n",
" updatemenus=[dict(type=\"buttons\", x=0.02, y=0, xanchor=\"left\", yanchor=\"bottom\",\n",
" direction=\"left\", buttons=[\n",
" dict(label=\"▶\", method=\"animate\", args=[None, {\"fromcurrent\": True,\n",
" \"frame\": {\"duration\": 600, \"redraw\": True}, \"transition\": {\"duration\": 50}}]),\n",
" dict(label=\"⏸\", method=\"animate\", args=[[None], {\"mode\": \"immediate\",\n",
" \"frame\": {\"duration\": 0, \"redraw\": False}}])])])\n",
" fig.update_xaxes(title_text=\"w₁\", range=[W0[0] - 10, W0[0] + 10], row=1, col=1)\n",
" fig.update_yaxes(title_text=\"w₂\", range=[W0[1] - 10, W0[1] + 10],\n",
" scaleanchor=\"x\", scaleratio=1, row=1, col=1)\n",
" fig.update_xaxes(title_text=\"λ\", row=1, col=2)\n",
" return fig\n",
"\n",
"regularization_figure(ridge_solution, \"ridge\").show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"regularization_figure(lasso_solution, \"lasso\").show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. **Least Squares $\\equiv$ Maximum Likelihood**\n",
"\n",
"Assume the target is the model output plus zero-mean Gaussian noise:\n",
"\n",
"$$\n",
"t = y(\\mathbf{x}, \\mathbf{w}) + \\varepsilon,\\qquad \\varepsilon \\sim \\mathcal{N}(0, \\sigma^2)\n",
"\\qquad\\Longrightarrow\\qquad\n",
"p(t \\mid \\mathbf{x}, \\mathbf{w}, \\sigma^2) = \\mathcal{N}\\big(t \\,\\big|\\, y(\\mathbf{x},\\mathbf{w}),\\, \\sigma^2\\big)\n",
"$$\n",
"\n",
"With i.i.d. data, the log-likelihood is\n",
"\n",
"$$\n",
"\\ln p(\\mathbf{t}\\mid\\mathbf{w},\\sigma^2)\n",
"= -\\frac{1}{2\\sigma^2}\\underbrace{\\sum_{n}\\big(t_n - \\mathbf{w}^\\top\\boldsymbol{\\phi}(\\mathbf{x}_n)\\big)^2}_{2\\,E(\\mathbf{w})}\n",
"\\;-\\;\\frac{N}{2}\\ln\\sigma^2 \\;-\\;\\frac{N}{2}\\ln 2\\pi\n",
"$$\n",
"\n",
"Only the first term depends on $\\mathbf{w}$, and it is $-E(\\mathbf{w})/\\sigma^2$. Maximizing the likelihood in $\\mathbf{w}$ is minimizing the sum-of-squares error, whatever $\\sigma^2$ is."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:51.531786Z",
"iopub.status.busy": "2026-09-15T03:39:51.531527Z",
"iopub.status.idle": "2026-09-15T03:39:51.543250Z",
"shell.execute_reply": "2026-09-15T03:39:51.542287Z"
}
},
"outputs": [],
"source": [
"rng2 = np.random.default_rng(0)\n",
"N, sigma_true = 80, 1.0\n",
"xm = np.linspace(0, 5, N)\n",
"Phi_m = np.column_stack([np.ones_like(xm), xm]) # [1, x] -> w = [intercept, slope]\n",
"w_true = np.array([1.0, 2.0])\n",
"tm = Phi_m @ w_true + rng2.normal(0, sigma_true, N)\n",
"\n",
"def sse(w): return float(np.sum((tm - Phi_m @ w) ** 2))\n",
"def loglik(w, s2): return -0.5 * sse(w) / s2 - 0.5 * N * np.log(s2) - 0.5 * N * np.log(2 * np.pi)\n",
"\n",
"w_ls = np.linalg.lstsq(Phi_m, tm, rcond=None)[0]\n",
"print(\"least-squares estimate w_hat :\", np.array2string(w_ls, precision=4))\n",
"print(\"true parameters w* :\", np.array2string(w_true, precision=4))\n",
"print(\"\\nSSE at w_hat :\", round(sse(w_ls), 4))\n",
"print(\"log-lik at w_hat :\", round(loglik(w_ls, sigma_true ** 2), 4))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:51.579274Z",
"iopub.status.busy": "2026-09-15T03:39:51.578997Z",
"iopub.status.idle": "2026-09-15T03:39:52.735153Z",
"shell.execute_reply": "2026-09-15T03:39:52.733955Z"
}
},
"outputs": [],
"source": [
"# Grid search the log-likelihood and confirm its argmax is the least-squares solution.\n",
"g0 = np.linspace(w_true[0] - 1.5, w_true[0] + 1.5, 301)\n",
"g1 = np.linspace(w_true[1] - 0.6, w_true[1] + 0.6, 301)\n",
"G0, G1 = np.meshgrid(g0, g1)\n",
"R = tm[None, None, :] - (G0[..., None] * Phi_m[:, 0] + G1[..., None] * Phi_m[:, 1])\n",
"SSE = np.sum(R ** 2, axis=-1)\n",
"LL = -0.5 * SSE / sigma_true**2 - 0.5 * N * np.log(sigma_true**2) - 0.5 * N * np.log(2 * np.pi)\n",
"\n",
"i_sse = np.unravel_index(SSE.argmin(), SSE.shape)\n",
"i_ll = np.unravel_index(LL.argmax(), LL.shape)\n",
"print(\"argmin SSE on the grid :\", (round(G0[i_sse], 4), round(G1[i_sse], 4)))\n",
"print(\"argmax logL on the grid:\", (round(G0[i_ll], 4), round(G1[i_ll], 4)))\n",
"print(\"same grid cell?\", i_sse == i_ll)\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(13, 5))\n",
"for ax, Z, name, cm in [(axs[0], SSE, \"$E(\\\\mathbf{w})$ (minimize)\", \"Blues_r\"),\n",
" (axs[1], LL, \"$\\\\ln p(\\\\mathbf{t}\\\\mid\\\\mathbf{w})$ (maximize)\", \"Blues\")]:\n",
" ax.contourf(G0, G1, Z, levels=40, cmap=cm)\n",
" ax.contour(G0, G1, Z, levels=18, colors=\"white\", linewidths=0.5, alpha=0.6)\n",
" ax.plot(*w_ls, marker=\"*\", ms=22, color=C_ALT, mec=\"k\", zorder=5)\n",
" ax.set(xlabel=\"$w_0$\", ylabel=\"$w_1$\", title=name)\n",
"plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 3.1 The noise variance comes out too\n",
"\n",
"Maximizing over $\\sigma^2$ as well gives the residual variance around the fitted function:\n",
"\n",
"$$\n",
"\\sigma^2_{\\text{ML}} = \\frac{1}{N}\\sum_{n=1}^{N}\\big(t_n - \\hat{\\mathbf{w}}^\\top\\boldsymbol{\\phi}(\\mathbf{x}_n)\\big)^2\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:52.737538Z",
"iopub.status.busy": "2026-09-15T03:39:52.737166Z",
"iopub.status.idle": "2026-09-15T03:39:52.892837Z",
"shell.execute_reply": "2026-09-15T03:39:52.891388Z"
}
},
"outputs": [],
"source": [
"sigma2_ml = np.mean((tm - Phi_m @ w_ls) ** 2)\n",
"print(f\"sigma^2_ML = {sigma2_ml:.4f}\")\n",
"print(f\"true sigma^2 = {sigma_true**2:.4f}\")\n",
"\n",
"grid_s2 = np.linspace(0.3, 2.5, 400)\n",
"ll_s2 = np.array([loglik(w_ls, s2) for s2 in grid_s2])\n",
"print(f\"argmax over the sigma^2 grid = {grid_s2[ll_s2.argmax()]:.4f}\")\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(grid_s2, ll_s2, lw=3, color=C_DATA)\n",
"ax.axvline(sigma2_ml, color=C_ALT, lw=2.5, ls=\"--\", label=r\"$\\sigma^2_{ML}$\")\n",
"ax.set(xlabel=r\"$\\sigma^2$\", ylabel=r\"$\\ln p(\\mathbf{t}\\mid\\hat{\\mathbf{w}},\\sigma^2)$\")\n",
"ax.legend(); plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. **Choosing a Noise Model**\n",
"\n",
"The chain runs *noise distribution $\\to$ likelihood $\\to$ error function*. Squared error is what a Gaussian gives you. A Laplace noise model, which has heavier tails, gives absolute error:\n",
"\n",
"$$\n",
"\\varepsilon \\sim \\mathcal{N}(0,\\sigma^2) \\;\\Rightarrow\\; -\\ln p \\propto \\varepsilon^2 \\;\\Rightarrow\\; \\textbf{MSE}\n",
"\\qquad\\qquad\n",
"\\varepsilon \\sim \\text{Laplace}(0,b) \\;\\Rightarrow\\; -\\ln p \\propto |\\varepsilon| \\;\\Rightarrow\\; \\textbf{MAE}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:52.895187Z",
"iopub.status.busy": "2026-09-15T03:39:52.894918Z",
"iopub.status.idle": "2026-09-15T03:39:53.284231Z",
"shell.execute_reply": "2026-09-15T03:39:53.282835Z"
}
},
"outputs": [],
"source": [
"eps = np.linspace(-4, 4, 400)\n",
"b = 1 / np.sqrt(2) # Laplace(b) has variance 2b^2, matching N(0,1)\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))\n",
"ax[0].plot(eps, st.norm(0, 1).pdf(eps), lw=3, color=C_DATA, label=\"Gaussian $\\\\mathcal{N}(0,1)$\")\n",
"ax[0].plot(eps, st.laplace(0, b).pdf(eps), lw=3, color=C_ALT, label=\"Laplace (same variance)\")\n",
"ax[0].set(title=\"Noise densities\", xlabel=r\"$\\varepsilon$\", ylabel=\"density\"); ax[0].legend()\n",
"\n",
"ax[1].plot(eps, 0.5 * eps ** 2, lw=3, color=C_DATA, label=r\"$-\\ln p_{\\rm Gauss} \\propto \\varepsilon^2/2$\")\n",
"ax[1].plot(eps, np.abs(eps) / b, lw=3, color=C_ALT, label=r\"$-\\ln p_{\\rm Laplace} = |\\varepsilon|/b$\")\n",
"ax[1].set(title=\"Implied error functions\", xlabel=r\"$\\varepsilon$\", ylabel=\"error\", ylim=(0, 8))\n",
"ax[1].legend()\n",
"plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 4.1 Behaviour under a single outlier\n",
"\n",
"The quadratic arm of squared error means one bad measurement can move the whole fit."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:53.286850Z",
"iopub.status.busy": "2026-09-15T03:39:53.286613Z",
"iopub.status.idle": "2026-09-15T03:39:53.658285Z",
"shell.execute_reply": "2026-09-15T03:39:53.657156Z"
}
},
"outputs": [],
"source": [
"x_o = np.linspace(0, 10, 40)\n",
"t_o = 1 + 2 * x_o + np.random.default_rng(3).normal(0, 1.0, 40)\n",
"t_o_bad = t_o.copy(); t_o_bad[35] += 40.0 # one corrupted measurement\n",
"\n",
"Phi_o = np.column_stack([np.ones_like(x_o), x_o])\n",
"def fit_mse(t): return np.linalg.lstsq(Phi_o, t, rcond=None)[0]\n",
"def fit_mae(t): return minimize(lambda w: np.abs(t - Phi_o @ w).sum(),\n",
" fit_mse(t), method=\"Nelder-Mead\",\n",
" options=dict(xatol=1e-9, fatol=1e-9, maxiter=8000)).x\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(13, 5), sharey=True)\n",
"for ax, tt, name in [(axs[0], t_o, \"clean data\"),\n",
" (axs[1], t_o_bad, \"one corrupted measurement\")]:\n",
" ax.scatter(x_o, tt, s=55, color=C_DATA, alpha=0.75, zorder=3)\n",
" wm, wa = fit_mse(tt), fit_mae(tt)\n",
" ax.plot(x_o, Phi_o @ wm, lw=3, color=C_RESID, label=f\"MSE fit: slope {wm[1]:.2f}\")\n",
" ax.plot(x_o, Phi_o @ wa, lw=3, ls=\"--\", color=C_SPAN, label=f\"MAE fit: slope {wa[1]:.2f}\")\n",
" ax.set(xlabel=\"$x$\", title=name); ax.legend(loc=\"upper left\", fontsize=11)\n",
"axs[0].set_ylabel(\"$t$\")\n",
"plt.tight_layout(); plt.show()\n",
"\n",
"print(\"slope, clean data -> MSE %.3f | MAE %.3f\" % (fit_mse(t_o)[1], fit_mae(t_o)[1]))\n",
"print(\"slope, with outlier -> MSE %.3f | MAE %.3f\" % (fit_mse(t_o_bad)[1], fit_mae(t_o_bad)[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5. **Regularization as a Prior Belief**\n",
"\n",
"Lecture 06 introduced ridge as \"add a penalty so the weights stay small\". The probabilistic view says where that penalty comes from: a prior belief that the weights are small. Put $\\mathbf{w} \\sim \\mathcal{N}(\\mathbf{0}, \\tau^2 I)$ and maximize the posterior instead of the likelihood:\n",
"\n",
"$$\n",
"\\hat{\\mathbf{w}}_{\\text{MAP}} = \\arg\\max_{\\mathbf{w}} \\big[\\ln p(\\mathbf{t}\\mid\\mathbf{w},\\sigma^2) + \\ln p(\\mathbf{w})\\big]\n",
"= \\arg\\min_{\\mathbf{w}} \\left[\\tfrac{1}{2\\sigma^2}\\lVert\\mathbf{t}-\\Phi\\mathbf{w}\\rVert^2 + \\tfrac{1}{2\\tau^2}\\lVert\\mathbf{w}\\rVert^2\\right]\n",
"$$\n",
"\n",
"Multiply through by $\\sigma^2$ and this is ridge regression with\n",
"\n",
"$$\n",
"\\lambda = \\frac{\\sigma^2}{\\tau^2}\n",
"$$\n",
"\n",
"$\\lambda$ is a ratio of variances: how noisy the data is believed to be, relative to how large the weights are believed to be."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:53.660814Z",
"iopub.status.busy": "2026-09-15T03:39:53.660586Z",
"iopub.status.idle": "2026-09-15T03:39:53.670802Z",
"shell.execute_reply": "2026-09-15T03:39:53.669895Z"
}
},
"outputs": [],
"source": [
"sigma2, tau2 = 0.25, 0.10\n",
"lam_map = sigma2 / tau2\n",
"\n",
"raw = build_Phi(x_s, 7, bias=False)\n",
"Phi_map = (raw - raw.mean(0)) / raw.std(0)\n",
"t_c = t_s - t_s.mean() # centred, so no intercept to handle\n",
"\n",
"w_map = np.linalg.solve(Phi_map.T @ Phi_map / sigma2 + np.eye(Phi_map.shape[1]) / tau2,\n",
" Phi_map.T @ t_c / sigma2)\n",
"w_ridge = np.linalg.solve(Phi_map.T @ Phi_map + lam_map * np.eye(Phi_map.shape[1]),\n",
" Phi_map.T @ t_c)\n",
"\n",
"print(f\"lambda = sigma^2 / tau^2 = {sigma2} / {tau2} = {lam_map:.1f}\\n\")\n",
"print(\"MAP (posterior mode):\", np.array2string(w_map, precision=6))\n",
"print(\"Ridge (lambda = 2.5) :\", np.array2string(w_ridge, precision=6))\n",
"print(\"\\nmax |difference| =\", f\"{np.abs(w_map - w_ridge).max():.2e}\")\n",
"\n",
"print(\"\\nA tighter prior shrinks harder:\")\n",
"for tau2_i in [10.0, 1.0, 0.1, 0.01]:\n",
" li = sigma2 / tau2_i\n",
" wi = np.linalg.solve(Phi_map.T @ Phi_map + li * np.eye(Phi_map.shape[1]), Phi_map.T @ t_c)\n",
" print(f\" tau^2 = {tau2_i:>6} -> lambda = {li:>7.2f} -> ||w|| = {np.linalg.norm(wi):.3f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The same substitution with a **Laplace** prior on $\\mathbf{w}$ gives the lasso. The two regularizers we used mechanically last lecture are two choices of prior."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 6. **The Bias-Variance Decomposition**\n",
"\n",
"Take the expected squared error at an input $x$, averaged over training sets $\\mathcal{D}$:\n",
"\n",
"$$\n",
"\\mathbb{E}_{\\mathcal{D}}\\big[(t - y(x;\\mathcal{D}))^2\\big]\n",
"= \\underbrace{\\sigma^2}_{\\text{noise}}\n",
"+ \\underbrace{\\big(\\mathbb{E}_{\\mathcal{D}}[y(x;\\mathcal{D})] - f(x)\\big)^2}_{\\text{bias}^2}\n",
"+ \\underbrace{\\mathbb{E}_{\\mathcal{D}}\\big[(y(x;\\mathcal{D}) - \\mathbb{E}_{\\mathcal{D}}[y(x;\\mathcal{D})])^2\\big]}_{\\text{variance}}\n",
"$$\n",
"\n",
"- **Noise** is irreducible. No model removes it.\n",
"- **Bias** is how far the *average* fitted function sits from the truth: underfitting.\n",
"- **Variance** is how much the fitted function moves when the training set changes: overfitting.\n",
"\n",
"Every term is an expectation over training sets, so we can estimate all three by simulation: draw many training sets, fit each, and look at the spread of the fitted functions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:53.672671Z",
"iopub.status.busy": "2026-09-15T03:39:53.672458Z",
"iopub.status.idle": "2026-09-15T03:39:54.457278Z",
"shell.execute_reply": "2026-09-15T03:39:54.455931Z"
}
},
"outputs": [],
"source": [
"SIGMA_BV, NSETS = 0.30, 400\n",
"f_true = lambda x: np.sin(5 * x)\n",
"xg_bv = np.linspace(-0.9, 0.9, 120)\n",
"\n",
"def Phi_leg(x, deg):\n",
" \"Legendre basis: well conditioned, so what we measure is variance and not rounding error.\"\n",
" return L.legvander(np.asarray(x).ravel(), deg)\n",
"\n",
"def decompose(n_train, degrees, seed=189):\n",
" rng = np.random.default_rng(seed)\n",
" out = []\n",
" for d in degrees:\n",
" preds = np.empty((NSETS, len(xg_bv)))\n",
" for s in range(NSETS):\n",
" xt = rng.uniform(-0.9, 0.9, n_train)\n",
" tt = f_true(xt) + SIGMA_BV * rng.standard_normal(n_train)\n",
" w = np.linalg.lstsq(Phi_leg(xt, d), tt, rcond=None)[0]\n",
" preds[s] = Phi_leg(xg_bv, d) @ w\n",
" bias2 = float(np.mean((preds.mean(0) - f_true(xg_bv)) ** 2))\n",
" var = float(np.mean(preds.var(0)))\n",
" out.append((d, bias2, var, SIGMA_BV ** 2, bias2 + var + SIGMA_BV ** 2, preds))\n",
" return out\n",
"\n",
"degrees = list(range(0, 11))\n",
"res100 = decompose(100, degrees)\n",
"\n",
"print(f\"N_train = 100, {NSETS} training sets per degree\\n\")\n",
"print(f\"{'deg':>4} {'bias^2':>9} {'variance':>9} {'noise':>8} {'total':>9}\")\n",
"for d, b2, v, nz, tot, _ in res100:\n",
" print(f\"{d:>4} {b2:>9.4f} {v:>9.4f} {nz:>8.4f} {tot:>9.4f}\")\n",
"best = min(res100, key=lambda r: r[4])\n",
"print(f\"\\nlowest expected error at degree {best[0]} (total {best[4]:.4f})\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 6.1 Seeing bias and variance directly\n",
"\n",
"Plot the fitted function from 40 of the training sets, plus their average, against the truth. Bias is the gap between the average (dashed) and the truth (solid). Variance is the width of the grey band."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:54.459627Z",
"iopub.status.busy": "2026-09-15T03:39:54.459407Z",
"iopub.status.idle": "2026-09-15T03:39:55.079996Z",
"shell.execute_reply": "2026-09-15T03:39:55.077228Z"
}
},
"outputs": [],
"source": [
"show_degrees = [1, 3, 5, 9]\n",
"fig, axs = plt.subplots(1, 4, figsize=(17, 4.2), sharey=True)\n",
"for ax, d in zip(axs, show_degrees):\n",
" rec = [r for r in res100 if r[0] == d][0]\n",
" preds = rec[5]\n",
" for p in preds[:40]:\n",
" ax.plot(xg_bv, p, color=\"0.55\", lw=0.8, alpha=0.55)\n",
" ax.plot(xg_bv, f_true(xg_bv), lw=3, color=C_DATA, label=\"truth $f(x)$\")\n",
" ax.plot(xg_bv, preds.mean(0), lw=3, ls=\"--\", color=C_ALT, label=\"average fit\")\n",
" ax.set_title(f\"degree {d}\\nbias$^2$ = {rec[1]:.3f}\\nvariance = {rec[2]:.3f}\", fontsize=13)\n",
" ax.set(xlabel=\"$x$\", ylim=(-2.0, 2.0))\n",
"axs[0].set_ylabel(\"$y$\"); axs[0].legend(fontsize=10, loc=\"lower left\")\n",
"plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:55.082329Z",
"iopub.status.busy": "2026-09-15T03:39:55.082106Z",
"iopub.status.idle": "2026-09-15T03:39:55.409514Z",
"shell.execute_reply": "2026-09-15T03:39:55.408371Z"
}
},
"outputs": [],
"source": [
"b2 = np.array([r[1] for r in res100])\n",
"vr = np.array([r[2] for r in res100])\n",
"nz = np.array([r[3] for r in res100])\n",
"tot = b2 + vr + nz\n",
"\n",
"fig, ax = plt.subplots(figsize=(8.5, 5))\n",
"ax.plot(degrees, b2, \"o-\", lw=3, ms=8, color=C_DATA, label=\"bias$^2$\")\n",
"ax.plot(degrees, vr, \"s-\", lw=3, ms=8, color=C_RESID, label=\"variance\")\n",
"ax.plot(degrees, nz, \":\", lw=2.5, color=\"0.4\", label=\"noise $\\\\sigma^2$\")\n",
"ax.plot(degrees, tot, \"^-\", lw=3, ms=8, color=C_SPAN, label=\"total\")\n",
"ax.axvline(degrees[int(tot.argmin())], color=C_ALT, ls=\"--\", lw=2.5,\n",
" label=f\"minimum at degree {degrees[int(tot.argmin())]}\")\n",
"ax.set(xlabel=\"polynomial degree\", ylabel=\"expected squared error\", yscale=\"log\",\n",
" xticks=degrees)\n",
"ax.legend(fontsize=11); plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 6.2 More data supports more complexity\n",
"\n",
"Nothing about the tradeoff is a property of the model alone. Shrink the training set and the variance term grows at every degree, so the minimum moves left: with less data the best model is *simpler* than the truth."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:55.411716Z",
"iopub.status.busy": "2026-09-15T03:39:55.411499Z",
"iopub.status.idle": "2026-09-15T03:39:56.442961Z",
"shell.execute_reply": "2026-09-15T03:39:56.441827Z"
}
},
"outputs": [],
"source": [
"res25 = decompose(25, degrees)\n",
"tot25 = np.array([r[4] for r in res25])\n",
"\n",
"fig, ax = plt.subplots(figsize=(8.5, 5))\n",
"for res, tot_i, col, name in [(res25, tot25, C_RESID, \"N = 25\"),\n",
" (res100, tot, C_DATA, \"N = 100\")]:\n",
" ax.plot(degrees, tot_i, \"o-\", lw=3, ms=8, color=col, label=f\"{name}, total error\")\n",
" k = int(np.argmin(tot_i))\n",
" ax.plot(degrees[k], tot_i[k], \"*\", ms=22, color=C_ALT, mec=\"k\", zorder=5)\n",
"ax.set(xlabel=\"polynomial degree\", ylabel=\"expected squared error\", yscale=\"log\",\n",
" xticks=degrees)\n",
"ax.legend(fontsize=12); plt.tight_layout(); plt.show()\n",
"\n",
"print(f\"N = 25 -> minimum at degree {degrees[int(tot25.argmin())]} (total {tot25.min():.4f})\")\n",
"print(f\"N = 100 -> minimum at degree {degrees[int(tot.argmin())]} (total {tot.min():.4f})\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 6.3 Regularization moves along the same curve\n",
"\n",
"Model complexity is not only the number of basis functions. Hold the degree fixed at 9 and turn $\\lambda$ instead: large $\\lambda$ buys bias and sells variance, exactly as reducing the degree does."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-15T03:39:56.445278Z",
"iopub.status.busy": "2026-09-15T03:39:56.445027Z",
"iopub.status.idle": "2026-09-15T03:39:58.965904Z",
"shell.execute_reply": "2026-09-15T03:39:58.964353Z"
}
},
"outputs": [],
"source": [
"def decompose_ridge(n_train, deg, lambdas, seed=189):\n",
" rng = np.random.default_rng(seed)\n",
" out = []\n",
" for lam in lambdas:\n",
" preds = np.empty((NSETS, len(xg_bv)))\n",
" for s in range(NSETS):\n",
" xt = rng.uniform(-0.9, 0.9, n_train)\n",
" tt = f_true(xt) + SIGMA_BV * rng.standard_normal(n_train)\n",
" A = Phi_leg(xt, deg)\n",
" P = np.eye(deg + 1); P[0, 0] = 0.0 # intercept left unpenalized\n",
" w = np.linalg.solve(A.T @ A + lam * P, A.T @ tt)\n",
" preds[s] = Phi_leg(xg_bv, deg) @ w\n",
" bias2 = float(np.mean((preds.mean(0) - f_true(xg_bv)) ** 2))\n",
" var = float(np.mean(preds.var(0)))\n",
" out.append((lam, bias2, var, bias2 + var + SIGMA_BV ** 2))\n",
" return out\n",
"\n",
"lams = np.logspace(-4, 3, 22)\n",
"rr = decompose_ridge(100, 9, lams)\n",
"b2r = np.array([r[1] for r in rr]); vrr = np.array([r[2] for r in rr])\n",
"totr = np.array([r[3] for r in rr])\n",
"\n",
"fig, ax = plt.subplots(figsize=(8.5, 5))\n",
"ax.plot(lams, b2r, \"o-\", lw=3, ms=7, color=C_DATA, label=\"bias$^2$\")\n",
"ax.plot(lams, vrr, \"s-\", lw=3, ms=7, color=C_RESID, label=\"variance\")\n",
"ax.plot(lams, totr, \"^-\", lw=3, ms=7, color=C_SPAN, label=\"total\")\n",
"ax.axvline(lams[int(totr.argmin())], color=C_ALT, ls=\"--\", lw=2.5,\n",
" label=rf\"best $\\lambda$ = {lams[int(totr.argmin())]:.3g}\")\n",
"ax.set(xscale=\"log\", yscale=\"log\", xlabel=r\"$\\lambda$ (degree fixed at 9)\",\n",
" ylabel=\"expected squared error\")\n",
"ax.legend(fontsize=11); plt.tight_layout(); plt.show()\n",
"\n",
"print(f\"lambda -> 0 : bias^2 {b2r[0]:.4f}, variance {vrr[0]:.4f}\")\n",
"print(f\"best lambda : bias^2 {b2r[totr.argmin()]:.4f}, variance {vrr[totr.argmin()]:.4f}\"\n",
" f\" (lambda = {lams[int(totr.argmin())]:.3g})\")\n",
"print(f\"lambda -> large : bias^2 {b2r[-1]:.4f}, variance {vrr[-1]:.4f}\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}