Estrarre gli errori di un fit

Gli errori dei fit si estraggono dalla matrice di covarianza che viene "restituita" dal comando di fit. Vedere in fondo come.

Per i più curiosi: discussione moooooolto approfondita

In [10]:
import uproot
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt

# ho trasformato histo di paw in root (con h2root)
histos = uproot.open('gamma_020006.hbk.root')
# prendo il 2000
h2000 = histos['h2000']

# per comodità metto gli edge dei bin in "x" e i valori in "y"
x = h2000.bins[:,0]
y = h2000.values

# plot 
plt.plot(x,y, drawstyle='steps', color='k')
plt.title('histo 2000')
plt.show()
In [11]:
# definisco esponenziale
def expo(x,A,B):
    return np.exp(A + B*x)

# range da fittare
sel = (x > 510) & (x < 1000)

# inizializzazione parametri (per fit non banali serve)
init_pars = [10, -0.01]
p, cov = opt.curve_fit(expo,x[sel],y[sel], p0=init_pars)

In p e cov si intercettano i "risultati" del fit. In particolare:

  • p -> parametri
  • cov -> è la matrice di covarianza.

Gli errori sui parametri si possono calcolare dalla radice quadrata della diagonale della matrice di covarianza, cioè:

In [12]:
errp = np.sqrt(np.diag(cov))

# stampo per leggerli
for i in range(2):
    print('Parametro',p[i],' -> errore',errp[i])

Parametro 12.219523133895741  -> errore 0.04049188230966112
Parametro -0.007034562443727151  -> errore 6.915133657987828e-05

Plot finale con risultato fit e relativi errori scritti sul plot

In [13]:
plt.plot(x,y, label='dati', drawstyle='steps', color='k')
label_da_mettere = f'fit exp(A+Bx)\nconst = {p[0]:.2f} ± {errp[0]:.2}\ntau = {p[1]:.2E} ± {errp[1]:.2E}'
plt.plot(x[sel], expo(x[sel], *p), '--r', linewidth=2, label=label_da_mettere)
plt.legend(fontsize=12)
plt.show()