Problems: Linear Regression

Problems: Linear Regression#

Get this problem set

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.

Before you start

This problem set accompanies 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.

%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")

3567 points, 463.4-3806.6 cm^-1

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');
../_images/b7ed81a9f1302bca2efcdad2245d163459124f5fc2e401ead755eabd9744d903.png

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.

# Run this once. It defines `check`, used after each Part A answer below.
import hashlib

_EXPECTED = {
    'q1': (1, ['887bc508ac93', 'd833e1b76c6e', '7af23019f7b7', 'a573c201482e', '51c5722120e6']),
    'q2': (3, ['eda665e625f7', '17d8aff0af34', 'fd6b59ebdc0b', '3d1df7e03fd4', '885056be2e0f', 'f6c1a2f88e57', '365eee4e8b61', 'dc1dd5d51876', '7e6beeac81b3', '23c200bc0fbe', 'fd29cbda40f8', '5b64dda48afb', 'd56953b73b41', '618ae9d0216d', 'cb1ac17bc882', '9545ae7d6017', '91fec1918875', 'b07c03b41378', '00dd3025974c', '6bd3a06c3285', '0f90e9f39d66']),
    'q3': (3, ['fa1a484807cf', '57f95f0cdb0f', '36ffd9ff38c4', 'a065cea9a80c', '8db02c74254f']),
}

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. 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 1

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.

# YOUR CODE HERE
peak_wavenumber = ...

check("q1", peak_wavenumber)

A2. Fit one Gaussian to the band (10 pts)#

Exercise 2

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.

# YOUR CODE HERE
w_gauss = ...

check("q2", w_gauss)

A3. Score the fit (10 pts)#

Exercise 3

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.

# YOUR CODE HERE
sse_single = ...

check("q3", sse_single)

Part B — Visualization (35 pts)#

Exercise 4

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, twice over:

  • once on the raw wavenumbers, and

  • once on centered wavenumbers, x_band - x_band.mean().

order counts columns, not degree

vandermonde(x, order) returns order columns holding \(x^0, x^1, \ldots, x^{order-1}\), so order = 5 gives \(1, x, x^2, x^3, x^4\) — five columns spanning polynomials of degree up to four. Watch out for the clash with Linear Algebra, whose “2nd-order Vandermonde matrix” demonstration counts the other way and produces three columns. Here, and in the tables you are about to build, order is the column count.

For each fit record the SSE and the condition number of \(\bar{\bar{X}}^T\bar{\bar{X}}\) (np.linalg.cond, from 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. For each N, allocate a design matrix of N + 1 columns with np.zeros((len(x_band), N + 1)), fill columns \(0 \ldots N-1\) with Gaussians of your chosen centers \(c_j\) and widths \(\sigma_j\), and set the last column to 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 inspecting the spectrum until the fit tracks the data as well as you can get it, then record SSE and the condition number. Your choices will not match anyone else’s, and they are not supposed to — aim for a reasonable representation of the features you can see. Give one short description of what each new Gaussian is capturing.

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. Three numbers from your own tables, then two sentences. Every number here is read off the results you produced — there is no target value to hit.

  1. Centering. At the largest order you ran, divide the raw-wavenumber condition number by the centered one. (Both are the same model on the same data; only the origin moved.)

  2. Basis. The largest condition number across all your Gaussian fits, next to the largest across your polynomial fits. Compare orders of magnitude, not digits.

  3. Saturation. For each N, the ratio \(\mathrm{SSE}(N) / \mathrm{SSE}(N+1)\). Which single step is the biggest?

Then, in one or two sentences: double precision carries about 16 significant digits, so which of your polynomial fits can still have returned meaningful weights — and what is it about the Gaussian columns that keeps their condition numbers where they are?

Label your axes with units.

# 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:

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 5

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. Three counts, then two sentences. “Number of components” can mean three different things, so count all three:

    • how many resolved maxima you can pick out in the band by eye;

    • how many columns your fits actually justify — the N past which SSE stops improving materially, by the same ratio test you used in Part B;

    • how many of your fitted amplitudes are within a factor of about 3 of the largest one.

    Then, in one or two sentences: do your three numbers agree, and which one would you report if someone asked how many components this band has?

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.

# YOUR CODE HERE

Additional Reading#

  1. NIST Chemistry WebBook, SRD 69 — methanol IR spectrum, Coblentz Society collection no. 8791.

  2. numpy.linalg.cond — the diagnostic behind Part B’s raw-versus-centered comparison.

  3. scipy.special.voigt_profile, used here only to evaluate a fixed shape.

Contents