Problems: Generative Models

Get this problem set

Download everything (Topic5.4-Generative_Models_Problems.zip) — the notebook and xps_carbon_monitoring.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 Generative 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 carbon XPS spectra from the Clustering problem set. You do not need to have done that set first; what you need from it is below.

X-ray photoelectron spectroscopy measures the binding energy of core electrons, which shifts as the atom’s chemical environment changes. For carbon, the C 1s line moves to higher binding energy as the carbon becomes more oxidized:

component	binding energy	environment
C–C / C–H	284.8 eV	hydrocarbon
C–O	286.3 eV	alcohol, ether
C=O	287.8 eV	carbonyl
O–C=O	289.0 eV	carboxyl, ester

data/xps_carbon_monitoring.csv holds 520 simulated C 1s spectra from 280 to 295 eV in 0.05 eV steps, in two batches:

• batch == 'library' — 400 reference spectra, 100 from each of four surface chemistries, with the true class in label (0–3) and chemistry. Every spectrum is a mixture: a dominant component plus adventitious hydrocarbon and some of the neighboring oxidation states, so the four classes overlap.

• batch == 'new' — 120 spectra from a later campaign on the same instrument. Nobody has looked at them yet, so label and chemistry are blank. They are the subject of Part C.

Parts A and B use only the library.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import KernelDensity
try:
    plt.style.use('../settings/plot_style.mplstyle')   # available inside the book
except OSError:
    pass                                              # downloaded notebook: use defaults

df = pd.read_csv('data/xps_carbon_monitoring.csv')
ecols = [c for c in df.columns if c.startswith('E_')]
energy = np.array([float(c[2:]) for c in ecols])
lib = df[df['batch'] == 'library']
new = df[df['batch'] == 'new']
y_lib = lib['label'].astype(int).values
print(f"library: {len(lib)} spectra;  new batch: {len(new)} spectra;  {len(ecols)} channels")

library: 400 spectra;  new batch: 120 spectra;  301 channels

The preprocessing is the one from the Clustering set, and is given here. It subtracts a linear background drawn between the means of the first and last 20 channels, then scales each spectrum to sum to 100, so that each channel is a percentage of that sample’s total signal. The density models in Parts A and B are fit to the first three PCA scores, Z_lib, rather than to all 301 channels. The Clustering set’s A3 asks why three is enough.

def preprocess(raw):
    """Linear background subtraction, then scale every spectrum to sum to 100."""
    bg = np.linspace(raw[:, :20].mean(axis=1), raw[:, -20:].mean(axis=1), raw.shape[1]).T
    sub = raw - bg
    return 100.0 * sub / sub.sum(axis=1)[:, None]

X_lib = preprocess(lib[ecols].values)
X_new = preprocess(new[ecols].values)

pca = PCA(n_components=3, random_state=0).fit(X_lib)   # fit to the library only
Z_lib = pca.transform(X_lib)
Z_new = pca.transform(X_new)
print(f"Z_lib {Z_lib.shape}, Z_new {Z_new.shape}")

Z_lib (400, 3), Z_new (120, 3)

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': (6, ['41b7574f6abc', '228cbaaacbc7', '128452833d44', '1430ee8a8510', '243fe7404b25', '47e7b4509d40', '7c21e1e8fed4', '32ab71fb401f', '829ed3e6ac20', 'ad829fee4e62', '84d1bc2bee5c']),
    'q2': (4, ['6ed02f37bc7b', '813fa72c2d1b', 'c7a2b60c3538', '296a6bdba4e6', '2a2429d9178e', '04c5418d2192', 'cff78b96728d', '961500042fc7', '9024a2d603ed', '17f95c430532', '2b247686f82c', '6ee13507013b', 'fee77879547b', 'f54fefb41403', 'e824af8e78c8', '69087a95b4d3', '9b152c5291e1']),
    'q3': (1, ['d05b348c3e33', '3eb63dc3d942', 'c2c0132ab5b2', 'e6ef08438bd7', 'a0fd81aca156', 'dc03145f08fd', '12a2f0b59673']),
}

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)

The chapter’s two density models are a Gaussian mixture, with a number of components you choose, and a kernel density estimate, with a bandwidth you choose. Each has one setting that controls how detailed it is, and the chapter gives each a principled way to set it: the BIC for the mixture, cross-validated log-likelihood for the KDE. A1 and A2 use those on the library. A3 then turns the mixture into the thing Part C needs, a rule for deciding whether a new spectrum looks like the library.

A1. How many Gaussians does the library need? (10 pts)

Exercise 1

For n = 1, 2, …, 8, fit

GaussianMixture(n, covariance_type='full', n_init=10, random_state=0)

to Z_lib and compute its BIC on Z_lib.

Assign the number of components with the lowest BIC to n_bic.

# YOUR CODE HERE
n_bic = ...

check("q1", n_bic)

A2. Choose a KDE bandwidth by cross-validation (10 pts)

Exercise 2

Use GridSearchCV to choose the bandwidth of a Gaussian KernelDensity model of Z_lib, searching

np.logspace(-1, 1, 41)

with cv=5, as in the chapter’s Bandwidth Selection section.

Assign the selected bandwidth to h_cv.

# YOUR CODE HERE
h_cv = ...

check("q2", h_cv)

A3. An alarm threshold from the library (10 pts)

The chapter’s anomaly-detection section scores every sample with score_samples, which returns its log-probability \(\ln P(\mathbf{x})\) under the fitted density, and flags the lowest 1% as unusual. Applied to a reference set, that gives a threshold: a new sample whose log-probability falls below it is less typical than 99% of the library.

Exercise 3

Fit GaussianMixture(n_bic, covariance_type='full', n_init=10, random_state=0) to Z_lib, using your answer to A1 for the number of components. Compute the log-probability of every library spectrum with score_samples, and take its 1st percentile with np.percentile.

Assign that value to ll_limit.

# YOUR CODE HERE
ll_limit = ...

check("q3", ll_limit)

---

Part B — Visualization (35 pts)

The clustering set ended with an unanswered question. k-means recovered the four chemistries with a purity well short of 1 at \(k = 4\), and it was not clear whether to blame the clustering or the data. A generative classifier can settle it. Fit a density to each chemistry separately, using the labels this time, and combine them with Bayes’ theorem as in the chapter’s not-so-naive Bayes section. If even that classifier cannot do much better than the clustering did, the classes overlap in the data, and no clustering algorithm could have separated them.

The second question is about the other thing a generative model is for. A3’s mixture is a model of the library’s density, and you can draw new samples from it. Do the samples look like XPS spectra?

Exercise 4

1. What would labels have bought? Split the library with

train_test_split(Z_lib, y_lib, test_size=0.3, random_state=0)

and classify the test spectra with the chapter’s not-so-naive Bayes classifier: a Gaussian KDE fit to each chemistry’s training spectra, plus the log of that chemistry’s prior. Do this for 21 bandwidths, np.logspace(-1, 1, 21), and plot test accuracy against bandwidth on a logarithmic x-axis. On the same axes, draw horizontal lines for:

• the test accuracy of GaussianNB on the same split;

• the purity of KMeans(4, n_init=50, random_state=0) fit to X_lib, scored against y_lib as in the Clustering set’s A2 (the sum of the confusion matrix’s column maxima, divided by the number of spectra).

Mark h_cv from A2. Report the best accuracy and the bandwidth it occurs at.

2. Generate spectra. Draw 400 samples from A3’s mixture with .sample(400), and map them back to spectra with pca.inverse_transform. Make two panels:

• the PC1–PC2 scores of the real library and of the samples, overlaid;

• five real and five synthetic spectra, from 282 to 292 eV.

3. Interpret. Answer each in one or two sentences:

• Compare your best accuracy from task 1 with the purity. Was the clustering’s purity a failure of k-means, or a property of the data?

• In the chapter, Gaussian naive Bayes lost to the not-so-naive classifier by 16 points on the digits. Compare the gap you found here, and explain the difference.

• Name one way the synthetic spectra differ from real ones, and explain where the difference comes from.

Label all axes, with units where there are any.

# YOUR CODE HERE

---

Part C — Open Ended (35 pts)

Now the new batch. The instrument’s log shows three things went wrong at some point during the campaign, and nobody recorded which spectra were affected:

• Sample charging. An insulating sample charged up, shifting its whole spectrum to higher binding energy by an unknown amount.

• Fluoropolymer contamination. PTFE tape was used on the sample holder for part of the campaign. Fluorinated carbon (CF\(_2\)) appears near 292 eV, where no library spectrum has a peak.

• A new oxidation treatment. A few samples came from a harsher treatment than anything in the library.

Your job is to find the affected spectra using only the library as a reference. The mixture from A3 is a natural first tool, but it has a blind spot: it only ever sees the three PCA scores. Anything a spectrum does outside those three directions is removed by the projection before the density model sees it. The squared prediction error from the Dimensionality Reduction problem set measures exactly that part:

\[ \text{SPE} = \lVert \mathbf{x} - \hat{\mathbf{x}} \rVert^2, \qquad \hat{\mathbf{x}} = \texttt{pca\_k.inverse\_transform}(\texttt{pca\_k.transform}(\mathbf{x})) \]

where pca_k is a PCA with \(k\) components fit to X_lib. So the two alarms ask different questions: the density alarm asks whether a spectrum sits in a typical place inside the library’s subspace, and the residual alarm asks whether it sits in that subspace at all.

Exercise 5

1. Density alarm. Score Z_new with A3’s mixture and flag every new spectrum whose log-probability is below ll_limit. Plot histograms of the library’s and the new batch’s log-probabilities on one set of axes with the threshold marked, and report how many new spectra are flagged.

2. Residual alarm. For \(k = 3, 4, \ldots, 10\), fit PCA(k, random_state=0) to X_lib and set the SPE limit at the 99th percentile of the library’s SPE. Plot the limit against \(k\), choose a \(k\), and explain the choice in one sentence. Then report how many new spectra are above the limit at your \(k\).

3. Diagnose. Sort the flagged spectra into three groups: flagged by the density alarm only, by the residual alarm only, and by both. For each flagged spectrum, decide which of the three incidents explains it, or whether it is a false alarm, using evidence from the spectrum itself. Report the count for each diagnosis. Plotting a few of each group against the library’s mean spectrum is a good place to start.

4. Recommend. In three sentences, tell the lab how to screen future batches. Say which alarm or combination of alarms to use, which incident your screen still misses and why, and how many false alarms they should expect per 100 ordinary spectra. Use your numbers.

There is more than one defensible diagnosis method. You are graded on whether your evidence supports each diagnosis.

# YOUR CODE HERE

---

Additional Reading

1. NIST X-ray Photoelectron Spectroscopy Database, SRD 20 Version 5.0 — the source of the C 1s chemical shifts used to build this dataset.

2. scikit-learn documentation, Gaussian mixture models and density estimation.

3. J. E. Jackson and G. S. Mudholkar, “Control procedures for residuals associated with principal component analysis,” Technometrics 21, 341–349 (1979) — the origin of the residual (Q, or SPE) statistic used in Part C.