{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "h1-00",
   "metadata": {},
   "source": [
    "# Hands-On 1: Βαθμονόμηση συστοιχίας αισθητήρων\n",
    "**Επεξεργασία Σημάτων · ECE-Y523 · Πανεπιστήμιο Πατρών**<br>\n",
    "Κωνσταντίνος Χατζηλυγερούδης · 29 Οκτωβρίου 2026\n",
    "\n",
    "Δέκα αισθητήρες μετρούν την ίδια θερμοκρασία. Δίνουν ακατέργαστες τιμές ADC\n",
    "(counts) με διαφορετικές συστηματικές αποκλίσεις και συσχετισμένο θόρυβο.\n",
    "Βαθμονομήστε τους, δώστε μία θερμοκρασία με αβεβαιότητα και διερευνήστε μια\n",
    "συστηματική διαφωνία ανάμεσα στους αισθητήρες.\n",
    "\n",
    "![Δέκα αισθητήρες σε δεξαμενή αναφοράς με κοινό ADC, ανεμιστήρα και θερμόμετρο αναφοράς· θερμοκρασίες των τριών πειραμάτων](figs/h1_setup_gr.png)\n",
    "\n",
    "Αριστερά φαίνονται η συστοιχία και η δεξαμενή αναφοράς. Οι τρεις διαδρομές\n",
    "θερμοκρασίας δεξιά συνδέουν τη φυσική διάταξη με τα αρχεία δεδομένων.\n",
    "\n",
    "| Ενότητα | Ερώτημα |\n",
    "|---|---|\n",
    "| 1. Βαθμονόμηση | Πώς μετατρέπονται τα counts σε °C; |\n",
    "| 2. Διόρθωση κακής βαθμονόμησης | Τι διορθώνουν η ridge και μια εκ των προτέρων κατανομή από τους υπόλοιπους αισθητήρες; |\n",
    "| 3. Δομή του θορύβου | Ποια σφάλματα είναι κοινά; |\n",
    "| 4. Μία θερμοκρασία | Ποια αβεβαιότητα παραμένει μετά τον συνδυασμό; |\n",
    "| 5. Διάγνωση βλάβης | Είναι συστηματική η διαφωνία και ποιος αισθητήρας την εξηγεί; |\n",
    "| 6. Επικύρωση | Συμφωνεί η τελική αναφορά με το θερμόμετρο αναφοράς; |\n",
    "\n",
    "**Αρχεία.** `calibration.csv`: 400 λήψεις στους 10, 20, 30, 40 °C, δέκα αισθητήρες.\n",
    "`recal_quick.csv`, `recal_full.csv`: δύο έλεγχοι του αισθητήρα 3.\n",
    "`operations.csv`: νέο πείραμα 300 λήψεων. Στα αρχεία της συστοιχίας οι στήλες είναι\n",
    "`snapshot, T_ref, probe_1, …, probe_10`, ενώ στα αρχεία επανελέγχου είναι\n",
    "`T_ref, probe_3`. Κάθε λήψη διαβάζει ταυτόχρονα όλη τη συστοιχία. Στη βαθμονόμηση,\n",
    "μια αναδευόμενη δεξαμενή αναφοράς κρατά όλους τους αισθητήρες στην ακριβώς γνωστή `T_ref`.\n",
    "\n",
    "**Η διαδρομή.** Από τα ακατέργαστα counts φτάνουμε σε βαθμονομημένες\n",
    "θερμοκρασίες, συνδυασμένη εκτίμηση με αβεβαιότητα και διάγνωση συστηματικής διαφωνίας.\n",
    "Η αναφορά λειτουργίας αποκαλύπτεται στην ενότητα 6.\n",
    "\n",
    "**Εκτέλεση του notebook.** Εκτέλεση από `Hands-On 1` με NumPy, SciPy και Matplotlib.\n",
    "Το φυλλάδιο ακολουθεί τις ίδιες έξι ενότητες. Προαιρετικές ασκήσεις:\n",
    "προαιρετικό notebook. Αποδείξεις:\n",
    "μαθηματικά παραρτήματα.\n",
    "\n",
    "**Τι σημαίνουν οι μετρήσεις.** Οι αισθητήρες μοιράζονται την τάση αναφοράς του\n",
    "ADC, την τροφοδοσία και τη ροή αέρα κατά μήκος της συστοιχίας. Επομένως μπορούν\n",
    "να κάνουν παρόμοια σφάλματα την ίδια στιγμή. Μια κοινή ηλεκτρική διαταραχή\n",
    "μετακινεί τις περισσότερες ενδείξεις μαζί· μια αλλαγή ροής αέρα επηρεάζει\n",
    "διαφορετικά τα δύο άκρα. Κάθε αισθητήρας έχει επίσης δικό του θόρυβο. Η\n",
    "βαθμονόμηση διορθώνει μετατόπιση και κέρδος, αλλά δεν αφαιρεί τα τυχαία σφάλματα\n",
    "μελλοντικών ενδείξεων. Δέκα αισθητήρες δεν ισοδυναμούν απαραίτητα με δέκα\n",
    "ανεξάρτητες μετρήσεις.\n",
    "\n",
    "Η δεξαμενή αναφοράς δίνει ακριβώς γνωστή είσοδο στη βαθμονόμηση. Στην κανονική\n",
    "λειτουργία αυτή η αναφορά απουσιάζει, άρα θερμοκρασία και διάγνωση πρέπει να\n",
    "προκύπτουν από τη συστοιχία. Κρατώντας αρχικά κρυφή την αναφορά επικύρωσης,\n",
    "διακρίνουμε όσα συμπεραίνονται από τους αισθητήρες από όσα ελέγχονται μόνο\n",
    "όταν γνωρίζουμε την αληθινή θερμοκρασία.\n",
    "\n",
    "**Σημειογραφία.** Οι πίνακες έχουν έντονα κεφαλαία σύμβολα, όπως\n",
    "$\\boldsymbol{X}$ και $\\boldsymbol{S}$· τα διανύσματα έντονα πεζά, όπως\n",
    "$\\boldsymbol{w}$ και $\\boldsymbol{r}$. Οι βαθμωτές ποσότητες $T_k$, $a_j$, $b_j$\n",
    "παραμένουν απλές. Διαστάσεις και μονάδες συνδέουν τους τύπους με τις πράξεις NumPy."
   ]
  },
  {
   "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(\"el\")\n",
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-02",
   "metadata": {},
   "source": [
    "## 1. Βαθμονόμηση: counts → °C\n",
    "*Διαλέξεις 1 και 5*\n",
    "\n",
    "Η θερμοκρασία αναφοράς είναι ακριβώς γνωστή· ο θόρυβος βρίσκεται στα counts.\n",
    "Προσαρμόζουμε $y_{kj}=a_j+b_j T_k+n_{kj}$, με $\\boldsymbol{X}_{k,:}=[1,T_k]$.\n",
    "Ο ίδιος πίνακας σχεδιασμού $N\\times2$ εξυπηρετεί όλους τους αισθητήρες:\n",
    "$\\boldsymbol{Y}=\\boldsymbol{X}\\boldsymbol{\\Theta}+\\text{θόρυβος}$, όπου ο $\\boldsymbol{\\Theta}$ είναι $2\\times C$.\n",
    "Χρησιμοποιούμε `lstsq(X, Y)` και αντιστρέφουμε κάθε προσαρμοσμένη ευθεία.\n",
    "\n",
    "**Από έναν αισθητήρα στην εξίσωση πινάκων.** Για τον αισθητήρα $j$, συγκεντρώνουμε\n",
    "τα counts στο $\\boldsymbol{y}_j\\in\\mathbb{R}^N$ και τους δύο συντελεστές στο\n",
    "$\\boldsymbol{\\theta}_j=[a_j,b_j]^{\\top}$. Κάθε γραμμή του πίνακα σχεδιασμού έχει\n",
    "τη σταθερά για τη μετατόπιση και τη θερμοκρασία αναφοράς για το κέρδος:\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",
    "Όλοι οι αισθητήρες βλέπουν τις ίδιες θερμοκρασίες, οπότε τα δεδομένα τους\n",
    "γίνονται στήλες του $\\boldsymbol{Y}$ και οι συντελεστές στήλες του\n",
    "$\\boldsymbol{\\Theta}$. Η `lstsq` δέχεται πίνακα δεξιών μελών. Λύνει τις δέκα\n",
    "προσαρμογές μαζί χωρίς να επιβάλλει ίσα κέρδη ή μετατοπίσεις. Η γνωστή έκφραση\n",
    "$(\\boldsymbol{X}^{\\top}\\boldsymbol{X})^{-1}\\boldsymbol{X}^{\\top}\\boldsymbol{Y}$\n",
    "περιγράφει τη λύση όταν οι στήλες είναι ανεξάρτητες. Η `lstsq` την υπολογίζει\n",
    "χωρίς να σχηματίζει ρητά αυτόν τον αντίστροφο.\n",
    "\n",
    "**Γιατί αυτή η κατεύθυνση προσαρμογής;** Τα ελάχιστα τετράγωνα θεωρούν γνωστό\n",
    "τον πίνακα σχεδιασμού και τοποθετούν τον θόρυβο στην έξοδο. Η ακριβής\n",
    "θερμοκρασία της δεξαμενής ανήκει στον $\\boldsymbol{X}$ και τα θορυβώδη counts\n",
    "στον $\\boldsymbol{Y}$. Η απευθείας προσαρμογή θερμοκρασίας ως προς counts\n",
    "αλλάζει την υπόθεση. Προσαρμόζουμε πρώτα τη φυσική απόκριση και την\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": [
    "**Ανάγνωση του αποτελέσματος.** Η πρώτη γραμμή του `Theta` περιέχει μετατοπίσεις\n",
    "σε counts και η δεύτερη κέρδη σε counts/°C. Η NumPy τις αποδίδει στα `a`, `b`.\n",
    "Αφαιρώντας τα προσαρμοσμένα counts παίρνουμε τον πίνακα υπολοίπων\n",
    "$\\boldsymbol{R}$, διαστάσεων $N\\times C$. Η διαίρεση των στηλών του με τα\n",
    "αντίστοιχα κέρδη δίνει τον $\\boldsymbol{E}$ σε °C. Το broadcasting εφαρμόζει\n",
    "κάθε κέρδος στον σωστό αισθητήρα.\n",
    "\n",
    "Το αριστερό διάγραμμα ελέγχει τη γραμμική περιγραφή της απόκρισης. Το δεξί\n",
    "εξετάζει αν απομένει συστηματική δομή στις τέσσερις θερμοκρασίες. Καμπυλότητα ή\n",
    "αποκλίσεις που εξαρτώνται από τη δεξαμενή θα αμφισβητούσαν το μοντέλο. Η\n",
    "τυπική απόκλιση υπολοίπων διαιρεί με $N-p$, όπου $p=2$, επειδή εκτιμήθηκαν δύο\n",
    "συντελεστές. Εκτιμά τον θόρυβο, όχι την αβεβαιότητα του προσαρμοσμένου κέρδους.\n",
    "Η σύγκριση μίας λήψης δείχνει την αφαίρεση συστηματικού σφάλματος· η διασπορά\n",
    "που απομένει δείχνει τον θόρυβο που δεν εξαλείφει η βαθμονόμηση.\n",
    "\n",
    "**Συζήτηση.** Πώς διαφέρουν η ονομαστική $(a,b)=(200,40)$ και η προσαρμοσμένη\n",
    "βαθμονόμηση για την ίδια λήψη; Τι αλλάζει: συστηματική απόκλιση, διασπορά ή και τα δύο;\n",
    "\n",
    "**Σημείο ελέγχου.** Τα counts έγιναν °C. Ο θόρυβος κάθε αισθητήρα είναι περίπου\n",
    "0.2 °C. Μικρές τάσεις στα υπόλοιπα στηρίζουν το μοντέλο· δεν αποδεικνύουν κανονικότητα.\n",
    "Τα συσχετισμένα σφάλματα θα μας απασχολήσουν στην ενότητα 3."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-05",
   "metadata": {},
   "source": [
    "## 2. Διόρθωση κακής βαθμονόμησης\n",
    "*Διαλέξεις 2, 5 και 6*\n",
    "\n",
    "Η σύντομη επαναβαθμονόμηση του αισθητήρα 3 γίνεται μόνο στους 25.0 και 25.2 °C. Οι θερμοκρασίες\n",
    "αυτές δίνουν ελάχιστο μοχλοβραχίονα για την εκτίμηση του κέρδους. Συγκρίνουμε:\n",
    "\n",
    "- **LS:** ελαχιστοποίηση του $\\|\\boldsymbol{y}-\\boldsymbol{X}_q\\boldsymbol{\\theta}\\|^2$.\n",
    "- **Απλή ridge:** πρόσθεση του $\\lambda\\|\\boldsymbol{\\theta}\\|^2$, με $\\lambda=1$.\n",
    "- **MAP:** κανονική εκ των προτέρων κατανομή από τους άλλους εννέα αισθητήρες,\n",
    "  με διαφορετική κλίμακα για μετατόπιση και κέρδος. Για τυπική απόκλιση θορύβου\n",
    "  $\\sigma$ και εκ των προτέρων ακρίβεια $\\boldsymbol{\\Lambda}_0$, λύνουμε\n",
    "  $(\\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",
    "Είναι μια ελεγχόμενη σύγκριση. Στα επόμενα βήματα διατηρούμε την πλήρη\n",
    "βαθμονόμηση των τεσσάρων θερμοκρασιών από την ενότητα 1.\n",
    "\n",
    "**Ο μοχλοβραχίονας που λείπει.** Οι δύο στήλες του $\\boldsymbol{X}_q$ είναι\n",
    "σχεδόν παράλληλες: κοντά στους 25 °C η στήλη θερμοκρασίας είναι περίπου 25 φορές\n",
    "η στήλη μονάδων. Αύξηση του $b$ και μείωση του $a$ μπορούν να αφήσουν σχεδόν\n",
    "ίδιο το $a+25b$. Έτσι πολύ διαφορετικές ευθείες προσαρμόζονται στο στενό εύρος.\n",
    "Για ανεξάρτητο θόρυβο counts με διακύμανση $\\sigma^2$,\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",
    "Η πλήρης βαθμονόμηση έχει $s_{TT}=50\\,000$ °C², ενώ η σύντομη επαναβαθμονόμηση μόνο\n",
    "$0.2$ °C². Με παρόμοιο θόρυβο η τυπική απόκλιση του κέρδους γίνεται περίπου\n",
    "500 φορές μεγαλύτερη. Μια δεύτερη θερμοκρασία βοηθά ουσιαστικά όταν απέχει\n",
    "αρκετά από την πρώτη.\n",
    "\n",
    "**Τι υποθέτει η κανονικοποίηση.** Η απλή ridge τιμωρεί εξίσου το αριθμητικό\n",
    "μέγεθος των δύο συντελεστών και συρρικνώνει προς το μηδέν. Όμως μετατόπιση σε\n",
    "counts και κέρδος σε counts/°C έχουν διαφορετική κλίμακα και φυσική σημασία.\n",
    "Η MAP χρησιμοποιεί τους άλλους αισθητήρες για μια εύλογη μέση τιμή και\n",
    "ξεχωριστές διασπορές. Ο `Lambda0` είναι η αντίστροφη εκ των προτέρων\n",
    "συνδιακύμανση. Μεγαλύτερες διακυμάνσεις σημαίνουν μικρότερη ακρίβεια, άρα\n",
    "περισσότερη επιρροή των δεδομένων. Πρόκειται για προσεγγιστική εμπειρική εκ των προτέρων\n",
    "κατανομή από εννέα παρόμοιους αισθητήρες, όχι βεβαιότητα για τον επισκευασμένο. Η κλίμακα\n",
    "θορύβου counts προέρχεται από την πλήρη βαθμονόμηση."
   ]
  },
  {
   "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": [
    "**Ανάγνωση της σύγκρισης.** Η μεγέθυνση κοντά στους 25 °C δείχνει γιατί τα\n",
    "υπόλοιπα της σύντομης επαναβαθμονόμησης μπορεί να παραπλανούν. Και οι τρεις μέθοδοι\n",
    "εξηγούν τα διαθέσιμα counts, ενώ το πλήρες εύρος αποκαλύπτει διαφορετικές\n",
    "κλίσεις. Στους 40 °C η LS διαφέρει από την πλήρη επαναβαθμονόμηση κατά περίπου +11 °C\n",
    "σε ισοδύναμες μονάδες θερμοκρασίας· η απλή ridge περιορίζει τη διαφορά σε\n",
    "+2.4 °C και η MAP με εκ των προτέρων κατανομή από τους υπόλοιπους αισθητήρες σε −0.27 °C. Οι διαφορές αφορούν προσαρμοσμένες\n",
    "ευθείες απόκρισης: διαφορά counts διαιρεμένη με το κέρδος της πλήρους επαναβαθμονόμησης.\n",
    "\n",
    "Η πλήρης επαναβαθμονόμηση είναι μέτρηση σύγκρισης, όχι ακριβώς γνωστή αλήθεια παραμέτρων.\n",
    "Αλλαγή του `prior_scale` από 1 σε 5 αυξάνει κάθε εκ των προτέρων διακύμανση\n",
    "κατά 25. Η MAP βασίζεται περισσότερο στα ασθενή δεδομένα και κινείται προς\n",
    "την LS. Η σύγκριση εξηγεί την κανονικοποίηση· δεν αντικαθιστά τους πλήρεις\n",
    "συντελεστές της συστοιχίας στα επόμενα βήματα.\n",
    "\n",
    "**Συζήτηση.** Τι αλλάζει όταν οι δύο εκ των προτέρων τυπικές αποκλίσεις\n",
    "πολλαπλασιάζονται επί πέντε; Πώς θα βοηθούσαν πιο απομακρυσμένες θερμοκρασίες αναφοράς;\n",
    "\n",
    "**Σημείο ελέγχου.** Μια καλή τοπική προσαρμογή μπορεί να δίνει πολύ κακή\n",
    "εξωπαρεμβολή. Η ridge σταθεροποιεί τη λύση· η επιλογή της εκ των\n",
    "προτέρων κατανομής καθορίζει πού θα σταθεροποιηθεί. Τα παραρτήματα και το\n",
    "προαιρετικό notebook αναπτύσσουν την SVD και τη σχέση μεροληψίας–διακύμανσης."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-08",
   "metadata": {},
   "source": [
    "## 3. Δομή του θορύβου\n",
    "*Διαλέξεις 4 και 5*\n",
    "\n",
    "Η προσαρμογή χρησιμοποίησε $p=2$ παραμέτρους ανά αισθητήρα. Εκτιμούμε τη\n",
    "συνδιακύμανση σε °C² ως $\\boldsymbol{S}=\\boldsymbol{E}^{\\top} \\boldsymbol{E}/(N-p)$. Τα ιδιοδιανύσματα περιγράφουν μοτίβα\n",
    "θορύβου και οι ιδιοτιμές δίνουν τις διακυμάνσεις τους. Παραγοντοποιούμε\n",
    "$\\boldsymbol{S}=\\boldsymbol{L}\\boldsymbol{L}^{\\top}$ και λύνουμε $\\boldsymbol{L}\\boldsymbol{q}=\\boldsymbol{e}$ για να λευκοποιήσουμε ένα διάνυσμα σφάλματος.\n",
    "\n",
    "**Τι σημαίνει ένα στοιχείο συνδιακύμανσης.** Το διαγώνιο $S_{jj}$ εκτιμά τη\n",
    "διακύμανση του αισθητήρα $j$. Το εκτός διαγωνίου $S_{ij}$ καταγράφει πώς\n",
    "μεταβάλλονται μαζί τα σφάλματα των $i,j$ στην *ίδια λήψη*. Η συσχέτιση αφαιρεί\n",
    "τις διαφορετικές κλίμακες των δύο αισθητήρων και διευκολύνει τη σύγκριση ζευγών.\n",
    "Αυτή η εξάρτηση αφορά αισθητήρες· η ανεξαρτησία διαφορετικών λήψεων παραμένει\n",
    "χωριστή υπόθεση του μοντέλου.\n",
    "\n",
    "Η διάσπαση $\\boldsymbol{S}=\\boldsymbol{U}\\boldsymbol{\\Lambda}\\boldsymbol{U}^{\\top}$\n",
    "δίνει ορθοκανονικές διευθύνσεις $\\boldsymbol{u}_i$. Η βαθμωτή συντεταγμένη\n",
    "$p_i=\\boldsymbol{u}_i^{\\top}\\boldsymbol{e}$ έχει διακύμανση $\\lambda_i$ στο\n",
    "εκτιμημένο μοντέλο. Μεγάλη ιδιοτιμή σημαίνει έντονη μεταβολή κατά το συγκεκριμένο\n",
    "μοτίβο. Τα ιδιοδιανύσματα περιγράφουν μοτίβα, όχι μεμονωμένους αισθητήρες.\n",
    "Αντιστροφή προσήμου δεν αλλάζει τη συνδιακύμανση ή την ερμηνεία τους.\n",
    "\n",
    "**Τι κάνει η επίλυση.** Για ένα διάνυσμα στήλη, η λευκοποίηση λύνει\n",
    "$\\boldsymbol{L}\\boldsymbol{q}=\\boldsymbol{e}$. Επειδή τα σφάλματα αποθηκεύονται\n",
    "σε γραμμές, η `E.T` τα μετατρέπει σε στήλες και η τελική αναστροφή αποκαθιστά\n",
    "τη διάταξη $N\\times C$. Ο πίνακας $\\boldsymbol{Q}$ περιέχει τις νέες\n",
    "συντεταγμένες ανά γραμμή. Οι συντεταγμένες αναμειγνύουν αισθητήρες· μια μεγάλη\n",
    "τιμή δεν ονομάζει απαραίτητα έναν ελαττωματικό αισθητήρα."
   ]
  },
  {
   "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": [
    "**Ανάγνωση των σχημάτων.** Ισχυρή θετική συσχέτιση σημαίνει ότι μια υψηλή\n",
    "ένδειξη συχνά συνοδεύεται από υψηλή ένδειξη άλλου αισθητήρα. Το πρώτο μοτίβο\n",
    "είναι σχεδόν ομοιόμορφο και το δεύτερο μεταβάλλεται κατά μήκος της συστοιχίας.\n",
    "Τα μεγάλα κλάσματα διακύμανσης εξηγούν γιατί το νέφος δύο σφαλμάτων είναι\n",
    "στενό και κεκλιμένο. Η λευκοποίηση αφαιρεί κλίμακα και κλίση σύμφωνα με το\n",
    "εκτιμημένο μοντέλο.\n",
    "\n",
    "Ο έλεγχος ταυτότητας γίνεται σκόπιμα στα δεδομένα βαθμονόμησης για να φανεί\n",
    "ο περιορισμός του. Ο ίδιος $\\boldsymbol{E}$ κατασκεύασε τον $\\boldsymbol{S}$,\n",
    "άρα η ταυτότητα προκύπτει αλγεβρικά ακόμη και με μη κανονικά δεδομένα. Ο\n",
    "έλεγχος μοντέλου χρειάζεται νέα δεδομένα. Στη συνέχεια, τα ανεξάρτητα υγιή\n",
    "σφάλματα επικύρωσης θα ελέγξουν αν η κλίμακα και η δομή μεταφέρονται στη λειτουργία.\n",
    "\n",
    "**Συζήτηση.** Ποιο ιδιοδιάνυσμα είναι σχεδόν ομοιόμορφο και ποιο μεταβάλλεται\n",
    "από μπροστά προς τα πίσω; Ποια φυσική διαταραχή εξηγεί το καθένα; Το πρόσημο είναι αυθαίρετο.\n",
    "\n",
    "**Σημείο ελέγχου.** Δύο μοτίβα συγκεντρώνουν περίπου το 94% της διακύμανσης.\n",
    "Η λευκοποίηση των ίδιων δεδομένων δίνει την ταυτότητα από κατασκευή. Ανεξάρτητος\n",
    "έλεγχος γίνεται στην ενότητα 6· η ταυτότητα από μόνη της δεν επικυρώνει το μοντέλο."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-11",
   "metadata": {},
   "source": [
    "## 4. Μία θερμοκρασία με αβεβαιότητα\n",
    "*Διάλεξη 5*\n",
    "\n",
    "Για ενδείξεις $\\boldsymbol{t}=T\\boldsymbol{1}+\\boldsymbol{e}$, τα σταθμισμένα ελάχιστα τετράγωνα δίνουν\n",
    "$\\boldsymbol{w}=\\boldsymbol{S}^{-1}\\boldsymbol{1}/(\\boldsymbol{1}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{1})$ και\n",
    "$\\hat T=\\boldsymbol{w}^{\\top}\\boldsymbol{t}$, με $\\mathrm{sd}(\\hat T)=(\\boldsymbol{1}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{1})^{-1/2}$.\n",
    "Η συνάρτηση επιστρέφει και την τετραγωνική απόσταση του υπολοίπου\n",
    "$\\boldsymbol{r}=\\boldsymbol{t}-\\hat T\\boldsymbol{1}$, που θα χρησιμοποιηθεί στην ενότητα 5. Αρχικά κρατάμε όλους τους αισθητήρες.\n",
    "\n",
    "**Το πρόβλημα εκτίμησης.** Κάθε βαθμονομημένος αισθητήρας παρατηρεί την ίδια\n",
    "άγνωστη βαθμωτή θερμοκρασία. Ο πίνακας σχεδιασμού αυτού του νέου προβλήματος\n",
    "είναι η στήλη $\\boldsymbol{1}$, όχι ο πίνακας δύο στηλών της βαθμονόμησης.\n",
    "Η WLS επιλέγει την κοινή στάθμη που ελαχιστοποιεί\n",
    "\n",
    "$$\n",
    "J(T)=(\\boldsymbol{t}-T\\boldsymbol{1})^{\\top}\\boldsymbol{S}^{-1}\n",
    "     (\\boldsymbol{t}-T\\boldsymbol{1}).\n",
    "$$\n",
    "\n",
    "Η παραγώγιση δίνει\n",
    "$\\hat T=\\boldsymbol{1}^{\\top}\\boldsymbol{S}^{-1}\\boldsymbol{t}/\n",
    "(\\boldsymbol{1}^{\\top}\\boldsymbol{S}^{-1}\\boldsymbol{1})$.\n",
    "Τα βάρη αθροίζουν στη μονάδα, άρα αναπαράγεται ακριβώς μια κοινή θερμοκρασία\n",
    "χωρίς θόρυβο. Το `Si1` υπολογίζει $\\boldsymbol{S}^{-1}\\boldsymbol{1}$ με\n",
    "επίλυση· η διαίρεση με το άθροισμα κανονικοποιεί τα βάρη.\n",
    "\n",
    "**Ο υπολογισμός αβεβαιότητας.** Για σταθερά βάρη ισχύει\n",
    "$\\operatorname{Var}(\\boldsymbol{w}^{\\top}\\boldsymbol{e})=\n",
    "\\boldsymbol{w}^{\\top}\\boldsymbol{S}\\boldsymbol{w}$.\n",
    "Ο απλός μέσος έχει διακύμανση $\\boldsymbol{1}^{\\top}\\boldsymbol{S}\\boldsymbol{1}/C^2$.\n",
    "Η υπόθεση ανεξαρτησίας αφαιρεί όλα τα εκτός διαγωνίου στοιχεία. Ένα κοινό σφάλμα\n",
    "εμφανίζεται σε κάθε ένδειξη και παραμένει σε κάθε αμερόληπτο συνδυασμό.\n",
    "Περισσότεροι αισθητήρες μειώνουν το ανεξάρτητο μέρος, όχι το κοινό."
   ]
  },
  {
   "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": [
    "**Ανάγνωση βαρών και στηλών.** Αρνητικά βάρη επιτρέπονται: η WLS είναι\n",
    "γραμμικός συνδυασμός, όχι ψηφοφορία με θετικά ποσοστά. Ένας αισθητήρας μπορεί\n",
    "να βοηθά στην αναίρεση διαταραχής που βλέπουν οι υπόλοιποι. Οι ακριβείς τιμές\n",
    "εξαρτώνται όμως και από την εκτιμημένη συνδιακύμανση· δεν αποτελούν από μόνες\n",
    "τους φυσική κατάταξη ποιότητας αισθητήρων.\n",
    "\n",
    "Οι στήλες συγκρίνουν *τυπικές αποκλίσεις*, όχι εύρη διαστημάτων 95%. Στο\n",
    "κανονικό μοντέλο ένα ονομαστικό διάστημα 95% είναι περίπου $\\hat T\\pm1.96s$.\n",
    "Εδώ WLS και σωστά αξιολογημένος απλός μέσος έχουν παρόμοια αβεβαιότητα. Το\n",
    "μεγάλο σφάλμα θα ήταν να συνδεθεί ο μέσος με περιθώριο ανεξάρτητων αισθητήρων.\n",
    "Το υπόλοιπο που επιστρέφει η `combine` αφαιρεί την εκτιμημένη κοινή θερμοκρασία\n",
    "και θα αποκαλύψει διαφωνίες μεταξύ αισθητήρων στην επόμενη ενότητα.\n",
    "\n",
    "**Συζήτηση.** Πώς αλλάζει η αβεβαιότητα αν μηδενιστούν τα εκτός διαγωνίου\n",
    "στοιχεία του $\\boldsymbol{S}$; Ποια υπόθεση αλλάζει;\n",
    "\n",
    "**Σημείο ελέγχου.** Ο κοινός θόρυβος παραμένει μετά τον μέσο όρο. Η WLS και ο\n",
    "απλός μέσος έχουν εδώ αβεβαιότητα κοντά στους 0.19 °C· η υπόθεση ανεξαρτησίας\n",
    "δίνει περίπου 0.07 °C. Αρνητικά βάρη μπορούν να αναιρούν συσχετισμένες διαταραχές.\n",
    "Οι αβεβαιότητες βασίζονται στο μοντέλο· απομένουν διάγνωση και επικύρωση."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-14",
   "metadata": {},
   "source": [
    "## 5. Διάγνωση συστηματικής διαφωνίας\n",
    "*Διάλεξη 4 και σύντομη εισαγωγή στη χ²*\n",
    "\n",
    "Για ανεξάρτητες Gaussian μεταβλητές $z_i$ με μηδενική μέση τιμή και μοναδιαία\n",
    "διακύμανση, $z_i\\sim\\mathcal N(0,1)$, το άθροισμα\n",
    "$\\sum_{i=1}^{\\nu}z_i^2$ ακολουθεί $\\chi^2_\\nu$. Το λευκοποιημένο υπόλοιπο έχει\n",
    "**$C-1$ βαθμούς ελευθερίας**: η εκτίμηση μίας θερμοκρασίας αφαιρεί μία διεύθυνση.\n",
    "Συγκρίνουμε το $d^2=\\boldsymbol{r}^{\\top} \\boldsymbol{S}^{-1}\\boldsymbol{r}$ με `stats.chi2.ppf(0.99, C-1)`.\n",
    "Με εκτιμημένο $\\boldsymbol{S}$, το κατώφλι δίνει προσεγγιστικά 1% συναγερμούς ανά λήψη.\n",
    "\n",
    "**Ένας συναγερμός δεν καταδικάζει έναν αισθητήρα.** Ακόμη και ένα υγιές πείραμα\n",
    "300 λήψεων αναμένεται να δώσει περίπου τρεις συναγερμούς. Πρώτα ελέγχουμε αν το\n",
    "*πλήθος* είναι ασυνήθιστα μεγάλο. Η `binom.sf(m-1, n, alpha)` δίνει την πιθανότητα\n",
    "τουλάχιστον `m` συναγερμών σε `n` ανεξάρτητες υγιείς λήψεις με ρυθμό `alpha`.\n",
    "\n",
    "Αν υπάρχουν υπερβολικά πολλοί συναγερμοί, αφαιρούμε κάθε αισθητήρα και\n",
    "επανεκτιμούμε τη θερμοκρασία (τώρα $C-2$ βαθμοί ελευθερίας). Για την επίδειξη με\n",
    "**μία μόνιμη βλάβη**, αποκλείουμε έναν αισθητήρα μόνο αν είναι ο\n",
    "μοναδικός που εξηγεί τους περισσότερους συναγερμούς και η υπόλοιπη συστοιχία\n",
    "παύει να δίνει υπερβολικά πολλούς. Αλλιώς διατηρούμε τους αισθητήρες και δηλώνουμε\n",
    "αβεβαιότητα στη διάγνωση. Ο κανόνας αναλύει όλο το πείραμα μαζί· δεν εγγυάται\n",
    "τη διάγνωση. Οι έλεγχοι προϋποθέτουν κανονικό θόρυβο, σταθερή συνδιακύμανση και\n",
    "ανεξάρτητες λήψεις.\n",
    "\n",
    "**Γιατί δεν χρειάζεται αναφορά.** Η αφαίρεση του $\\hat T\\boldsymbol{1}$\n",
    "απαλείφει την κοινή στάθμη, διατηρώντας τις διαφωνίες μεταξύ αισθητήρων.\n",
    "Στις λευκοποιημένες συντεταγμένες η εκτίμηση στάθμης αφαιρεί μία διεύθυνση,\n",
    "άρα η τετραγωνική απόσταση έχει $C-1$ βαθμούς ελευθερίας. Για $C=10$, το\n",
    "ονομαστικό κατώφλι 99% είναι 21.67. Μεγάλη απόσταση σημαίνει διαφωνία σε\n",
    "διεύθυνση με μικρό αναμενόμενο θόρυβο. Μια κοινή κίνηση συμβατή με αληθινή\n",
    "αλλαγή θερμοκρασίας απορροφάται κυρίως στο $\\hat T$.\n",
    "\n",
    "**Από τον συναγερμό στη διάγνωση.** Το `alarm_p` είναι πιθανότητα ουράς για\n",
    "το *πλήθος συναγερμών υπό το υγιές μοντέλο*. Δεν είναι πιθανότητα βλάβης\n",
    "αισθητήρα. Αναμενόμενοι μεμονωμένοι συναγερμοί δεν δικαιολογούν αποκλεισμό.\n",
    "Αν το πλήθος είναι ασυνήθιστα υψηλό, οι έλεγχοι αφαίρεσης εξετάζουν ποιος\n",
    "αισθητήρας εξηγεί τη διαφωνία. Χρησιμοποιείται ο υποπίνακας συνδιακύμανσης\n",
    "των αισθητήρων που παραμένουν και επανεκτιμάται η κοινή θερμοκρασία. Ο\n",
    "υποπίνακας του πλήρους πίνακα ακρίβειας θα ήταν διαφορετική πράξη.\n",
    "\n",
    "Το `blame[j]` είναι το κλάσμα αρχικών συναγερμών που παύει χωρίς τον αισθητήρα\n",
    "$j$. Ένας μοναδικός ισχυρός υποψήφιος αποκλείεται μόνο αν και τα υπόλοιπα\n",
    "δεδομένα περάσουν τον έλεγχο πλήθους. Η δεύτερη συνθήκη αποτρέπει τη σύγχυση\n",
    "μερικής εξήγησης με πλήρη διόρθωση. Πολλές βλάβες, αλλαγή μοντέλου θορύβου\n",
    "ή αμφίσημες βαθμολογίες μπορεί να αφήσουν αβέβαιη τη διάγνωση."
   ]
  },
  {
   "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": [
    "**Ανάγνωση χρονικής εξέλιξης και αφαιρέσεων.** Το αριστερό διάγραμμα δίνει\n",
    "την απόσταση κάθε λήψης· τα κόκκινα σημεία υπερβαίνουν το κατώφλι. Η ομάδα\n",
    "όψιμων συναγερμών είναι πιο πληροφοριακή από μία πρώιμη υπέρβαση. Το δεξί\n",
    "διάγραμμα εξετάζει ποια αφαίρεση παύει τη διαφωνία. Μια υψηλή στήλη αποκτά\n",
    "νόημα μαζί με το υπερβολικό πλήθος και τον έλεγχο της υπόλοιπης συστοιχίας,\n",
    "όχι ανεξάρτητα από αυτά.\n",
    "\n",
    "Οι πρώτες 100 λήψεις δείχνουν τη διάκριση: ένα υγιές τμήμα μπορεί να έχει\n",
    "υπερβάσεις χωρίς στοιχεία μόνιμης βλάβης. Οι σημασμένες λήψεις διατηρούνται\n",
    "στην έξοδο ακόμη και αν δεν αποκλειστεί αισθητήρας. Η χρονική εξέλιξη\n",
    "υποδεικνύει πιθανή αλλαγή συμπεριφοράς, αλλά ο υπολογισμός πάνω σε όλο το πείραμα\n",
    "δεν εκτιμά χρόνο έναρξης ούτε υλοποιεί κανόνα διακοπής σε πραγματικό χρόνο.\n",
    "\n",
    "**Συζήτηση.** Πού γίνεται συστηματική η διαφωνία στο διάγραμμα; Δικαιολογεί το\n",
    "`diagnose(T_probes[:100], S)` αποκλεισμό αισθητήρα; Γιατί δεν αρκούν μερικοί συναγερμοί;\n",
    "\n",
    "**Σημείο ελέγχου.** Ο ανιχνευτής δεν χρησιμοποίησε θερμοκρασίες αναφοράς,\n",
    "δοσμένο χρόνο βλάβης ή αριθμό αισθητήρα. Το διάγραμμα στηρίζει τη διερεύνηση\n",
    "μιας αλλαγής· δεν υλοποιεί εκτιμητή χρονικής στιγμής αλλαγής. Η διάγνωση\n",
    "προϋποθέτει μία μόνιμη βλάβη και μπορεί να παραμείνει αβέβαιη."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-17",
   "metadata": {},
   "source": [
    "## 6. Επικύρωση του αποτελέσματος\n",
    "\n",
    "Συνδυάστε ξανά τους αισθητήρες που διατηρούνται. **Τώρα** αποκαλύπτουμε τη\n",
    "θερμοκρασία αναφοράς και ελέγχουμε τα περιθώρια σφάλματος. Ελέγχουμε επίσης\n",
    "τη λευκοποίηση σε νέα σφάλματα. Μόνο για την τελική επικύρωση, η γνωστή αλήθεια\n",
    "των συνθετικών δεδομένων είναι ότι ο αισθητήρας 7 αποκτά μετατόπιση +0.30 °C\n",
    "από τη λήψη 200. Χρησιμοποιούμε τις πρώτες 200 λήψεις για έλεγχο υγιούς\n",
    "λειτουργίας· οι ετικέτες αυτές δεν συμμετείχαν στη διάγνωση.\n",
    "\n",
    "**Τρεις διαφορετικοί έλεγχοι.** Πρώτα συγκρίνουμε την τελική θερμοκρασία με\n",
    "την αναφορά: μεροληψία και διασπορά περιγράφουν τα σφάλματα της εκτίμησης.\n",
    "Έπειτα ελέγχουμε την κάλυψη. Για ονομαστικό επίπεδο $q$ χρησιμοποιούμε\n",
    "$z_q=\\Phi^{-1}((1+q)/2)$ και μετράμε τις λήψεις που ικανοποιούν\n",
    "$|\\hat T_k-T_{\\mathrm{ref},k}|\\le z_qs$. Σωστά βαθμονομημένα διαστήματα\n",
    "ακολουθούν τη διαγώνιο στο δεξί διάγραμμα. Τέλος λευκοποιούμε νέα υγιή σφάλματα\n",
    "με τον *παλιό* παράγοντα βαθμονόμησης $\\boldsymbol{L}$. Η συνδιακύμανσή τους\n",
    "πρέπει να είναι κοντά στην ταυτότητα αν μεταφέρεται η δομή θορύβου. Νέα\n",
    "εκτίμηση μετασχηματισμού στα ίδια σφάλματα θα ακύρωνε αυτόν τον έλεγχο.\n",
    "\n",
    "Οι έλεγχοι απαντούν διαφορετικά ερωτήματα. Εύλογα μέσα σφάλματα δεν\n",
    "εγγυώνται σωστή αβεβαιότητα, ενώ σωστή μέση κάλυψη δεν ταυτοποιεί βλάβη.\n",
    "Ο γνωστός χρόνος και ο δείκτης βλάβης χρησιμοποιούνται τώρα για διαχωρισμό\n",
    "υγιών και ελαττωματικών λήψεων επικύρωσης. Δεν ήταν διαθέσιμοι στη διάγνωση.\n",
    "Ο τελικός συνδυασμός χρησιμοποιεί τον υποπίνακα των διατηρούμενων αισθητήρων\n",
    "και επανυπολογίζει τα βάρη· δεν διαγράφει απλώς ένα παλιό βάρος."
   ]
  },
  {
   "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": [
    "**Τελική συζήτηση.**\n",
    "1. Τι στηρίζει το μοντέλο βαθμονόμησης;\n",
    "2. Γιατί ήταν αναξιόπιστο το κέρδος της σύντομης επαναβαθμονόμησης και τι πρόσφερε η εκ των προτέρων κατανομή;\n",
    "3. Τι έδειξαν το πλήθος συναγερμών, οι αφαιρέσεις αισθητήρων και η χρονική εξέλιξη;\n",
    "4. Ποια αβεβαιότητα θερμοκρασίας δίνετε και ποια κάλυψη πέτυχε;\n",
    "5. Ποιες υποθέσεις μπορεί να αποτύχουν σε νέο πείραμα;\n",
    "\n",
    "**Αναμενόμενη επικύρωση.** Όλο το πείραμα λειτουργίας στηρίζει τον αισθητήρα 7 ως ύποπτο.\n",
    "Ο έλεγχος χωρίς αναφορά σημαίνει περίπου 1.5% των υγιών λήψεων και 70% μετά τη\n",
    "βλάβη. Με αποκλεισμό του αισθητήρα 7, η τυπική απόκλιση είναι περίπου 0.190 °C\n",
    "και η κάλυψη 95% περίπου 92%. Η μικρή υποκάλυψη και οι αποκλίσεις στη λευκοποίηση\n",
    "νέων δεδομένων θυμίζουν ότι ο $\\boldsymbol{S}$ είναι εκτιμημένος. Η αβεβαιότητα των συντελεστών\n",
    "βαθμονόμησης είναι μικρή για την πλήρη βαθμονόμηση (παράρτημα F) και εδώ παραλείπεται.\n",
    "Ο αποκλεισμός επιλέχθηκε από το ίδιο πείραμα· η επικύρωση είναι έλεγχος λυμένης\n",
    "εφαρμογής, όχι ανεξάρτητη αξιολόγηση της επίδοσης ταξινόμησης βλαβών.\n",
    "\n",
    "**Προαιρετική συνέχεια.** Το προαιρετικό notebook καλύπτει στοίβαξη Kronecker,\n",
    "αβεβαιότητα συντελεστών, επαναλαμβανόμενους ελέγχους, SVD/ridge, καμπύλες\n",
    "εκ των υστέρων κατανομής, γεωμετρία χ² και ισχύ ανίχνευσης. Δοκιμάστε λιγότερες\n",
    "λήψεις βαθμονόμησης, άλλη εκ των προτέρων κατανομή ή πολλές βλάβες. Κρατήστε\n",
    "κρυφή τη θερμοκρασία αναφοράς μέχρι να ολοκληρώσετε τη διάγνωση."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
