creazione esperimento molle
This commit is contained in:
413
statlib.py
Normal file
413
statlib.py
Normal file
@@ -0,0 +1,413 @@
|
||||
|
||||
# Heavy lifting
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from scipy import stats
|
||||
|
||||
# Mostrare i dati
|
||||
import ipysheet
|
||||
from IPython.display import display
|
||||
import matplotlib.pyplot as plt
|
||||
import seaborn as sns
|
||||
|
||||
# Calcolo simbolico
|
||||
import sympy as sp
|
||||
|
||||
def compat(x1, x2, u1,u2):
|
||||
|
||||
diff = x1 - x2
|
||||
k = np.abs(diff) / np.sqrt(u1**2 + u2**2)
|
||||
return k, diff
|
||||
|
||||
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
|
||||
|
||||
# Probabilità
|
||||
def p_t_student(valori, err_strumentale):
|
||||
|
||||
n = len(valori)
|
||||
GdL = n - 1
|
||||
|
||||
media = valori.mean()
|
||||
s = valori.std(ddof=1)
|
||||
|
||||
s = np.sqrt(s**2 + err_strumentale**2)
|
||||
return 2 * (1 - stats.t.cdf(np.abs(valori - media) / s, df=GdL))
|
||||
|
||||
|
||||
# Indice del peggiore outlier o None
|
||||
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
|
||||
|
||||
|
||||
# Rimozione outlier
|
||||
def rimuovi_outlier(valori, err_strumentale):
|
||||
|
||||
rimossi = []
|
||||
campione = valori.copy()
|
||||
|
||||
while len(campione) > 2:
|
||||
idx = trova_outlier(campione, err_strumentale)
|
||||
|
||||
if idx is None: # nessun outlier: stop
|
||||
break
|
||||
|
||||
rimossi.append(campione[idx])
|
||||
campione = np.delete(campione, idx)
|
||||
|
||||
return campione, rimossi
|
||||
|
||||
# Analisi statistica con criterio di Chauvenet
|
||||
def stat_chauv(self, prefissi, err_strumentale):
|
||||
|
||||
self.analisi_stat = pd.DataFrame()
|
||||
campione = self.campione
|
||||
|
||||
|
||||
if isinstance(prefissi, str):
|
||||
prefissi = [prefissi]
|
||||
if isinstance(err_strumentale, (int, float)):
|
||||
err_strumentale = [err_strumentale] * len(prefissi)
|
||||
|
||||
for prefix, err_strum in zip(prefissi, err_strumentale):
|
||||
cols = [col for col in campione.columns if col.startswith(prefix)]
|
||||
|
||||
for i, row in campione.iterrows():
|
||||
valori = row[cols].dropna().values.astype(float)
|
||||
campione, rimossi = rimuovi_outlier(valori, err_strum)
|
||||
|
||||
n = len(campione)
|
||||
s = campione.std(ddof=1)
|
||||
s = s / np.sqrt(n)
|
||||
u = np.sqrt(s**2 + err_strum**2)
|
||||
|
||||
self.analisi_stat.at[i, prefix] = campione.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
|
||||
|
||||
|
||||
# Aggiunta a classe Data
|
||||
Data.stat_chauv = stat_chauv
|
||||
|
||||
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
|
||||
|
||||
# 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()
|
||||
|
||||
# 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
|
||||
|
||||
|
||||
# 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², ∞)={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()
|
||||
|
||||
|
||||
# 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()
|
||||
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
|
||||
|
||||
# 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
|
||||
|
||||
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()
|
||||
Reference in New Issue
Block a user