From 4a3432180044e476e09f3f1194d04465b9973205 Mon Sep 17 00:00:00 2001 From: Jacopo Basso Ricci Date: Sun, 27 Sep 2026 21:51:48 +0200 Subject: [PATCH] Aggiunta la libreria base per la statistica --- .gitignore | 4 + .python-version | 1 + README.md | 2 + pyproject.toml | 23 ++ statistica.py | 638 +++++++++++++++++++++++++++++++++++++++++++++ test_statistica.py | 295 +++++++++++++++++++++ 6 files changed, 963 insertions(+) create mode 100644 .gitignore create mode 100644 .python-version create mode 100644 README.md create mode 100644 pyproject.toml create mode 100644 statistica.py create mode 100644 test_statistica.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7265069 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.venv/ +uv.lock +__marimo__/ +__pycache__/ diff --git a/.python-version b/.python-version new file mode 100644 index 0000000..6324d40 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.14 diff --git a/README.md b/README.md new file mode 100644 index 0000000..1c17231 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# LOFM +Tutto il codice che useremo per LOFM diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..4989bfa --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,23 @@ +[project] +name = "lofm" +version = "0.1.0" +description = "Add your description here" +readme = "README.md" +requires-python = ">=3.14" +dependencies = [ + "marimo>=0.25.0", + "matplotlib>=3.11.2", + "numpy>=2.5.3", + "pandas>=3.0.6", + "scipy>=1.18.1", + "seaborn>=0.13.2", + "sympy>=1.14.0", + "tabulate>=0.10.0", +] + +[project.scripts] +lofm = "lofm:main" + +[build-system] +requires = ["uv_build>=0.12.19,<0.13.0"] +build-backend = "uv_build" diff --git a/statistica.py b/statistica.py new file mode 100644 index 0000000..4ad75b5 --- /dev/null +++ b/statistica.py @@ -0,0 +1,638 @@ +import marimo + +__generated_with = "0.25.0" +app = marimo.App() + + +@app.cell +def _(): + import marimo as mo + + return (mo,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Import librerie + """) + return + + +@app.cell +def _(): + # Heavy lifting + import numpy as np + import pandas as pd + from scipy import stats + + # Mostrare i dati + import matplotlib.pyplot as plt + import seaborn as sns + + # Calcolo simbolico + import sympy as sp + + return np, pd, plt, sns, sp, stats + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Compatibilità + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + **compat( x1, x2, u1, u2 )** + + *input:* + - x1, x2: quantità + - u1, u2: incertezze associate + + *output:* + - k: fattore di compatibilità (scarto normalizzato) + - diff: differenza x1 - x2 + """) + return + + +@app.cell +def _(np): + def compat(x1, x2, u1, u2): + diff = x1 - x2 + k = np.abs(diff) / np.sqrt(u1**2 + u2**2) + return k, diff + + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Classe Data + + *attributi:* + + - campione: dati grezzi + - analisi_stat: + - valori medi + - incertezze strumentali, statistiche e complessive + - lunghezza campione + - outlier + - residui con incertezza + - param_reg: + - A, B + - uA, uB, covAB + - chi², P + """) + return + + +@app.class_definition +class Data: + def __init__(self, + campione, + analisi_stat=None, param_reg=None): + + self.campione = campione + self.analisi_stat = analisi_stat + self.param_reg = param_reg + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Analisi con criterio di Chauvenet + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + **data.stat_chauv( prefissi, err_strumentale )** + + *input:* + - prefix: nome/lista di nomi delle variabili da analizzare + - err_strumentale: incertezza/lista di incertezze da risoluzione strumentale + + *output data.analisi_stat:* + - \: media del campione + - u\_strum: incertezza da risoluzione strumentale + - u\_stat: incertezza statistica (s/sqrt(N)) + - u\: incertezza complessiva + - n\: lunghezza del campione + - out\: lista outlier + """) + return + + +@app.cell +def _(np, pd, stats): + def p_t_student(valori, err_strumentale): + n = len(valori) + GdL = n - 1 + media = valori.mean() + s = np.sqrt(valori.std(ddof=1)**2 + err_strumentale**2) + return 2 * (1 - stats.t.cdf(np.abs(valori - media) / s, df=GdL)) + + + def trova_outlier(valori, err_strumentale): + n = len(valori) + soglia = 1 / (2 * n) + p = p_t_student(valori, err_strumentale) + idx_min = np.argmin(p) + if p[idx_min] < soglia: + return idx_min + return None + + + def rimuovi_outlier(valori, err_strumentale): + rimossi = [] + campione = valori.copy() + while len(campione) > 2: + idx = trova_outlier(campione, err_strumentale) + if idx is None: + break + rimossi.append(campione[idx]) + campione = np.delete(campione, idx) + return campione, rimossi + + + def stat_chauv(self, prefissi, err_strumentale): + campione_df = self.campione # riferimento al DataFrame originale (non sovrascritto) + + if isinstance(prefissi, str): + prefissi = [prefissi] + if isinstance(err_strumentale, (int, float)): + err_strumentale = [err_strumentale] * len(prefissi) + + # Pre-crea il DataFrame con indice corretto e colonne object per le liste + self.analisi_stat = pd.DataFrame(index=campione_df.index) + for prefix in prefissi: + self.analisi_stat[f"out{prefix}"] = pd.Series(dtype=object, + index=campione_df.index) + + for prefix, err_strum in zip(prefissi, err_strumentale): + cols = [col for col in campione_df.columns if col.startswith(prefix)] + + for i, row in campione_df.iterrows(): + valori = row[cols].dropna().values.astype(float) + valori_puliti, rimossi = rimuovi_outlier(valori, err_strum) + + n = len(valori_puliti) + s = valori_puliti.std(ddof=1) / np.sqrt(n) + u = np.sqrt(s**2 + err_strum**2) + + self.analisi_stat.at[i, prefix] = valori_puliti.mean() + self.analisi_stat.at[i, f"u{prefix}_strum"] = err_strum + self.analisi_stat.at[i, f"u{prefix}_stat"] = s + self.analisi_stat.at[i, f"u{prefix}"] = u + self.analisi_stat.at[i, f"n{prefix}"] = n + self.analisi_stat.at[i, f"out{prefix}"] = rimossi + + return self + + Data.stat_chauv = stat_chauv + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Calcolo variabile + Calcola una nuova variabile a partire dai valori medi e propagando le incertezze in data.analisi_stat + + **calc_var( funz, nome_risultato )** + + *input:* + - funz: funzione simbolica della variabile da calcolare + - nome_risultato: stringa con nome della variabile da calcolare + + *output data.analisi_stat:* + - \: valore della variabile in ogni punto sperimentale + - u\: incertezza + - %u\_\: contributo percentuale di ogni incertezza propagata + """) + return + + +@app.cell +def _(np, sp): + def calc_var(self, funz, nome_risultato): + # Estrazione e ordinamento variabili simboliche + vars_list = sorted(list(funz.free_symbols), key=lambda s: s.name) + + # Creazione incertezze e contributi simbolici + u = {var: sp.Symbol(f'u{var.name}', positive=True) for var in vars_list} + contr = [sp.diff(funz, var)*u[var] for var in vars_list] + + # Calcolo simbolico incertezza propagata + u_prop = 0 + for contributo in contr: + u_prop += contributo**2 + u_prop = sp.sqrt(u_prop) + + # Creazione lista incertezze simboliche + u_list = [u[var] for var in vars_list] + + # Funzioni numeriche + funz_fn = sp.lambdify(vars_list, funz, 'numpy') + u_prop_fn = sp.lambdify((vars_list, u_list), u_prop, 'numpy') + contr_fn = sp.lambdify((vars_list, u_list), contr, 'numpy') + + # Estrazione dati da analisi_stat + vars_num = [self.analisi_stat[var.name].values.astype(float) for var in vars_list] + u_num = [self.analisi_stat[f'u{var.name}'].values.astype(float) for var in vars_list] + + # Calcolo numerico risultati + risultato_val = funz_fn(*vars_num) + risultato_u = u_prop_fn(vars_num, u_num) + risultato_contr = (np.array(contr_fn(vars_num, u_num)) / risultato_u)**2 * 100 + + # Aggiunta risultati ad analisi_stat + self.analisi_stat[nome_risultato] = risultato_val + self.analisi_stat[f"u{nome_risultato}"] = risultato_u + for i, var in enumerate(vars_list): + self.analisi_stat[f"%u{nome_risultato}_{var.name}"] = risultato_contr[i] + + return self + + # Aggiunta a classe Data + Data.calc_var = calc_var + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Regressione lineare + y = Ax + By + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + **data.reg_lin( stampa_param=True, plot_regressione=True, calc_residui=True,** +        **x_label="\", y_label="\", r_label="\",** +        **titolo_reg="\", titolo_residui="\",** +        **cd_A=4, cd_B=4, scala_barre=1 )** + + *input data.analisi_stat:* + - x, y: campioni x e y + - ux, uy: relative incertezze + + *input:* + - stampa_param: opzione stampa parametri regressione + - plot_regressione: opzione plot regressione lineare + - residui: opzione calcolo e plot residui + - x_label, y_label, r_label: etichette assi plot regressione e residui + - titolo_reg, titolo_residui: titolo plot regressione e residui + - cd_A, cd_B: cifre decimali visualizzazione parametri + - scala barre: scala ingrandimento barre di errore nella regressione + + *output param_reg:* + - A, B: parametri regressione + - uA, uB, covAB: relative incertezze e covarianze + - chi², P: chi quadro e relativa probabilità + + *output analisi_stat:* + - r, ur: residui e relativa incertezza + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ### Residui + """) + return + + +@app.cell +def _(plt, sns): + # Plot residui + + def plot_residui(x, r, + ux, ur, + x_label, r_label, titolo): + fig, ax = plt.subplots(figsize=(8, 5)) + # Residui con barre d'errore + ax.errorbar( + x, r, + xerr=ux, yerr=ur, + fmt='o', color=sns.color_palette()[0], + ecolor='gray', elinewidth=1, capsize=3, + markersize=5, label="Residui" + ) + + # Linea dello zero + ax.axhline(0, color='red', linestyle='--', linewidth=1) + + # Estetica + sns.despine(ax=ax) + ax.set_xlabel(x_label) + ax.set_ylabel(r_label) + ax.set_title(titolo) + ax.legend() + ax.grid(True, linestyle=':', linewidth=0.5, alpha=0.7) + + plt.tight_layout() + plt.show() + + return (plot_residui,) + + +@app.cell +def _(sp): + # Calcolo residui + + def residui(analisi_stat, + A_num, B_num, + uA_num, uB_num, covAB_num): + + # Variabili simboliche + A, B, x, y = sp.symbols('A B x y', real=True) + uA, uB, ux, uy = sp.symbols('uA uB ux uy', positive=True) + covAB = sp.symbols('covAB', real=True) + + # Residuo: r = y - (Ax + B) + r = y - (A*x + B) + + # Propagazione errore (con covarianza A,B) + dr_dA = sp.diff(r, A) + dr_dB = sp.diff(r, B) + dr_dx = sp.diff(r, x) + dr_dy = sp.diff(r, y) + + u_r = sp.sqrt( + (dr_dA * uA)**2 + + (dr_dB * uB)**2 + + (dr_dx * ux)**2 + + (dr_dy * uy)**2 + + 2 * dr_dA * dr_dB * covAB + ) + + # Funzioni numeriche + r_fn = sp.lambdify((x , y , A , B ), r, 'numpy') + u_r_fn = sp.lambdify( + (x , y , ux , uy , A , B , uA , uB , covAB ), + u_r , 'numpy' + ) + + # Calcolo numerico + analisi_stat["r"] = r_fn( + analisi_stat["x"], + analisi_stat["y"], + A_num, + B_num + ) + + analisi_stat["ur"] = u_r_fn( + analisi_stat["x"], analisi_stat["y"], + analisi_stat["ux"], analisi_stat["uy"], + A_num, B_num, + uA_num, uB_num, covAB_num + ) + + return analisi_stat + + return (residui,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ### Regressione + """) + return + + +@app.cell +def _(np, plt, sns): + # Plot regressione + + def plot_reg(x, y, ux, uy, + A, B, uA, uB, P, + x_label, y_label, titolo, + cd_A=4, cd_B=4, scala_barre=1): + + fig, ax = plt.subplots(figsize=(8, 5)) + + x_fit = np.linspace(x.min(), x.max(), 300) + y_fit = A * x_fit + B + + ax.errorbar( + x, y, + xerr = scala_barre * ux, + yerr = scala_barre * uy, + fmt='o', color=sns.color_palette()[0], + ecolor='gray', elinewidth=1, capsize=3, + markersize=5, label=f"Dati (barre errore x{scala_barre})" + ) + ax.plot( + x_fit, y_fit, + color='red', linewidth=1.5, + label = f"$A={A:.{cd_A}f}\\pm{uA:.{cd_A}f}$\n" + f"$B={B:.{cd_B}f}\\pm{uB:.{cd_B}f}$\n" + f"$P(\\chi^2, \\infty)={P:.4f}$\n" + ) + + sns.despine(ax=ax) + ax.set_xlabel(x_label) + ax.set_ylabel(y_label) + ax.set_title(titolo) + ax.legend(fontsize=9) + ax.grid(True, linestyle=':', linewidth=0.5, alpha=0.7) + + plt.tight_layout() + plt.show() + + return (plot_reg,) + + +@app.cell +def _(np, pd, plot_reg, plot_residui, residui, stats): + # Funzioni regressione + + def uy_equiv(x, y, ux, uy): + # Stima iniziale di A con sola uy + sy2 = uy**2 + Sw = np.sum(1 / sy2) + Sx = np.sum(x / sy2) + Sxx = np.sum(x**2 / sy2) + Sy = np.sum(y / sy2) + Sxy = np.sum(x * y / sy2) + delta = Sxx * Sw - Sx**2 + A_est = (Sxy * Sw - Sx * Sy) / delta + + # Propagazione + u_eq = np.sqrt(uy**2 + A_est**2 * ux**2) + return u_eq + + + def reg_lin(self, + stampa_param=True, plot_regressione=True, calc_residui=True, + x_label="", y_label="", r_label="", + titolo_reg="", titolo_residui="", + cd_A=4, cd_B=4, scala_barre=1): + + x = np.asarray(self.analisi_stat["x"], dtype=float) + y = np.asarray(self.analisi_stat["y"], dtype=float) + ux = np.asarray(self.analisi_stat["ux"], dtype=float) + uy = np.asarray(self.analisi_stat["uy"], dtype=float) + + # Propagazione errore x → y + uy_eq = uy_equiv(x, y, ux, uy) + + # Somme pesate + w = 1.0 / uy_eq**2 + Sw = np.sum(w) + Sx = np.sum(w * x) + Sxx = np.sum(w * x**2) + Sy = np.sum(w * y) + Sxy = np.sum(w * x * y) + delta = Sxx * Sw - Sx**2 + + # Parametri + A = (Sxy * Sw - Sx * Sy) / delta + B = (Sxx * Sy - Sxy * Sx) / delta + uA = np.sqrt(Sw / delta) + uB = np.sqrt(Sxx / delta) + covAB = -Sx / delta + + # Chi quadro + x2 = np.sum((y - A * x - B)**2 / uy_eq**2) + dof = len(x) - 2 + P = stats.chi2.sf(x2, dof) + + # Raccolta parametri + self.param_reg = pd.DataFrame() + # Vecchi modo di caricare i parametri + # self.param_reg["A"] = A + # self.param_reg["B"] = B + # self.param_reg["uA"] = uA + # self.param_reg["uB"] = uB + # self.param_reg["covAB"] = covAB + # self.param_reg["x2"] = x2 + # self.param_reg["P"] = P + self.param_reg = pd.Series({ + "A": A, + "B": B, + "uA": uA, + "uB": uB, + "covAB": covAB, + "x2": x2, + "P": P + }) + + # Stampa + if stampa_param == True: + print("Ax + B : ") + print(f"A = {A:.{cd_A}f} ± {uA:.{cd_A}f}") + print(f"B = {B:.{cd_B}f} ± {uB:.{cd_B}f}") + print(f"covAB = {covAB:.6f}") + print(f"chi² = {x2:.2f}") + print(f"P(chi², ∞) = {P:.2f}") + + # Plot + if plot_regressione == True: + plot_reg(x, y, ux, uy, + A, B, uA, uB, P, + x_label, y_label, titolo_reg, + cd_A, cd_B, scala_barre) + + # Residui + if calc_residui == True: + self.analisi_stat = residui(self.analisi_stat, A, B, uA, uB, covAB) + plot_residui(x, self.analisi_stat["r"], + ux, self.analisi_stat["ur"], + x_label, r_label, titolo_residui) + + return self + + + # Aggiunta a classe Data + Data.reg_lin = reg_lin + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Plot gaussiane + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + **plot_gauss( gaussiane )** + + *input gaussiane:* + [ (mi, sigma, colore, label),... ] + - mi: valore medio (misura) + - sigma: deviazione standard (incertezza) + - colore: indice colore nella palette + - label: etichetta per la didascalia + """) + return + + +@app.cell +def _(np, plt, sns, stats): + def plot_gauss(gaussiane): + # Creazione figura + plt.figure(figsize=(12, 7)) + + # Creazione asse x + xMin = float('inf') + xMax = float('-inf') + + for mu, sigma, _, _ in gaussiane: + minimoLocale = mu - 4 * sigma + massimoLocale = mu + 4 * sigma + + if minimoLocale < xMin: + xMin = minimoLocale + + if massimoLocale > xMax: + xMax = massimoLocale + + x = np.linspace(xMin, xMax, 1000) + + # Ciclo gaussiane + for mu, sigma, colore, etichetta in gaussiane: + y = stats.norm.pdf(x, mu, sigma) + plt.plot(x, y, color=sns.color_palette()[colore], linewidth=1, label=etichetta) + + puntiLinee = [mu - sigma, mu, mu + sigma] + for px in puntiLinee: + py = stats.norm.pdf(px, mu, sigma) + plt.vlines(x=px, + ymin=0, ymax=py, + colors=sns.color_palette()[colore], + linestyles='dashed', linewidth=1) + + # Dettagli estetici finali + plt.ylim(bottom=0) + plt.title('Confronto dati') + plt.xlabel('k') + plt.legend() + + # Mostriamo il grafico + plt.show() + + return + + +if __name__ == "__main__": + app.run() diff --git a/test_statistica.py b/test_statistica.py new file mode 100644 index 0000000..4cb02f8 --- /dev/null +++ b/test_statistica.py @@ -0,0 +1,295 @@ +import marimo + +__generated_with = "0.25.0" +app = marimo.App(width="medium") + + +@app.cell +def _(): + import marimo as mo + + return (mo,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Test del modulo `statistica.py` + + Questo notebook importa `statistica.py` (notebook marimo) e testa **ogni + funzione** definita al suo interno, verificando che i risultati siano + corretti e segnalando eventuali errori. + + Poiché in `statistica.py` le funzioni sono definite dentro alle celle + `@app.cell` (non con `@app.function`), per poterle riutilizzare qui + eseguiamo l'intero notebook con `app.run()`, che restituisce un + dizionario con tutte le definizioni (`defs`). Solo la classe `Data` + (definita con `@app.class_definition`) sarebbe importabile direttamente. + """) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Import e caricamento del modulo da testare + """) + return + + +@app.cell +def _(): + import matplotlib + + import numpy as np + import pandas as pd + import sympy as sp + + import statistica + + # Esegue tutte le celle del notebook importato e recupera le definizioni + _outputs, defs = statistica.app.run() + + Data = defs["Data"] + compat = defs["compat"] + stat_chauv = defs["stat_chauv"] + calc_var = defs["calc_var"] + reg_lin = defs["reg_lin"] + plot_gauss = defs["plot_gauss"] + trova_outlier = defs["trova_outlier"] + rimuovi_outlier = defs["rimuovi_outlier"] + return Data, compat, np, pd, plot_gauss, sp + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Helper di test + """) + return + + +@app.cell +def _(mo): + risultati_test = [] + + def esegui_test(nome, fn): + """Esegue fn(), registra esito e ritorna un mo.md con il risultato.""" + try: + dettaglio = fn() + risultati_test.append((nome, True, "")) + corpo = f"✅ **{nome}**: OK" + if dettaglio: + corpo += f"\n\n{dettaglio}" + return mo.md(corpo) + except Exception as e: + import traceback + tb = traceback.format_exc() + risultati_test.append((nome, False, str(e))) + return mo.md(f"❌ **{nome}**: ERRORE — `{e}`\n\n```\n{tb}\n```") + + return esegui_test, risultati_test + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Test 1 — `compat(x1, x2, u1, u2)` + """) + return + + +@app.cell +def _(compat, esegui_test, np): + def _test_compat(): + k, diff = compat(10.0, 10.5, 0.2, 0.3) + assert np.isclose(diff, -0.5), f"diff errato: {diff}" + atteso_k = 0.5 / np.sqrt(0.2**2 + 0.3**2) + assert np.isclose(k, atteso_k), f"k errato: {k} atteso {atteso_k}" + return f"k = {k:.4f}, diff = {diff:.4f} (attesi: k={atteso_k:.4f}, diff=-0.5)" + + esegui_test("compat", _test_compat) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Test 2 — `Data.stat_chauv(prefissi, err_strumentale)` + + Costruiamo un campione con 4 punti sperimentali, ciascuno con 8 misure + ripetute della grandezza `T`. Nel primo punto inseriamo un valore + anomalo (~10σ) che il criterio di Chauvenet dovrebbe individuare e + rimuovere. + """) + return + + +@app.cell +def _(Data, esegui_test, np, pd): + np.random.seed(0) + _righe = [np.random.normal(10, 0.05, size=8) for _ in range(4)] + _righe[0][0] = 10.5 # outlier moderato ma rilevabile + + campione_chauv = pd.DataFrame(_righe, columns=[f"T{i+1}" for i in range(8)]) + dato_chauv = Data(campione_chauv) + + def _test_stat_chauv(): + dato_chauv.stat_chauv("T", 0.02) + colonne_attese = {"T", "uT_strum", "uT_stat", "uT", "nT", "outT"} + mancanti = colonne_attese - set(dato_chauv.analisi_stat.columns) + assert not mancanti, f"colonne mancanti in analisi_stat: {mancanti}" + assert len(dato_chauv.analisi_stat.at[0, "outT"]) >= 1, ( + "l'outlier atteso non è stato rimosso" + ) + return dato_chauv.analisi_stat.to_markdown() + + esegui_test("Data.stat_chauv", _test_stat_chauv) + return campione_chauv, dato_chauv + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Test 3 — `Data.calc_var(funz, nome_risultato)` + + Aggiungiamo una seconda grandezza `P` allo stesso `Data` e calcoliamo + `Z = T * P` con propagazione simbolica dell'incertezza. + """) + return + + +@app.cell +def _(campione_chauv, dato_chauv, esegui_test, np, pd, sp): + def _test_calc_var(): + np.random.seed(1) + campione_p = pd.DataFrame( + np.random.normal(2, 0.02, size=(4, 8)), + columns=[f"P{i+1}" for i in range(8)], + ) + dato_chauv.campione = pd.concat([campione_chauv, campione_p], axis=1) + dato_chauv.stat_chauv(["T", "P"], [0.02, 0.01]) + + T, P = sp.symbols("T P", positive=True) + dato_chauv.calc_var(T * P, "Z") + + for col in ["Z", "uZ", "%uZ_T", "%uZ_P"]: + assert col in dato_chauv.analisi_stat.columns, f"colonna mancante: {col}" + + T_val = dato_chauv.analisi_stat["T"].astype(float) + P_val = dato_chauv.analisi_stat["P"].astype(float) + Z_val = dato_chauv.analisi_stat["Z"].astype(float) + assert np.allclose(Z_val, T_val * P_val), "Z != T*P: propagazione errata" + + return dato_chauv.analisi_stat[["T", "uT", "P", "uP", "Z", "uZ", "%uZ_T", "%uZ_P"]].to_markdown() + + esegui_test("Data.calc_var", _test_calc_var) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Test 4 — `Data.reg_lin(...)` (incl. `residui`, `plot_reg`, `plot_residui`) + + Generiamo dati sintetici y = 2x + 1 con rumore gaussiano e verifichiamo + che la regressione pesata recuperi correttamente A e B, e che vengano + prodotti anche i grafici e i residui. + """) + return + + +@app.cell +def _(Data, esegui_test, np, pd): + def _test_reg_lin(): + np.random.seed(2) + n = 12 + A_vero, B_vero = 2.0, 1.0 + x = np.linspace(0, 10, n) + ux = np.full(n, 0.1) + uy = np.full(n, 0.2) + y = A_vero * x + B_vero + np.random.normal(0, 0.2, n) + + dato_reg = Data(pd.DataFrame({"x": x, "y": y, "ux": ux, "uy": uy})) + dato_reg.analisi_stat = dato_reg.campione.copy() + + dato_reg.reg_lin( + stampa_param=True, + plot_regressione=True, + calc_residui=True, + x_label="x", + y_label="y", + r_label="residuo", + titolo_reg="Test regressione", + titolo_residui="Test residui", + ) + + assert dato_reg.param_reg is not None, "param_reg non impostato" + for chiave in ["A", "B", "uA", "uB", "covAB", "x2", "P"]: + assert chiave in dato_reg.param_reg.index, f"manca {chiave} in param_reg" + + errore_A = abs(dato_reg.param_reg["A"] - A_vero) + assert errore_A < 0.3, f"A stimato troppo lontano dal vero: {dato_reg.param_reg['A']}" + assert "r" in dato_reg.analisi_stat.columns + assert "ur" in dato_reg.analisi_stat.columns + + p = dato_reg.param_reg + return ( + f"A = {p['A']:.4f} ± {p['uA']:.4f} (vero: {A_vero}) \n" + f"B = {p['B']:.4f} ± {p['uB']:.4f} (vero: {B_vero}) \n" + f"chi² = {p['x2']:.3f}, P = {p['P']:.3f}" + ) + + esegui_test("Data.reg_lin", _test_reg_lin) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Test 5 — `plot_gauss(gaussiane)` + """) + return + + +@app.cell +def _(esegui_test, plot_gauss): + def _test_plot_gauss(): + gaussiane = [ + (0, 1, 0, "N(0,1)"), + (2, 0.5, 1, "N(2, 0.5)"), + ] + plot_gauss(gaussiane) + return "Grafico generato senza errori per 2 gaussiane." + + esegui_test("plot_gauss", _test_plot_gauss) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Riepilogo finale + """) + return + + +@app.cell +def _(mo, risultati_test): + _n_ok = sum(1 for _, ok, _ in risultati_test if ok) + _n_tot = len(risultati_test) + _righe = [f"| {nome} | {'✅ OK' if ok else '❌ ERRORE: ' + err} |" for nome, ok, err in risultati_test] + _tabella = "| Funzione | Esito |\n|---|---|\n" + "\n".join(_righe) + + _titolo = ( + f"### ✅ Tutti i test sono passati ({_n_ok}/{_n_tot})" + if _n_ok == _n_tot + else f"### ⚠️ {_n_tot - _n_ok} test falliti su {_n_tot}" + ) + + mo.md(f"{_titolo}\n\n{_tabella}") + return + + +if __name__ == "__main__": + app.run()