--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: Generalized Linear Models :::{admonition} Get this problem set :class: seealso {download}`Download everything (Topic3.2-Generalized_Linear_Models_Problems.zip) ` — the notebook and `flammability_calibration.csv`, in a folder that is ready to run as-is. The Download badge at the top of the page will also give you the notebook on its own. On Vocareum everything is already set up for you. ::: :::{admonition} Before you start :class: tip This problem set accompanies {doc}`Generalized Linear Models `. It is worth **100 points**: - **Part A — Skill Checks (30 pts)** — short answers, auto-graded, resubmit as often as you like until they pass. - **Part B — Visualization (35 pts)** — plots plus written interpretation, peer graded. - **Part C — Open Ended (35 pts)** — one synthesis problem, peer graded. The parts build on each other: Part A works out the syntax you need for Part B, and Part B produces the evidence you argue from in Part C. Do them in order. ::: ## Setup This set uses the methane ignition tests from {doc}`the Topic 3.1 problem set `, plus a second campaign. **Campaign A** is the same 800 tests as before, all run at a controlled 25 °C: a mixture of methane, added nitrogen and air is sparked, and `ignited` records whether a flame propagated. **Campaign B** is 600 further tests run on the same rig at temperatures from 25 to 150 °C. Heating a mixture widens the flammable range at both ends and lowers the oxygen concentration a flame can survive on, so campaign B's flammable region is not campaign A's. **Campaign B is not for Parts A and B.** Everything up to Part C uses campaign A only. Topic 3.1 established what the region looks like: bounded below in fuel concentration, bounded above, and closed off again once enough nitrogen has been added. A detector built there out of three separate threshold conditions reached an F1 of about 0.83. This set asks a different question — what does **one fitted model** do with the same region — and the answer turns out to depend far more on which features you give it than on which loss function you minimize. Three things to know about the fitting below. 1. **Labels come in two codings.** The loss functions want $y \in \{-1, +1\}$; the chapter's `acc_prec_recall` wants $\{0, 1\}$. Both are defined in the Setup cell. 2. **The features are standardized**, using means and standard deviations computed from the **training rows only**. Squared features on raw volume percents would range over thousands and the optimizer would struggle. 3. **Every fit starts from `w = 0`** and uses `minimize` with its default settings, so your numbers will match the reference values. ```{code-cell} ipython3 %matplotlib inline import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import minimize from sklearn.model_selection import train_test_split try: plt.style.use('../settings/plot_style.mplstyle') # available inside the book except OSError: pass # downloaded notebook: use defaults df = pd.read_csv('data/flammability_calibration.csv') A = df[df['campaign'] == 'A'].reset_index(drop=True) fuel = A['fuel_vol_pct'].values n2 = A['n2_added_vol_pct'].values y01 = A['ignited'].values # {0, 1} — for acc_prec_recall y = y01 * 2 - 1 # {-1, +1} — for the loss functions idx = np.arange(len(A)) train_idx, test_idx = train_test_split(idx, test_size=0.3, random_state=0, stratify=y01) # Standardizers from the TRAINING rows only. mu_f, sd_f = fuel[train_idx].mean(), fuel[train_idx].std() mu_n, sd_n = n2[train_idx].mean(), n2[train_idx].std() z1 = (fuel - mu_f) / sd_f z2 = (n2 - mu_n) / sd_n X_lin = np.column_stack([z1, z2]) # the two raw features X_quad = np.column_stack([z1, z2, z1**2, z2**2, z1 * z2]) # full quadratic basis print(f"campaign A: {len(A)} tests campaign B: {int((df['campaign'] == 'B').sum())} tests") print(f"{len(train_idx)} training rows, {len(test_idx)} test rows " f"({int(y01[test_idx].sum())} of the test rows ignited)") print(f"X_lin {X_lin.shape} X_quad {X_quad.shape}") ``` These are the chapter's functions, reproduced so the notebook runs on its own. ```{code-cell} ipython3 def add_intercept(X): """Prepend a column of ones to X for the bias term.""" return np.append(np.ones((X.shape[0], 1)), X, axis=1) def linear_classifier(X, w): """Return a boolean array: True where the model predicts class 1.""" return np.dot(add_intercept(X), w) > 0 def max_cost(w, X, y): Xb = np.dot(add_intercept(X), w) return np.sum(np.maximum(0, -y * Xb)) def softmax_cost(w, X, y): Xb = np.dot(add_intercept(X), w) return np.sum(np.log(1 + np.exp(-y * Xb))) def margin_cost(w, X, y): Xb = np.dot(add_intercept(X), w) return np.sum(np.maximum(0, 1 - y * Xb)) def margin_cost_squared(w, X, y): Xb = np.dot(add_intercept(X), w) return np.sum(np.maximum(0, 1 - y * Xb) ** 2) def acc_prec_recall(y_model, y_actual): TP = np.sum(np.logical_and(y_model == y_actual, y_model == 1)) TN = np.sum(np.logical_and(y_model == y_actual, y_model == 0)) FP = np.sum(np.logical_and(y_model != y_actual, y_model == 1)) FN = np.sum(np.logical_and(y_model != y_actual, y_model == 0)) acc = (TP + TN) / (TP + TN + FP + FN) if TP == 0: prec, recall = 0, 0 else: prec = TP / (TP + FP) recall = TP / (TP + FN) return acc, prec, recall def f1_score(y_model, y_actual): """F1 for the positive class, built from the chapter's acc_prec_recall.""" _, p, r = acc_prec_recall(y_model, y_actual) return 0.0 if (p + r) == 0 else 2 * p * r / (p + r) ``` Two more helpers, specific to this set. `basis_lin` and `basis_quad` rebuild the two design matrices from raw volume percents, and `plot_boundary` draws a fitted model's decision boundary over the mixture plane so you do not have to do the contour bookkeeping by hand. ```{code-cell} ipython3 def basis_lin(f_pct, n_pct): """The two standardized features, built from raw volume percents.""" return np.column_stack([(f_pct - mu_f) / sd_f, (n_pct - mu_n) / sd_n]) def basis_quad(f_pct, n_pct): """The full quadratic basis, built from raw volume percents.""" z1 = (f_pct - mu_f) / sd_f z2 = (n_pct - mu_n) / sd_n return np.column_stack([z1, z2, z1**2, z2**2, z1 * z2]) def plot_boundary(ax, w, basis, f_max=20.0, n_max=50.0, n_grid=400, **kwargs): """Draw the decision boundary of weight vector `w` on the axes `ax`. `basis` is basis_lin or basis_quad, and must be the one `w` was fitted on. Returns True if a boundary was drawn. A model whose decision value never changes sign over the plane has no boundary to draw, and returns False. """ fg, ng = np.meshgrid(np.linspace(0, f_max, n_grid), np.linspace(0, n_max, n_grid)) dec = np.dot(add_intercept(basis(fg.ravel(), ng.ravel())), w).reshape(fg.shape) if dec.min() >= 0 or dec.max() <= 0: return False # no sign change: nothing to contour ax.contour(fg, ng, dec, levels=[0], **kwargs) return True # these rebuild exactly the matrices the Setup already made assert np.allclose(basis_lin(fuel, n2), X_lin) assert np.allclose(basis_quad(fuel, n2), X_quad) print("basis helpers agree with X_lin and X_quad") ``` Every question below reports **F1 for the flammable class** rather than accuracy. Only about a quarter of the tests ignite, so a model that never predicts ignition already scores 0.746 on accuracy, and several of the fits below do exactly that. F1 makes that failure visible. Run this once. It defines the `check` helper used after each Part A answer. It will not run inside the book — use the downloaded notebook. ```{code-cell} ipython3 :tags: [skip-execution] # Run this once. It defines `check`, used after each Part A answer below. import hashlib _EXPECTED = { 'q1': (4, ['855bd67d7574', '579811d583d7', 'abe9f752f166', '3e575a0e767f', 'e92857b87276', '8226ad2a43f9', '4861a0075db2']), 'q2': (4, ['4d7d761974e0', 'e3286f5874ed', '032948267672', '74cab0799db7', '5850cc3d7757', '9b1963da5d3d', '4cd6b956fd94', '585b2fedf0d6', '9a30dbe74781', '30527ef1acca', '4533fd474a64', 'fb15d8fea2bc', 'fe4c5a83d8eb', '38eedeb4a701', '27c3956911ff', '8a750d9bcd6d', 'f60dd8d8e6ea', 'dc5dc218dd5a', '131997e8a7a2', 'c2b8dfe1b52b', '2358c4d84b02']), 'q3': (4, ['eb7d66cb8e1c', '650e6b3ebb7f', 'ec3cf0aa426c', '61fa80d72ff5', 'bf1baeabe4f9', '47acb56e26ec', '1165ebacb6fc', '6f31900e731a', 'e9c740e7bf61', '0ae9a8b1ebb4', '58802acf83cb', '4ff0748233b0', '370b856c08c0', 'f8c49fcfb0e2', '5cb5c11ca76c', 'fa4724762cbc', 'ab05e8b66739', '778709960dde', '51086d7bf7b0', '6fe9a5747a19', '58c062ffa8d8']), } def check(qid, value): """Compare an answer against a stored digest. Practice only -- nothing is recorded.""" places, digests = _EXPECTED[qid] if value is Ellipsis: print(f"{qid}: not answered yet") return got = hashlib.sha256(f"{float(value):.{places}f}".encode()).hexdigest()[:12] print(f"{qid}: {'correct' if got in digests else 'not the expected value'}") ``` --- ## Part A — Skill Checks (30 pts) Three questions, 10 points each, and together they are one controlled experiment. A1 and A2 minimize the **same loss** on **different features**. A2 and A3 minimize **different losses** on the **same features**. Compare the sizes of the two effects when you have all three. Fit on `train_idx`, score on `test_idx`, start every fit from `w = 0`. ### A1. A linear model on the raw features (10 pts) :::{exercise} :label: pr-cls-linear-ceiling Minimize `margin_cost_squared` over `X_lin[train_idx]` with labels `y[train_idx]`, starting from a zero weight vector of the right length. Predict on `X_lin[test_idx]` with `linear_classifier`, convert the boolean predictions to integers, and score them against `y01[test_idx]` with `f1_score`. Assign the F1 score to `f1_linear`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE f1_linear = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q1", f1_linear) ``` ### A2. The same loss, on a quadratic basis (10 pts) :::{exercise} :label: pr-cls-quadratic-basis Repeat A1 exactly — same loss, same split, same starting point — but fit `X_quad` instead of `X_lin`. `X_quad` adds $z_1^2$, $z_2^2$ and $z_1 z_2$ to the two standardized features, so the model can express a boundary that curves. Assign the resulting F1 score to `f1_quadratic`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE f1_quadratic = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q2", f1_quadratic) ``` ### A3. A different loss, on the same basis (10 pts) :::{exercise} :label: pr-cls-logistic-quad Now change the loss and leave the features alone: minimize `softmax_cost` — the logistic-regression objective the chapter derives by smoothing the perceptron loss — over `X_quad[train_idx]`, again from a zero start. Assign the resulting F1 score to `f1_logistic`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE f1_logistic = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q3", f1_logistic) ``` --- ## Part B — Visualization (35 pts) Part A gave three numbers. This part shows what is behind the first of them, and what the third one buys. ::::{exercise} :label: pr-cls-loss-ceiling **1. Four losses, one ceiling.** Fit all four of the chapter's loss functions — `max_cost`, `softmax_cost`, `margin_cost` and `margin_cost_squared` — on `X_lin[train_idx]`, each from a zero start. Report a table of the fitted weights, the final loss value, the test accuracy and the test F1 for each. Plot the four decision boundaries on one scatter of campaign A in the (`fuel_vol_pct`, `n2_added_vol_pct`) plane, colored by `ignited`. Use the `plot_boundary` helper from the Setup with `basis_lin`; it returns `False` for a model that has no boundary to draw, so check the return value and say in your legend which losses produced one. In **one or two sentences**: say what the four have in common, and how the best of them compares to a model that never predicts ignition. **2. Why two of them report nothing.** Two of the four score F1 exactly 0. Evaluate `max_cost` and `softmax_cost` at $w = [0, 0, 0]$ — just evaluate them, do not fit anything — and report both values. Now walk toward the origin and watch. Take three starting points, say $w_a = [-10, -4, -10]$, $w_b = [0, 1, 1]$ and $w_c = [1, -1, 0.5]$, and for $t$ from 0 to 1.5 evaluate each loss at $t\,w$. Plot the three `max_cost` curves against $t$ on one panel and the three `softmax_cost` curves on a second, marking $t = 0$ on both. Finally, minimize `max_cost` from each of the three starts and report $\lVert w \rVert$ at each solution. In **two sentences**: say what the optimizer found, and why that is a property of the loss function rather than a failure of the optimizer. Your `max_cost` panel should make the reason visible. **3. What the quadratic basis recovered.** Take the `X_quad` fit from A2 and plot its decision boundary over the same scatter with `plot_boundary` — this time passing `basis_quad`, since the helper has to build the same features the weights were fitted on. Then read the boundary off numerically: at added nitrogen of 0, 15 and 30 vol %, find the lowest and highest fuel concentration your model calls flammable. Compare each against the published methane limits — a lower limit of **5.0 vol %**, and an upper limit that falls from **15.8 vol %** as nitrogen is added, reaching the lean limit at about 41 vol % ([Zlochower and Green, 2009](https://doi.org/10.1016/j.jlp.2009.03.006)). The published upper limit at a given added nitrogen $n$ is $$U(n) = 15.8 + (5.55 - 15.8)\,\frac{n}{41.47}$$ Report the six recovered numbers, the three published pairs, and the deviations. In **two sentences**: say how close the recovered lean limit is to the published value, and where in the region the model is least accurate. :::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE fig, ax = plt.subplots(figsize=(7, 5)) ``` --- ## Part C — Open Ended (35 pts) Part B ended with a boundary that follows the published limits closely. This part asks what regularization does to it — and then whether any of it survives contact with a campaign run at temperatures the model never saw. The chapter's support-vector machine adds an $L_2$ penalty on the weights, excluding the intercept, to the margin loss. Use this version, which penalizes the **squared** margin so that the sweep below is reproducible whichever optimizer you use: ``` def regularized_cost_squared(w, X, y, alpha=1.0): Xb = np.dot(add_intercept(X), w) cost = np.sum(np.maximum(0, 1 - y * Xb) ** 2) cost += alpha * np.linalg.norm(w[1:], 2) # intercept excluded return cost ``` Two quantities the chapter defines, both computed from a fitted $w$: - the **margin width**, $2 / \lVert w_{1:} \rVert_2$, excluding the intercept - the **support vectors**, the training points with $y_i (\mathbf{x}_i \cdot w) \le 1$ ::::{exercise} :label: pr-cls-margin-alpha **1. Sweep the penalty.** Fit `regularized_cost_squared` on `X_quad[train_idx]` from a zero start for $\alpha \in \{0, 0.01, 0.1, 1, 10, 100\}$. Report a table of $\alpha$, the margin width, the number of support vectors, and the campaign-A test F1. Then plot it: the support-vector count against $\alpha$ on a logarithmic $x$ axis, and on a second panel the campaign-A test F1 against $\alpha$ on the same axis. (Place $\alpha = 0$ at the left edge or omit it and say which you did.) In **two sentences**: say whether the margin and the test F1 move together as $\alpha$ grows, relate the support-vector count to the margin width using the chapter's ridge-regression analogy, and say what the model has become by $\alpha = 100$. **2. Open campaign B.** Now load the 600 campaign-B tests. Standardize them with the **campaign-A training** means and standard deviations — not their own — and build the same quadratic basis. Score every $\alpha$ from task 1 on campaign B. Split campaign B into three temperature bins — 25–50 °C, 50–100 °C and above 100 °C — report the number of tests in each, and report F1 in each bin for every $\alpha$. In **two sentences**: say which bin is hardest and whether the ordering across the three bins is clean, and say whether any $\alpha$ improved campaign-B performance over $\alpha = 0$. **3. Choose.** Pick one $\alpha$ to deploy. State the margin, the support-vector count, and the campaign-B precision and recall it gives you. Defend the choice in **two sentences** from your own table — including, if that is what your numbers show, the case for not regularizing at all. Then name one change to how campaign A was run that would make the model transfer better to campaign B. :::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ```