--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: Linear Regression :::{admonition} Get this problem set :class: seealso {download}`Download everything (Topic1.3-Linear_Regression_Problems.zip) ` — the notebook and `methanol_IR.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}`Linear Regression `. 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 The chapter fit basis functions to an **ethanol** infrared spectrum. Here you will work with a **methanol** spectrum measured under different conditions: gas phase, 70 mmHg of methanol in nitrogen to 600 mmHg total, 5 cm path length, 2 cm⁻¹ resolution. It comes from the Coblentz Society collection in the NIST Chemistry WebBook (spectrum 8791, measured at Dow Chemical in 1964), converted from transmittance to absorbance by $A = -\log_{10} T$. Methanol is the simplest alcohol, so its spectrum is a stripped-down version of ethanol's: a strong C–O stretch near 1030 cm⁻¹ and a cluster of C–H stretches between 2800 and 3050 cm⁻¹, without ethanol's extra methylene modes. That makes it a good place to ask how many basis functions a band actually needs. ```{code-cell} ipython3 %matplotlib inline import numpy as np import pandas as pd import matplotlib.pyplot as plt try: plt.style.use('../settings/plot_style.mplstyle') # available inside the book except OSError: pass # downloaded notebook: use defaults df = pd.read_csv('data/methanol_IR.csv') x_all = df['wavenumber [cm^-1]'].values y_all = df['absorbance'].values print(f"{len(x_all)} points, {x_all.min():.1f}-{x_all.max():.1f} cm^-1") ``` ```{code-cell} ipython3 fig, ax = plt.subplots(figsize=(8, 3)) ax.plot(x_all, y_all, lw=0.8) ax.set_xlabel('wavenumber [cm$^{-1}$]') ax.set_ylabel('absorbance'); ``` Run this once to load the autograder. It will not run inside the book — use the downloaded notebook. ```{code-cell} ipython3 :tags: [skip-execution] import otter grader = otter.Notebook() ``` --- ## Part A — Skill Checks (30 pts) Three questions, 10 points each. Each asks you to assign a single number to a named variable. Run the check cell after each one; there is no limit on attempts. ### A1. Locate the C–O stretch (10 pts) :::{exercise} :label: pr-nm-meoh-co-peak Restrict the spectrum to the window $950 \le \tilde\nu \le 1120$ cm⁻¹, which is the region every part of this problem set works in. Find the wavenumber at which absorbance is largest in that window. Two NumPy tools you need here. A **boolean array** built by comparing an array against a value can be used to index it, so `band = (x_all >= 950) & (x_all <= 1120)` and then `x_all[band]` selects the window. And `np.argmax` returns the *index* of the largest element, which you can use to look up the corresponding wavenumber. Assign the wavenumber in cm⁻¹ to `peak_wavenumber`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE peak_wavenumber = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] grader.check("q1") ``` ### A2. Fit one Gaussian to the band (10 pts) :::{exercise} :label: pr-nm-meoh-gauss-fit Model the same window as a single Gaussian sitting on a flat baseline: $$ A(\tilde\nu) \;=\; w_0 \exp\!\left[-\frac{(\tilde\nu - c)^2}{2\sigma^2}\right] + w_1 $$ with the center $c$ fixed at your answer from A1 and $\sigma = 25$ cm⁻¹ (the same width the chapter used). Build the two-column design matrix $\bar{\bar{X}}$ — one column for the Gaussian, one of ones for the baseline — and solve the normal equations $\bar{\bar{X}}^T\bar{\bar{X}}\vec{w} = \bar{\bar{X}}^T\vec{y}$ with `np.linalg.solve`. Assign the Gaussian coefficient $w_0$ to `w_gauss`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE w_gauss = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] grader.check("q2") ``` ### A3. Score the fit (10 pts) :::{exercise} :label: pr-nm-meoh-fit-sse Score the A2 fit **over the same window** with the sum of squared errors the chapter uses: $$ \mathrm{SSE} \;=\; \sum_i \left(y_i - \hat{y}_i\right)^2 \qquad\text{where}\qquad \hat{y} = \bar{\bar{X}}\vec{w} $$ Assign it to `sse_single`. For scale: predicting nothing but the mean absorbance gives SSE = 33.35, so that is the number to beat. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE sse_single = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] grader.check("q3") ``` --- ## Part B — Visualization (35 pts) :::{exercise} :label: pr-nm-meoh-basis-tune One Gaussian left SSE at 11.8. Stay in the **same 950–1120 cm⁻¹ window** and find out what it takes to do better, using the two basis families from the chapter. Everything here is an ordinary linear least-squares fit: build $\bar{\bar{X}}$, then solve $\bar{\bar{X}}^T\bar{\bar{X}}\vec{w} = \bar{\bar{X}}^T\vec{y}$ with `np.linalg.solve`. **1. Polynomials, and why the origin matters.** Use the chapter's `vandermonde(x, order)` for `order = 3, 5, 7, 9, 11, 13` columns, twice over: - once on the raw wavenumbers, and - once on **centered** wavenumbers, `x_band - x_band.mean()`. For each fit record the SSE **and** the condition number of $\bar{\bar{X}}^T\bar{\bar{X}}$ (`np.linalg.cond`, from {doc}`Linear Algebra `). The two versions describe the same family of curves — centering only moves the origin — so any difference between them is not about the model. **2. Gaussians, placed by hand.** Extend the chapter's two-column construction: allocate `np.zeros((len(x_band), N + 1))`, fill column $j$ with a Gaussian of your chosen center $c_j$ and width $\sigma_j$, and make the last column ones for the baseline. Only the amplitudes come from the solve — the centers and widths are yours to choose. Work from `N = 1` to `N = 6`. For each `N`, move and widen the Gaussians by eye until the fit tracks the data as well as you can get it, then record SSE and the condition number. Say in one line what you were trying to capture with each new component. **3. Plot.** The data with your `N = 1`, `N = 3` and `N = 6` fits overlaid, and beneath it, on shared x-axes, the differences $y - \hat{y}$ for those three. **4. Interpret**, in **4–6 sentences**: - Raw-wavenumber polynomials stop improving past about 7 columns, while centered ones keep improving. Explain this using your condition numbers. What is a condition number of $10^{40}$ telling you about the weights that came back? - The Gaussian condition numbers stay below about 40 no matter how many components you add, which is some thirty orders of magnitude better. What is different about those columns? - Your Gaussian SSE should drop sharply somewhere and then nearly stop. Where, and what does that say about the band? Label your axes with units. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE fig, ax = plt.subplots() ``` --- ## Part C — Open Ended (35 pts) So far every basis function has been a Gaussian. That was a choice, not a necessity — the general linear model works with any fixed set of columns, and a different column shape may describe the same data with fewer of them. A widely used alternative is the **Voigt profile**. Where a Gaussian falls off as $e^{-x^2}$, the Voigt profile decays much more slowly away from its center, so it has noticeably heavier tails for the same central width. It is controlled by two width parameters instead of one: $\sigma$ sets the Gaussian-like core and $\gamma$ sets how heavy the tails are, with $\gamma = 0$ recovering a Gaussian exactly. Use this helper: ```{code-cell} ipython3 from scipy.special import voigt_profile def voigt(x, center, sigma, gamma): """One Voigt profile. gamma = 0 gives a Gaussian; larger gamma gives heavier tails.""" return voigt_profile(x - center, sigma, gamma) ``` The important point for this course is that **nothing about the fitting changes.** Fix each center, $\sigma$ and $\gamma$ by hand as you did in Part B, put one profile in each column of $\bar{\bar{X}}$, and the amplitudes are still the solution of the same normal equations. The shape parameters are chosen, not fitted. ::::{exercise} :label: pr-nm-meoh-count-components Use a Voigt basis on the same 950–1120 cm⁻¹ window to answer the question this problem set has been building toward: **How many distinct components does this band contain?** Support your answer with: 1. Voigt fits for `N = 1, 2, 3` and at least one larger `N`, each with hand-chosen centers, $\sigma$ and $\gamma$, reporting SSE and the condition number. Compare against your Part B Gaussian results at the same `N`. 2. The fitted amplitudes of your best model. They are not all alike, and the pattern matters more than the individual values. 3. Evidence that your chosen `N` is not merely the largest you tried: show what one more component buys, and show the difference plot $y - \hat{y}$ has no obvious structure left. 4. A short written argument (**one paragraph**). "Number of components" can mean the number of visible maxima, or the number of columns the data actually justifies, or the number of underlying physical features. These need not be the same number, and saying so with evidence is a better answer than picking one. :::{note} Domain knowledge of spectroscopy may be helpful here, but it is **not required** — the question can be answered entirely from the fits, the SSE and the difference plots. If you do know some spectroscopy, say what it adds; if you do not, do not go looking for it. ::: There is more than one defensible answer. You are graded on the reasoning and the evidence. :::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ``` --- ## Summary - Part A built a two-column design matrix by hand and solved the normal equations, exactly as the chapter does, and scored the result with the sum of squared errors: 11.8 against a mean-only baseline of 33.3. - Part B set polynomial and Gaussian bases against each other on the same band. Polynomials on raw wavenumbers stall once the condition number passes what double precision can carry; centering the same model recovers thirty orders of magnitude of conditioning. A localized Gaussian basis stays well conditioned at any size, and its SSE stops improving once there is one component per resolved maximum. - Part C changed the shape of the basis functions rather than their number, and found that three Voigt profiles describe the band better than thirteen polynomial columns or six Gaussians — then asked what "three components" actually means. ## Additional Reading 1. NIST Chemistry WebBook, SRD 69 — [methanol IR spectrum](https://webbook.nist.gov/cgi/cbook.cgi?ID=C67561&Type=IR-SPEC), Coblentz Society collection no. 8791. 2. [`numpy.linalg.cond`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.cond.html) — the diagnostic behind Part B's raw-versus-centered comparison. 3. [`scipy.special.voigt_profile`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.voigt_profile.html), used here only to evaluate a fixed shape.