{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "h1-ext-00",
   "metadata": {},
   "source": [
    "# Hands-On 1: προαιρετικές ασκήσεις\n",
    "Εκτελέστε τις μετά τις έξι ενότητες του hands-on. Οι αριθμητικές αυτές ασκήσεις\n",
    "συμπληρώνουν τα μαθηματικά παραρτήματα και την κύρια συνεδρία. Κάθε άσκηση\n",
    "εξετάζει ένα πρόσθετο ερώτημα. Η αρχικοποίηση επανυπολογίζει την πλήρη βαθμονόμηση και τη\n",
    "σύγκριση της σύντομης επαναβαθμονόμησης, ώστε το notebook να εκτελείται αυτόνομα."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-01",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import session_plotting as plots\n",
    "from session_plotting import tr\n",
    "plots.configure(\"el\")\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "\n",
    "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",
    "R = Y - X @ Theta\n",
    "E = R / b\n",
    "s = np.sqrt((R**2).sum(axis=0) / (N - 2))\n",
    "S = E.T @ E / (N - 2)\n",
    "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",
    "sigma = s[2]\n",
    "cov_ls = sigma**2 * np.linalg.inv(Xq.T @ Xq)\n",
    "others = np.arange(C) != 2\n",
    "mu0 = Theta[:, others].mean(axis=1)\n",
    "sd0 = Theta[:, others].std(axis=1, ddof=1)\n",
    "Lambda0 = np.diag(1 / sd0**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",
    "S_post = np.linalg.inv(P)\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",
    "rng = np.random.default_rng(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-02",
   "metadata": {},
   "source": [
    "## Άσκηση 1. Γιατί μία επίλυση πίνακα ισοδυναμεί με δέκα ανεξάρτητες προσαρμογές;\n",
    "Κατασκευάστε το πρόβλημα με μπλοκ διαγώνιο πίνακα και συγκρίνετε τη λύση με την\n",
    "επίλυση πολλών δεξιών μελών. Η απόδειξη WLS βρίσκεται στο παράρτημα A.\n",
    "Εξετάστε και την αβεβαιότητα των συντελεστών: γιατί μετρά η απόσταση των θερμοκρασιών;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_big = np.kron(np.eye(C), X)\n",
    "theta_big = np.linalg.lstsq(X_big, Y.T.reshape(-1), rcond=None)[0]\n",
    "print(\"same fit:\", np.allclose(theta_big.reshape(C, 2).T, Theta))\n",
    "XtX_inv = np.linalg.inv(X.T @ X)\n",
    "sd_theta = np.sqrt(np.outer(np.diag(XtX_inv), s**2))\n",
    "print(\"sd [a; b]:\", sd_theta)\n",
    "print(\"S_TT:\", ((T_cal - T_cal.mean())**2).sum())\n",
    "\n",
    "z = (R / s).ravel()\n",
    "grid = np.linspace(-4, 4, 300)\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "ax.hist(z, bins=50, density=True, alpha=0.5)\n",
    "ax.plot(grid, stats.norm.pdf(grid))\n",
    "ax.set(xlabel=\"r / s\", ylabel=tr(\"density\", \"πυκνότητα\"))\n",
    "plt.tight_layout(); plt.show()\n",
    "print(\"P(|r/s| > 2):\", np.mean(np.abs(z) > 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-04",
   "metadata": {},
   "source": [
    "## Άσκηση 2. Πώς φαίνεται πειστική η τοπική προσαρμογή ενώ το κέρδος είναι αναξιόπιστο;\n",
    "Σχεδιάστε πιθανές ευθείες, προσομοιώστε επαναλαμβανόμενους ελέγχους και εξετάστε\n",
    "την ασθενή διεύθυνση της SVD. Σχεδιάστε τους συντελεστές φίλτρου της ridge.\n",
    "Γιατί ένα μεγάλο εύρος κανονικοποίησης δίνει σχεδόν ίδιους συντελεστές;\n",
    "Δείτε το παράρτημα B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "one_bath = T_q == 25.0\n",
    "theta_one = np.linalg.lstsq(Xq[one_bath], y_q[one_bath], rcond=None)[0]\n",
    "print(\"one-bath rank / minimum-norm fit:\", np.linalg.matrix_rank(Xq[one_bath]), theta_one)\n",
    "fan = rng.multivariate_normal(theta_ls, cov_ls, size=40)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.4))\n",
    "for ax, limits in zip(axes, [(24.9, 25.3), (0, 45)]):\n",
    "    t = np.array(limits)\n",
    "    for theta in fan:\n",
    "        ax.plot(t, theta[0] + theta[1] * t, color=\"red\", alpha=0.2, lw=0.6)\n",
    "    ax.plot(t, theta_full[0] + theta_full[1] * t, \"k--\")\n",
    "    ax.plot(T_q, y_q, \"k.\")\n",
    "    ax.set(xlabel=\"T [°C]\", ylabel=\"counts\", xlim=limits)\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "Y_sim = (Xq @ Theta[:, 2])[:, None] + sigma * rng.standard_normal((len(T_q), 5000))\n",
    "b_sim = np.linalg.lstsq(Xq, Y_sim, rcond=None)[0][1]\n",
    "print(\"gain: mean / sd / P(b<0):\", b_sim.mean(), b_sim.std(), np.mean(b_sim < 0))\n",
    "_, sv, Vt = np.linalg.svd(Xq, full_matrices=False)\n",
    "print(\"singular values / weak direction:\", sv, Vt[1])\n",
    "\n",
    "def ridge(lam):\n",
    "    return np.linalg.solve(Xq.T @ Xq + lam * np.eye(2), Xq.T @ y_q)\n",
    "\n",
    "lams = np.logspace(-7, 5, 120)\n",
    "path = np.array([ridge(lam) for lam in lams])\n",
    "fig, axes = plt.subplots(1, 3, figsize=(11, 3.2))\n",
    "for i in range(2):\n",
    "    axes[0].semilogx(lams, sv[i]**2 / (sv[i]**2 + lams), label=f\"f{i+1}\")\n",
    "    axes[i+1].semilogx(lams, path[:, i])\n",
    "    axes[i+1].axhline(theta_full[i], color=\"black\", ls=\"--\")\n",
    "axes[0].legend(); axes[1].set_ylabel(\"a [counts]\"); axes[2].set_ylabel(\"b [counts/°C]\")\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"λ\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-06",
   "metadata": {},
   "source": [
    "## Άσκηση 3. Τι προσφέρουν η εκ των προτέρων κατανομή και η πιθανοφάνεια;\n",
    "Το διάγραμμα υπολογίζει απευθείας τη λογαριθμική πιθανοφάνεια και τις εκ των\n",
    "προτέρων και εκ των υστέρων κατανομές. Οι καμπύλες αβεβαιότητας και οι\n",
    "επαναλαμβανόμενοι έλεγχοι διακρίνουν την εκ των υστέρων αβεβαιότητα από τη\n",
    "μεταβλητότητα δειγματοληψίας. Αλλάξτε τη μέση τιμή και τις κλίμακες της εκ των\n",
    "προτέρων κατανομής. Μικρή αβεβαιότητα δεν συνεπάγεται μικρή μεροληψία.\n",
    "Αποδείξεις στο παράρτημα B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-07",
   "metadata": {},
   "outputs": [],
   "source": [
    "A, B = np.meshgrid(np.linspace(-700, 500, 500), np.linspace(32, 72, 500))\n",
    "residual = y_q[:, None, None] - A[None] - B[None] * T_q[:, None, None]\n",
    "loglik = -0.5 * (residual**2).sum(axis=0) / sigma**2\n",
    "logprior = -0.5 * (((A-mu0[0])/sd0[0])**2 + ((B-mu0[1])/sd0[1])**2)\n",
    "logpost = loglik + logprior\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "for values, color, label in [(loglik, \"red\", \"likelihood\"),\n",
    "                              (logprior, \"blue\", \"prior\"), (logpost, \"purple\", \"posterior\")]:\n",
    "    ax.contour(A, B, values - values.max(), levels=[-4.5, -2, -0.5], colors=color)\n",
    "    ax.plot([], [], color=color, label=label)\n",
    "ax.plot(*theta_map, \"ko\")\n",
    "ax.set(xlabel=\"a [counts]\", ylabel=\"b [counts/°C]\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "\n",
    "T_line = np.linspace(0, 50, 201)\n",
    "XT = np.column_stack([np.ones_like(T_line), T_line])\n",
    "curves = {\"LS\": cov_ls, \"prior\": np.diag(sd0**2), \"MAP\": S_post,\n",
    "          \"full\": sigma**2 * np.linalg.inv(Xf.T @ Xf), \"rack\": sigma**2 * XtX_inv}\n",
    "fig, ax = plt.subplots(figsize=(8, 3.4))\n",
    "for label, covariance in curves.items():\n",
    "    sd_line = np.sqrt(np.einsum(\"ti,ij,tj->t\", XT, covariance, XT)) / theta_full[1]\n",
    "    ax.semilogy(T_line, sd_line, label=label)\n",
    "ax.set(xlabel=\"T [°C]\", ylabel=\"sd [°C]\"); ax.legend()\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "fits = {\"LS\": np.linalg.lstsq(Xq, Y_sim, rcond=None)[0],\n",
    "        \"ridge\": np.linalg.solve(Xq.T @ Xq + np.eye(2), Xq.T @ Y_sim),\n",
    "        \"MAP\": np.linalg.solve(P, Xq.T @ Y_sim / sigma**2 + (Lambda0 @ mu0)[:, None])}\n",
    "for label, fitted in fits.items():\n",
    "    err = (np.array([1,40]) @ (fitted - Theta[:, 2, None])) / b[2]\n",
    "    print(label, \"bias / sd / RMSE:\", err.mean(), err.std(), np.sqrt((err**2).mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-08",
   "metadata": {},
   "source": [
    "## Άσκηση 4. Γιατί η διδιάστατη «έλλειψη 2 σίγμα» περιέχει μόνο 86%;\n",
    "Προσομοιώστε χ² σε διαφορετικές διαστάσεις και σχεδιάστε ελλείψεις στις αρχικές\n",
    "συντεταγμένες. Συγκρίνετε λευκοποίηση Cholesky και συμμετρική λευκοποίηση:\n",
    "διαφορετικές συντεταγμένες, ίδιες τετραγωνικές αποστάσεις. Παραρτήματα D και E."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-09",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(11, 3.1))\n",
    "for ax, df in zip(axes, [1, 2, 10]):\n",
    "    draws = (rng.standard_normal((20000, df))**2).sum(axis=1)\n",
    "    ax.hist(draws, bins=70, density=True, alpha=0.5)\n",
    "    grid = np.linspace(0.01, stats.chi2.ppf(0.999, df), 300)\n",
    "    ax.plot(grid, stats.chi2.pdf(grid, df))\n",
    "    ax.axvline(stats.chi2.ppf(0.99, df), color=\"red\", ls=\"--\")\n",
    "    ax.set_title(f\"χ²({df})\")\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "S2 = S[:2, :2]\n",
    "values, vectors = np.linalg.eigh(S2)\n",
    "t = np.linspace(0, 2*np.pi, 200)\n",
    "fig, ax = plt.subplots(figsize=(4.5, 4))\n",
    "ax.plot(E[:, 0], E[:, 1], \".\", alpha=0.4)\n",
    "for c in [1, 2, 3]:\n",
    "    ellipse = vectors @ (c*np.sqrt(values)[:, None] * np.array([np.cos(t), np.sin(t)]))\n",
    "    ax.plot(*ellipse, label=f\"c={c}, P={stats.chi2.cdf(c*c, 2):.0%}\")\n",
    "ax.set(xlabel=\"e1 [°C]\", ylabel=\"e2 [°C]\", aspect=\"equal\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "L = np.linalg.cholesky(S)\n",
    "values, vectors = np.linalg.eigh(S)\n",
    "W = (vectors * values**-0.5) @ vectors.T\n",
    "q_chol = np.linalg.solve(L, E.T).T\n",
    "q_sym = E @ W\n",
    "print(\"same distances:\", np.allclose((q_chol**2).sum(1), (q_sym**2).sum(1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-10",
   "metadata": {},
   "source": [
    "## Άσκηση 5. Γιατί η συσχέτιση διευκολύνει την ανίχνευση μικρής βλάβης ενός αισθητήρα;\n",
    "Εδώ **χρησιμοποιούμε τη γνωστή αλήθεια** των δεδομένων επικύρωσης: ο αισθητήρας 7\n",
    "αποκτά μετατόπιση +0.30 °C στη λήψη 200. Με τη θερμοκρασία αναφοράς συγκρίνουμε\n",
    "τα z-scores με την απόσταση χ²(C). Το κύριο hands-on διαγιγνώσκει άγνωστες\n",
    "θερμοκρασίες με C−1 βαθμούς ελευθερίας. Το παράρτημα E αναπτύσσει τη δεσμευμένη\n",
    "διακύμανση και την ισχύ μέσω της μη κεντρικής χ²."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "h1-ext-11",
   "metadata": {},
   "outputs": [],
   "source": [
    "ops = np.loadtxt(\"data/operations.csv\", delimiter=\",\", skiprows=1)\n",
    "errors = (ops[:, 2:] - a) / b - ops[:, 1, None]\n",
    "q = np.linalg.solve(L, errors.T).T\n",
    "d2 = (q**2).sum(axis=1)\n",
    "threshold = stats.chi2.ppf(0.99, C)\n",
    "precision = np.linalg.inv(S)\n",
    "j = 6\n",
    "sd_total = np.sqrt(S[j, j])\n",
    "sd_conditional = 1 / np.sqrt(precision[j, j])\n",
    "print(\"total / conditional sd:\", sd_total, sd_conditional)\n",
    "delta = np.linspace(0, 0.8, 161)\n",
    "zcrit = stats.norm.ppf(0.995)\n",
    "power_z = stats.norm.sf(zcrit - delta/sd_total) + stats.norm.cdf(-zcrit - delta/sd_total)\n",
    "power_m = stats.ncx2.sf(threshold, C, delta**2 * precision[j, j])\n",
    "fig, ax = plt.subplots(figsize=(8, 3.4))\n",
    "ax.plot(delta, power_z, label=\"z\")\n",
    "ax.plot(delta, power_m, label=\"Mahalanobis\")\n",
    "ax.set(xlabel=\"δ [°C]\", ylabel=tr(\"alarm probability\", \"πιθανότητα συναγερμού\"))\n",
    "ax.legend(); plt.tight_layout(); plt.show()\n",
    "print(\"observed before / after:\", (d2[:200] > threshold).mean(), (d2[200:] > threshold).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "h1-ext-12",
   "metadata": {},
   "source": [
    "## Πρόσθετες ασκήσεις\n",
    "- Εκτιμήστε τη συνδιακύμανση με μόνο 25 λήψεις από κάθε θερμοκρασία βαθμονόμησης.\n",
    "  Πώς αλλάζει ο ρυθμός ψευδών συναγερμών σε νέα δεδομένα; Το παράρτημα C εξηγεί\n",
    "  τη μεροληψία της αντίστροφης εκτιμημένης συνδιακύμανσης.\n",
    "- Δημιουργήστε καθαρό πείραμα λειτουργίας πριν αντικαταστήσετε τη βηματική βλάβη\n",
    "  με βαθμιαία ολίσθηση ή πρόσθετο θόρυβο. Μην προσθέσετε τη νέα βλάβη πάνω στην υπάρχουσα.\n",
    "- Δοκιμάστε δύο ελαττωματικούς αισθητήρες. Γιατί ο κανόνας του κύριου hands-on\n",
    "  μπορεί να διατηρήσει όλους τους αισθητήρες και να δηλώσει αβέβαιη διάγνωση;\n",
    "- Αξιολογήστε κάθε τροποποιημένο διαγνωστικό κανόνα σε χωριστό νέο προσομοιωμένο πείραμα."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}