{ "cells": [ { "cell_type": "markdown", "id": "h1-00", "metadata": {}, "source": [ "# Hands-On 1: Calibrating a sensor array\n", "**Signal Processing · ECE-Y523 · University of Patras**
\n", "Konstantinos Chatzilygeroudis · 29 October 2026\n", "\n", "Ten probes measure the same temperature. Each reports biased raw ADC counts,\n", "and their noise is correlated. Build a calibration, report one temperature with\n", "an uncertainty, and investigate a recurring disagreement between probes.\n", "\n", "![Ten probes in a reference bath with a shared ADC, fan and reference thermometer; temperature programmes for the three runs](figs/h1_setup_en.png)\n", "\n", "The rack and reference bath are shown on the left. The three temperature\n", "programmes on the right connect the physical setup to the data files.\n", "\n", "| Stage | Question |\n", "|---|---|\n", "| 1. Calibrate | How do counts become °C? |\n", "| 2. Rescue a poor calibration | What can ridge and a physical prior repair? |\n", "| 3. Understand the noise | Which errors are shared? |\n", "| 4. Report a temperature | What uncertainty survives averaging? |\n", "| 5. Diagnose faults | Is the disagreement persistent, and which probe explains it? |\n", "| 6. Validate | Does the final report agree with the reference? |\n", "\n", "**Files.** `calibration.csv`: 400 snapshots at 10, 20, 30, 40 °C, ten probes.\n", "`recal_quick.csv` and `recal_full.csv`: two checks of probe 3.\n", "`operations.csv`: a new 300-snapshot run. In the array files the columns are\n", "`snapshot, T_ref, probe_1, …, probe_10`; the recalibration files have `T_ref, probe_3`.\n", "A snapshot is one simultaneous reading of the rack. During calibration, a stirred\n", "reference bath keeps all probes at the accurately known temperature `T_ref`.\n", "\n", "**The walkthrough.** Raw counts become calibrated temperatures, a combined\n", "estimate with uncertainty, and a diagnosis of recurring disagreement.\n", "The operational reference is revealed only in stage 6.\n", "\n", "**Running the notebook.** Run from `Hands-On 1` with NumPy, SciPy and Matplotlib.\n", "The handout follows these six stages. Optional exercises are in\n", "the optional notebook; proofs are in\n", "the optional appendices.\n", "\n", "**What the measurements mean.** The probes share an ADC reference, a power supply\n", "and airflow along the rack. They can therefore make similar errors at the same\n", "time. A shared voltage disturbance moves most readings together; an airflow\n", "change can affect the two ends differently. Each probe also has its own noise.\n", "Calibration addresses offsets and gains. It does not remove random errors from\n", "future readings, and ten probes need not provide ten independent measurements.\n", "\n", "The reference bath makes the calibration input accurately known. During normal\n", "operation that reference is unavailable, so the final temperature and the fault\n", "diagnosis must come from the probes themselves. Keeping the validation reference\n", "hidden initially separates what can be inferred from the rack from what can be\n", "checked only when the true temperature is known.\n", "\n", "**Notation.** Matrices use bold uppercase symbols, such as $\\boldsymbol{X}$ and\n", "$\\boldsymbol{S}$; vectors use bold lowercase symbols, such as $\\boldsymbol{w}$\n", "and $\\boldsymbol{r}$. Scalars such as $T_k$, $a_j$ and $b_j$ remain plain.\n", "Array shapes and units connect the formulas to the NumPy operations below." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-01", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import session_plotting as plots\n", "\n", "plots.configure(\"en\")\n", "np.set_printoptions(precision=3, suppress=True)" ] }, { "cell_type": "markdown", "id": "h1-02", "metadata": {}, "source": [ "## 1. Calibrate: counts → °C\n", "*Lectures 1 and 5*\n", "\n", "The reference temperature is accurately known; the counts contain the noise.\n", "Fit $y_{kj}=a_j+b_j T_k+n_{kj}$, with $\\boldsymbol{X}_{k,:}=[1,T_k]$.\n", "The same $N\\times2$ design matrix serves all probes:\n", "$\\boldsymbol{Y}=\\boldsymbol{X}\\boldsymbol{\\Theta}+\\text{noise}$, where $\\boldsymbol{\\Theta}$ is $2\\times C$.\n", "Use `lstsq(X, Y)`, then invert each fitted line.\n", "\n", "**From one probe to the matrix equation.** For probe $j$, collect all its counts\n", "in $\\boldsymbol{y}_j\\in\\mathbb{R}^N$ and its two coefficients in\n", "$\\boldsymbol{\\theta}_j=[a_j,b_j]^{\\top}$. Each row of the design matrix contains\n", "a constant for the offset and the reference temperature for the gain:\n", "\n", "$$\n", "\\boldsymbol{X}=\\begin{bmatrix}1&T_1\\\\ \\vdots&\\vdots\\\\1&T_N\\end{bmatrix},\n", "\\qquad \\boldsymbol{y}_j=\\boldsymbol{X}\\boldsymbol{\\theta}_j+\\boldsymbol{n}_j.\n", "$$\n", "\n", "All probes see the same reference temperatures, so their data are columns of\n", "$\\boldsymbol{Y}$ and their coefficients are columns of $\\boldsymbol{\\Theta}$.\n", "`lstsq` accepts this matrix right-hand side. It solves the ten fits together;\n", "it does not force their offsets or gains to be equal. The familiar expression\n", "$(\\boldsymbol{X}^{\\top}\\boldsymbol{X})^{-1}\\boldsymbol{X}^{\\top}\\boldsymbol{Y}$\n", "describes the solution when the columns are independent. `lstsq` computes it\n", "without explicitly forming that inverse.\n", "\n", "**Why this direction of fit?** Least squares treats the design matrix as known\n", "and puts the measurement noise in the output. Here the accurately measured bath\n", "temperature belongs in $\\boldsymbol{X}$ and noisy counts belong in\n", "$\\boldsymbol{Y}$. Fitting temperature directly against noisy counts changes that\n", "assumption. Instead, first fit the physical response, then invert it:\n", "$\\hat T_{kj}=(y_{kj}-\\hat a_j)/\\hat b_j$." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-03", "metadata": {}, "outputs": [], "source": [ "cal = np.loadtxt(\"data/calibration.csv\", delimiter=\",\", skiprows=1)\n", "T_cal, Y = cal[:, 1], cal[:, 2:]\n", "N, C = Y.shape\n", "X = np.column_stack([np.ones(N), T_cal])\n", "Theta = np.linalg.lstsq(X, Y, rcond=None)[0]\n", "a, b = Theta\n", "\n", "def to_temperature(counts, a, b):\n", " return (counts - a) / b\n", "\n", "R = Y - X @ Theta # counts; N x C\n", "E = R / b # °C; N x C\n", "p = X.shape[1]\n", "sigma_counts = np.sqrt((R**2).sum(axis=0) / (N - p))\n", "print(\"nominal snapshot 100 [°C]:\", (Y[100] - 200) / 40)\n", "print(\"calibrated snapshot 100 [°C]:\", to_temperature(Y[100], a, b))\n", "print(\"[a; b] =\", Theta, \"\\nsd [°C] =\", sigma_counts / b)\n", "plots.calibration(T_cal, Y, Theta, E)" ] }, { "cell_type": "markdown", "id": "h1-04", "metadata": {}, "source": [ "**Reading the result.** The first row of `Theta` contains offsets in counts;\n", "the second contains gains in counts/°C. NumPy unpacks these rows into `a` and\n", "`b`. Subtracting the fitted counts gives $\\boldsymbol{R}$, an $N\\times C$\n", "residual matrix. Dividing its columns by their gains produces\n", "$\\boldsymbol{E}$ in °C; broadcasting applies each gain to its own probe.\n", "\n", "The left panel checks whether straight lines describe the raw response. The\n", "right panel asks whether anything systematic remains at the four bath\n", "temperatures. A curved or bath-dependent pattern would call the model into\n", "question. The residual sd divides by $N-p$, with $p=2$, because the fit estimated\n", "two coefficients. This is a noise estimate, rather than an uncertainty on the\n", "fitted gain. The nominal-versus-calibrated snapshot illustrates bias removal;\n", "the remaining spread illustrates noise that calibration cannot eliminate.\n", "\n", "**Discussion.** How do the nominal $(a,b)=(200,40)$ and fitted calibrations\n", "differ for one snapshot? Which changes: bias, scatter or both?\n", "\n", "**Checkpoint.** Counts have become °C. Each probe's remaining noise is about\n", "0.2 °C. Small residual trends support the model; they do not establish Gaussianity.\n", "The repeated correlated errors will matter in stage 3." ] }, { "cell_type": "markdown", "id": "h1-05", "metadata": {}, "source": [ "## 2. Rescue a poor calibration\n", "*Lectures 2, 5 and 6*\n", "\n", "Probe 3 has a short recalibration at only 25.0 and 25.2 °C. Those temperatures give\n", "almost no lever arm for estimating a gain. Compare three fits:\n", "\n", "- **LS:** minimise $\\|\\boldsymbol{y}-\\boldsymbol{X}_q\\boldsymbol{\\theta}\\|^2$.\n", "- **Plain ridge:** add $\\lambda\\|\\boldsymbol{\\theta}\\|^2$ with $\\lambda=1$.\n", "- **MAP:** use a Gaussian prior from the other nine probes, with a separate\n", " scale for offset and gain. With noise sd $\\sigma$ and prior precision $\\boldsymbol{\\Lambda}_0$,\n", " solve $(\\boldsymbol{X}_q^{\\top} \\boldsymbol{X}_q/\\sigma^2+\\boldsymbol{\\Lambda}_0)\\boldsymbol{\\theta}=\\boldsymbol{X}_q^{\\top} \\boldsymbol{y}/\\sigma^2+\\boldsymbol{\\Lambda}_0\\boldsymbol{\\mu}_0$.\n", "\n", "This is a controlled comparison. The downstream rack calibration remains the\n", "full four-bath fit from stage 1.\n", "\n", "**The missing lever arm.** The two columns of $\\boldsymbol{X}_q$ are almost\n", "parallel: near 25 °C, the temperature column is almost 25 times the column of\n", "ones. Increasing $b$ and decreasing $a$ can leave $a+25b$ nearly unchanged.\n", "Consequently, many very different lines fit the narrow measurement range.\n", "For independent count noise of variance $\\sigma^2$, the gain variance is\n", "\n", "$$\n", "\\operatorname{Var}(\\hat b)=\\frac{\\sigma^2}{s_{TT}},\\qquad\n", "s_{TT}=\\sum_k(T_k-\\bar T)^2.\n", "$$\n", "\n", "The full rack calibration has $s_{TT}=50\\,000$ °C², whereas the short recalibration has\n", "only $0.2$ °C². Similar noise therefore gives a gain sd about 500 times larger.\n", "A second temperature helps only if it is sufficiently separated from the first.\n", "\n", "**What the regularisation claims.** Plain ridge penalises the numerical sizes\n", "of both coefficients equally and shrinks towards zero. Yet an offset in counts\n", "and a gain in counts/°C have different scales and physical meanings. The MAP\n", "comparison uses the other probes to estimate a plausible coefficient mean and\n", "separate spreads. `Lambda0` is the inverse prior covariance. Larger prior\n", "variances mean weaker precision, so the data have more influence. This is an\n", "approximate empirical prior from nine similar probes, not certain knowledge\n", "about the serviced probe. The count-noise scale comes from the full calibration." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-06", "metadata": {}, "outputs": [], "source": [ "quick = np.loadtxt(\"data/recal_quick.csv\", delimiter=\",\", skiprows=1)\n", "T_q, y_q = quick.T\n", "Xq = np.column_stack([np.ones(len(T_q)), T_q])\n", "theta_ls = np.linalg.lstsq(Xq, y_q, rcond=None)[0]\n", "theta_ridge = np.linalg.solve(Xq.T @ Xq + np.eye(2), Xq.T @ y_q)\n", "others = np.arange(C) != 2\n", "mu0 = Theta[:, others].mean(axis=1)\n", "prior_scale = 1.0 # comparison: 1.0 versus 5.0\n", "sd0 = prior_scale * Theta[:, others].std(axis=1, ddof=1)\n", "Lambda0 = np.diag(1 / sd0**2)\n", "sigma = sigma_counts[2]\n", "P = Xq.T @ Xq / sigma**2 + Lambda0\n", "theta_map = np.linalg.solve(P, Xq.T @ y_q / sigma**2 + Lambda0 @ mu0)\n", "\n", "full = np.loadtxt(\"data/recal_full.csv\", delimiter=\",\", skiprows=1)\n", "Xf = np.column_stack([np.ones(len(full)), full[:, 0]])\n", "theta_full = np.linalg.lstsq(Xf, full[:, 1], rcond=None)[0]\n", "fits = {\"LS\": theta_ls, \"ridge\": theta_ridge, \"MAP\": theta_map, \"full\": theta_full}\n", "print(\"cond(X), cond(Xq):\", np.linalg.cond(X), np.linalg.cond(Xq))\n", "for name, theta in fits.items():\n", " error = (np.array([1, 40]) @ (theta - theta_full)) / theta_full[1]\n", " print(f\"{name:6s}: [a,b] = {theta}, error(T=40) = {error:+.2f} °C\")\n", "plots.quick_check(T_q, y_q, fits, full[:, 0], full[:, 1])" ] }, { "cell_type": "markdown", "id": "h1-07", "metadata": {}, "source": [ "**Reading the comparison.** Zooming near 25 °C shows why the short-recalibration\n", "residuals alone are misleading. All three methods can explain the observed\n", "counts, while the full-range panel exposes their very different slopes.\n", "At 40 °C the LS line differs from the full recalibration by about +11 °C in equivalent\n", "temperature units; plain ridge reduces that discrepancy to about +2.4 °C and\n", "physical MAP to about −0.27 °C. These are differences between fitted response\n", "lines, obtained by dividing a count difference by the full-recalibration gain.\n", "\n", "The full recalibration is a comparison measurement, not exact parameter truth.\n", "Changing `prior_scale` from 1 to 5 makes each prior variance 25 times larger.\n", "The MAP estimate then relies more on the weak short-recalibration data and moves towards\n", "LS. This comparison explains regularisation; it does not replace the full-rack\n", "coefficients used in the remaining stages.\n", "\n", "**Discussion.** What changes when both prior standard deviations are multiplied\n", "by five? How would wider-spaced reference temperatures help?\n", "\n", "**Checkpoint.** A good local fit can extrapolate very badly. Ridge makes the\n", "solution stable; a physically meaningful prior determines where it stabilises.\n", "The appendix and the optional notebook develop the SVD and bias–variance explanations." ] }, { "cell_type": "markdown", "id": "h1-08", "metadata": {}, "source": [ "## 3. Understand the noise\n", "*Lectures 4 and 5*\n", "\n", "The fit used $p=2$ parameters per probe. Estimate the covariance in °C² with\n", "$\\boldsymbol{S}=\\boldsymbol{E}^{\\top} \\boldsymbol{E}/(N-p)$. Its eigenvectors describe noise patterns; its eigenvalues give\n", "their variances. Factor $\\boldsymbol{S}=\\boldsymbol{L}\\boldsymbol{L}^{\\top}$ and solve $\\boldsymbol{L}\\boldsymbol{q}=\\boldsymbol{e}$ to whiten an error vector.\n", "\n", "**What a covariance entry means.** The diagonal $S_{jj}$ estimates the noise\n", "variance of probe $j$. An off-diagonal $S_{ij}$ records how errors of probes\n", "$i$ and $j$ vary together in the *same snapshot*. The correlation divides out\n", "the two marginal scales, making this dependence easier to compare across\n", "pairs. This is dependence across probes; independence across different\n", "snapshots remains a separate modelling assumption.\n", "\n", "Writing $\\boldsymbol{S}=\\boldsymbol{U}\\boldsymbol{\\Lambda}\\boldsymbol{U}^{\\top}$\n", "gives orthonormal directions $\\boldsymbol{u}_i$. The scalar coordinate\n", "$p_i=\\boldsymbol{u}_i^{\\top}\\boldsymbol{e}$ has variance $\\lambda_i$ under the\n", "estimated model. A large eigenvalue means errors vary strongly along that\n", "pattern. The eigenvectors describe patterns rather than individual sensors.\n", "Their signs can reverse without changing the covariance or interpretation.\n", "\n", "**What the solve does.** For one column vector, whitening solves\n", "$\\boldsymbol{L}\\boldsymbol{q}=\\boldsymbol{e}$. The stored error vectors are rows,\n", "so `E.T` puts them into columns, and the final transpose restores the\n", "$N\\times C$ layout. The result is a matrix $\\boldsymbol{Q}$ whose rows contain\n", "whitened coordinates. These coordinates mix probes; a large coordinate need\n", "not name a particular faulty sensor." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-09", "metadata": {}, "outputs": [], "source": [ "S = E.T @ E / (N - p)\n", "eigenvalues, U = np.linalg.eigh(S)\n", "eigenvalues, U = eigenvalues[::-1], U[:, ::-1]\n", "L = np.linalg.cholesky(S)\n", "Q = np.linalg.solve(L, E.T).T\n", "print(\"variance fractions:\", eigenvalues / eigenvalues.sum())\n", "print(\"max |Cov(Q) - I|:\", np.max(np.abs(Q.T @ Q / (N - p) - np.eye(C))))\n", "plots.noise(S, eigenvalues, U)\n", "plots.whitening(E, Q)" ] }, { "cell_type": "markdown", "id": "h1-10", "metadata": {}, "source": [ "**Reading the figures.** Strong positive correlation means a high reading on\n", "one probe often accompanies a high reading on another. The first eigenpattern\n", "is almost uniform; the second changes along the rack. Their large variance\n", "fractions explain why the cloud of the first two probe errors is thin and\n", "tilted rather than round. Whitening removes this scale and tilt under the\n", "fitted covariance model.\n", "\n", "The identity check is deliberately done on the calibration residuals to show\n", "its limitation. The same $\\boldsymbol{E}$ constructed $\\boldsymbol{S}$, so\n", "identity follows algebraically even if the data are non-Gaussian. A model check\n", "needs new data. Later, the independent healthy validation errors will test\n", "whether the learned scale and cross-probe structure transfer to operation.\n", "\n", "**Discussion.** Which eigenvector is nearly uniform, and which varies from\n", "front to back? What physical disturbance could produce each? Eigenvector signs are arbitrary.\n", "\n", "**Checkpoint.** Two patterns carry about 94% of the variance. Whitening these\n", "same calibration data gives the identity by construction. An independent check\n", "comes in stage 6; this identity alone does not validate the noise model." ] }, { "cell_type": "markdown", "id": "h1-11", "metadata": {}, "source": [ "## 4. Report one temperature with uncertainty\n", "*Lecture 5*\n", "\n", "For readings $\\boldsymbol{t}=T\\boldsymbol{1}+\\boldsymbol{e}$, weighted least squares gives\n", "$\\boldsymbol{w}=\\boldsymbol{S}^{-1}\\boldsymbol{1}/(\\boldsymbol{1}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{1})$ and\n", "$\\hat T=\\boldsymbol{w}^{\\top}\\boldsymbol{t}$, with $\\mathrm{sd}(\\hat T)=(\\boldsymbol{1}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{1})^{-1/2}$.\n", "The function also returns the squared distance of the residual $\\boldsymbol{r}=\\boldsymbol{t}-\\hat T\\boldsymbol{1}$;\n", "we will use it in stage 5. All probes are included at this point.\n", "\n", "**The estimation problem.** Each calibrated probe is another observation of\n", "the same unknown scalar temperature. The design matrix for this new problem\n", "is the column $\\boldsymbol{1}$, not the two-column calibration matrix. WLS\n", "chooses the common level that minimises\n", "\n", "$$\n", "J(T)=(\\boldsymbol{t}-T\\boldsymbol{1})^{\\top}\\boldsymbol{S}^{-1}\n", " (\\boldsymbol{t}-T\\boldsymbol{1}).\n", "$$\n", "\n", "Its derivative gives\n", "$\\hat T=\\boldsymbol{1}^{\\top}\\boldsymbol{S}^{-1}\\boldsymbol{t}/\n", "(\\boldsymbol{1}^{\\top}\\boldsymbol{S}^{-1}\\boldsymbol{1})$.\n", "The weights in the formula above sum to one, so a noiseless common temperature\n", "is reproduced exactly. `Si1` computes $\\boldsymbol{S}^{-1}\\boldsymbol{1}$ through\n", "a solve; dividing by its sum normalises the weights.\n", "\n", "**The uncertainty calculation.** For any fixed weights,\n", "$\\operatorname{Var}(\\boldsymbol{w}^{\\top}\\boldsymbol{e})=\n", "\\boldsymbol{w}^{\\top}\\boldsymbol{S}\\boldsymbol{w}$.\n", "A plain mean has variance $\\boldsymbol{1}^{\\top}\\boldsymbol{S}\\boldsymbol{1}/C^2$;\n", "treating probes as independent drops every off-diagonal term. A shared error\n", "appears in every reading and survives any unbiased combination. More probes\n", "can reduce their independent noise, but cannot average away this common part." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-12", "metadata": {}, "outputs": [], "source": [ "def combine(T_probes, covariance):\n", " one = np.ones(covariance.shape[0])\n", " Si1 = np.linalg.solve(covariance, one)\n", " variance = 1 / (one @ Si1)\n", " T = T_probes @ Si1 * variance\n", " residual = T_probes - T[:, None]\n", " d2 = (residual * np.linalg.solve(covariance, residual.T).T).sum(axis=1)\n", " return T, np.sqrt(variance), d2\n", "\n", "ops = np.loadtxt(\"data/operations.csv\", delimiter=\",\", skiprows=1)\n", "k, counts = ops[:, 0], ops[:, 2:] # reference column reserved for stage 6\n", "T_probes = to_temperature(counts, a, b)\n", "T_all, sd_all, d2 = combine(T_probes, S)\n", "Si1 = np.linalg.solve(S, np.ones(C))\n", "w = Si1 / Si1.sum()\n", "sd_mean = np.sqrt(S.sum()) / C\n", "sd_naive = np.sqrt(np.diag(S).mean() / C)\n", "print(\"sum(w):\", w.sum(), \"sd WLS / mean / independent:\", sd_all, sd_mean, sd_naive)\n", "plots.weights(w, sd_all, sd_mean, sd_naive)" ] }, { "cell_type": "markdown", "id": "h1-13", "metadata": {}, "source": [ "**Reading the weights and bars.** Negative weights are allowed: WLS is a\n", "linear combination, not a vote constrained to positive fractions. A probe can\n", "help cancel a disturbance seen by other probes. However, the exact weights\n", "also depend on the estimated covariance, so their detailed values should not\n", "be overinterpreted as a physical ranking of sensor quality.\n", "\n", "The bars compare *standard deviations*, not 95% interval widths. Under the\n", "Gaussian model a nominal 95% interval uses approximately $\\hat T\\pm1.96s$.\n", "Here WLS and a correctly assessed plain mean have very similar uncertainties;\n", "the major mistake would be to attach an independent-probe error bar to the mean.\n", "The residual returned by `combine` removes the estimated common temperature\n", "and will reveal disagreements between probes in the next stage.\n", "\n", "**Discussion.** How does the uncertainty change when the off-diagonal entries\n", "of $\\boldsymbol{S}$ are set to zero? Which assumption has changed?\n", "\n", "**Checkpoint.** Shared noise survives averaging. Both WLS and the plain mean\n", "have uncertainty near 0.19 °C here; pretending independence gives about 0.07 °C.\n", "Negative WLS weights can cancel correlated disturbances. These are model-based\n", "uncertainties; validation and diagnosis are still needed." ] }, { "cell_type": "markdown", "id": "h1-14", "metadata": {}, "source": [ "## 5. Diagnose recurring disagreement\n", "*Lecture 4, plus a short introduction to χ²*\n", "\n", "For independent standard normal $z_i$, $\\sum_{i=1}^{\\nu}z_i^2$ has distribution\n", "$\\chi^2_\\nu$. A whitened residual has **$C-1$ degrees of freedom**: fitting one\n", "unknown temperature removes one direction. Thus compare\n", "$d^2=\\boldsymbol{r}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{r}$ with `stats.chi2.ppf(0.99, C-1)`. With estimated $\\boldsymbol{S}$, this is\n", "an approximate 1% per-snapshot threshold.\n", "\n", "**An alarm is not a faulty-probe verdict.** A healthy 300-snapshot run should\n", "produce about three alarms. First ask whether the *number* of alarms is unusually\n", "large. `binom.sf(m-1, n, alpha)` is the probability of at least `m` alarms in `n`\n", "independent healthy snapshots at rate `alpha`.\n", "\n", "For a run with excess alarms, omit each probe and refit the temperature (now\n", "$C-2$ degrees of freedom). For this **single recurring fault** demonstration, exclude\n", "one probe only if it uniquely clears most alarms and the remaining rack no\n", "longer has excess alarms. Otherwise retain the probes and report uncertainty\n", "about the diagnosis. This is an offline diagnostic rule, not a guaranteed fault\n", "classifier; the χ² and count checks rely on the Gaussian, constant-covariance,\n", "independent-snapshot model.\n", "\n", "**Why the reference is unnecessary.** Subtracting $\\hat T\\boldsymbol{1}$ removes\n", "the common level while retaining disagreements between probes. In whitened\n", "coordinates, fitting that level removes one direction. The remaining squared\n", "length has $C-1$ degrees of freedom. For $C=10$, the nominal 99% threshold is\n", "21.67. A large distance means the readings disagree along a direction with\n", "little expected noise. A common movement consistent with a real temperature\n", "change is largely absorbed into $\\hat T$.\n", "\n", "**From an alarm to a diagnosis.** `alarm_p` is a tail probability for the\n", "*number of alarms under the healthy model*. It is not a probability that a\n", "probe is faulty. Expected isolated alarms should not trigger automatic removal.\n", "If the count is unusually high, leave-one-out tests ask whether a particular\n", "probe explains the disagreement. Each test uses the covariance submatrix of\n", "the retained probes and refits the common temperature. Taking a submatrix of\n", "the full precision matrix would be a different operation.\n", "\n", "`blame[j]` records the fraction of original alarms cleared without probe $j$.\n", "A unique strong candidate is excluded only if the remaining run also passes\n", "the alarm-count check. This second condition prevents a partial explanation\n", "from being mistaken for a complete repair. Several faulty probes, a changed\n", "noise model or ambiguous scores can leave the cause unresolved." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-15", "metadata": {}, "outputs": [], "source": [ "def diagnose(T_probes, covariance, alpha=0.01, batch_alpha=0.01):\n", " n, channels = T_probes.shape\n", " _, _, distances = combine(T_probes, covariance)\n", " threshold = stats.chi2.ppf(1 - alpha, channels - 1)\n", " flags = distances > threshold\n", " alarm_p = stats.binom.sf(flags.sum() - 1, n, alpha)\n", " blame = np.zeros(channels)\n", " suspect = None\n", " remaining_p = None\n", "\n", " if alarm_p < batch_alpha:\n", " reduced_flags = []\n", " for j in range(channels):\n", " keep = np.arange(channels) != j\n", " _, _, reduced_d2 = combine(T_probes[:, keep], covariance[np.ix_(keep, keep)])\n", " reduced = reduced_d2 > stats.chi2.ppf(1 - alpha, channels - 2)\n", " reduced_flags.append(reduced)\n", " blame[j] = 1 - reduced[flags].mean()\n", " candidates = np.flatnonzero(blame > 0.5)\n", " if len(candidates) == 1:\n", " j = int(candidates[0])\n", " remaining_p = stats.binom.sf(reduced_flags[j].sum() - 1, n, alpha)\n", " if remaining_p >= batch_alpha:\n", " suspect = j\n", "\n", " return dict(d2=distances, flag=flags, threshold=threshold, blame=blame,\n", " alarm_p=alarm_p, remaining_p=remaining_p, suspect=suspect)\n", "\n", "diagnosis = diagnose(T_probes, S)\n", "suspect = diagnosis[\"suspect\"]\n", "print(\"alarms / snapshots:\", diagnosis[\"flag\"].sum(), \"/\", len(k))\n", "print(\"P(healthy alarm count >= observed):\", diagnosis[\"alarm_p\"])\n", "print(\"suspect:\", None if suspect is None else suspect + 1)\n", "print(\"first 100 snapshots, suspect:\", diagnose(T_probes[:100], S)[\"suspect\"])\n", "plots.detection(k, diagnosis[\"d2\"], diagnosis[\"threshold\"], diagnosis[\"flag\"], diagnosis[\"blame\"])" ] }, { "cell_type": "markdown", "id": "h1-16", "metadata": {}, "source": [ "**Reading the timeline and omission scores.** The left panel shows the residual\n", "distance for every snapshot; red points exceed the threshold. Its late cluster\n", "of alarms is more informative than an isolated early crossing. The right panel\n", "asks which omission clears that disagreement. A high bar is meaningful together\n", "with excess-alarm evidence and the remaining-rack check, not on its own.\n", "\n", "The first 100 snapshots illustrate this distinction: a healthy stretch can contain\n", "threshold crossings and still provide no recurring-fault evidence. The output\n", "keeps snapshot flags even when no probe is excluded. An alarm timeline can\n", "suggest when behaviour changed, but this batch calculation neither estimates\n", "an onset time nor implements an online shutdown rule.\n", "\n", "**Discussion.** Where does disagreement become persistent in the timeline?\n", "Does `diagnose(T_probes[:100], S)` justify excluding a probe? Why are a few alarms insufficient?\n", "\n", "**Checkpoint.** The detector has used no reference temperatures and no supplied\n", "fault time or probe number. Its timeline supports investigation of a change;\n", "it does not implement a change-point estimator. The diagnosis assumes one\n", "recurring probe fault and can leave the cause unresolved." ] }, { "cell_type": "markdown", "id": "h1-17", "metadata": {}, "source": [ "## 6. Validate the result\n", "\n", "Recombine the retained probes. **Now** reveal the reference temperature and\n", "check the error bars against it. Also check the calibration whitener on new\n", "errors. For this final validation only, the synthetic answer key says that\n", "probe 7 acquires a +0.30 °C offset from snapshot 200. Use the first 200 snapshots\n", "to check healthy behaviour; these labels played no part in diagnosis.\n", "\n", "**Three different checks.** First compare the final temperature with the\n", "reference: bias and scatter describe the errors of the reported estimate.\n", "Second check interval coverage. A nominal level $q$ uses\n", "$z_q=\\Phi^{-1}((1+q)/2)$ and counts snapshots satisfying\n", "$|\\hat T_k-T_{\\mathrm{ref},k}|\\le z_qs$. If those intervals are well calibrated,\n", "observed coverage follows the diagonal in the right panel. Third, whiten new\n", "healthy probe errors using the *old* calibration factor $\\boldsymbol{L}$.\n", "Their covariance should be near identity if the noise structure transfers.\n", "Re-estimating the whitener on those same new errors would defeat that check.\n", "\n", "These checks answer different questions. Reasonable mean errors do not ensure\n", "correct uncertainty, and correct average coverage does not identify a faulty\n", "probe. The known onset and probe index are now used to separate healthy and\n", "faulty validation snapshots. They were unavailable to the diagnostic function.\n", "The final combination uses the retained covariance submatrix and recomputes its\n", "weights, rather than deleting one old weight without renormalising." ] }, { "cell_type": "code", "execution_count": null, "id": "h1-18", "metadata": {}, "outputs": [], "source": [ "keep = np.ones(C, dtype=bool)\n", "if suspect is not None:\n", " keep[suspect] = False\n", "T_final, sd_final, _ = combine(T_probes[:, keep], S[np.ix_(keep, keep)])\n", "T_ref = ops[:, 1] # revealed after diagnosis\n", "T_mean = T_probes[:, keep].mean(axis=1)\n", "sd_naive_final = np.sqrt(np.diag(S)[keep].mean() / keep.sum())\n", "error = T_final - T_ref\n", "healthy = k < 200 # answer key: validation only\n", "Q_new = np.linalg.solve(L, (T_probes[healthy] - T_ref[healthy, None]).T).T\n", "print(\"held-out whitened variances:\", np.diag(np.cov(Q_new, rowvar=False)))\n", "print(\"alarm rate before / after:\", diagnosis[\"flag\"][healthy].mean(), diagnosis[\"flag\"][~healthy].mean())\n", "print(\"final sd [°C]:\", sd_final, \"95% coverage:\", np.mean(np.abs(error) <= 1.96 * sd_final))\n", "print(\"healthy-only suspect:\", diagnose(T_probes[healthy], S)[\"suspect\"])\n", "plots.report(k, T_ref, T_final, sd_final, T_mean, sd_naive_final)" ] }, { "cell_type": "markdown", "id": "h1-19", "metadata": {}, "source": [ "**Closing discussion.**\n", "1. What supports the calibration model?\n", "2. Why was the short-recalibration gain unreliable, and what did the prior contribute?\n", "3. What did the alarm count, leave-one-out scores and timeline show?\n", "4. What temperature uncertainty do you report, and what coverage did it achieve?\n", "5. Which assumptions could fail in a new run?\n", "\n", "**Expected validation.** The complete run supports probe 7 as the suspect.\n", "The reference-free test flags about 1.5% of healthy snapshots and 70% after the\n", "fault. With probe 7 excluded, sd is about 0.190 °C and 95% coverage is about 92%.\n", "The remaining undercoverage and held-out whitening departures remind us that\n", "$\\boldsymbol{S}$ is estimated. Calibration-coefficient uncertainty is small for this full\n", "calibration (appendix F), and has been omitted here. Exclusion is selected from\n", "this run; this validation is a worked check, not an independent evaluation of\n", "fault-classification performance.\n", "\n", "**Optional follow-up.** The optional notebook covers Kronecker stacking, coefficient\n", "uncertainty, repeated short recalibrations, SVD/ridge paths, posterior contours, χ²\n", "geometry and detector power. Further prompts: use fewer calibration snapshots,\n", "change the prior, or simulate several faults. Keep the reference hidden until\n", "you have produced a diagnosis." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }