--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: Generative Models :::{admonition} Get this problem set :class: seealso {download}`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. ::: :::{admonition} Before you start :class: tip This problem set accompanies {doc}`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 {doc}`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. ```{code-cell} ipython3 %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") ``` 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. ```{code-cell} ipython3 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}") ``` 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': (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} :label: pr-eda-xps-gmmbic 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE n_bic = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q1", n_bic) ``` ### A2. Choose a KDE bandwidth by cross-validation (10 pts) :::{exercise} :label: pr-eda-xps-kdebw 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE h_cv = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-xps-lllimit 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`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ll_limit = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] 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} :label: pr-eda-xps-labels **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. ::: ```{code-cell} ipython3 :tags: [skip-execution] # 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} :label: pr-eda-xps-monitor **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. :::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ``` --- ## Additional Reading 1. NIST X-ray Photoelectron Spectroscopy Database, [SRD 20 Version 5.0](https://srdata.nist.gov/xps/) — the source of the C 1s chemical shifts used to build this dataset. 2. scikit-learn documentation, [Gaussian mixture models](https://scikit-learn.org/stable/modules/mixture.html) and [density estimation](https://scikit-learn.org/stable/modules/density.html). 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.