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()