639 lines
17 KiB
Python
639 lines
17 KiB
Python
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:*
|
|
- \<prefix>: media del campione
|
|
- u\<prefix>_strum: incertezza da risoluzione strumentale
|
|
- u\<prefix>_stat: incertezza statistica (s/sqrt(N))
|
|
- u\<prefix>: incertezza complessiva
|
|
- n\<prexix>: lunghezza del campione
|
|
- out\<prefix>: 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:*
|
|
- \<nome_risultato>: valore della variabile in ogni punto sperimentale
|
|
- u\<nome_risultato>: incertezza
|
|
- %u\<nome_risultato>_\<var>: 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()
|