Grafici con matplotlib 3

matplotlib e numpy sempre presenti

Guida completa: matplotlib.pyplot.plot

Full docs (numpy è parte della suite scipy) numpy

Per chi viene da matlab: numpy for matlab users

FIT dei dati

Uso di polyfit o curve_fit

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

Creo un po' di dati e provo a fittare.

In [2]:
a = np.linspace(0,10,10)
b = a + np.random.normal(0,2,10) 
plt.plot(a,b,'o')
plt.show()

Fit lineare. Semplicissimo con polyfit di numpy

In [3]:
p = np.polyfit(a,b,1)
print('p è un ',type(p),', ha dimensione ',len(p),' e il contenuto è ',p)

p è un  <class 'numpy.ndarray'> , ha dimensione  2  e il contenuto è  [0.95297913 1.09643709]
In [4]:
plt.plot(a,b,'o')
plt.plot(a,p[1]+p[0]*a,'r--')
plt.show()

Accedo ai risultati del fit usando l'opzione full=True quando fitto

In [5]:
p = np.polyfit(a,b,1, cov=True)
print('p è un ',type(p),', ha dimensione ',len(p),' e il contenuto è ',p)

p è un  <class 'tuple'> , ha dimensione  2  e il contenuto è  (array([0.95297913, 1.09643709]), array([[ 0.01543672, -0.07718359],
       [-0.07718359,  0.5431438 ]]))

Accedo ai coefficienti e alla matrice di covarianza separatamente (intercetto il return value)

In [6]:
p, cov = np.polyfit(a,b,1, cov=True)
print(p)
print('---')
print(cov)

[0.95297913 1.09643709]
---
[[ 0.01543672 -0.07718359]
 [-0.07718359  0.5431438 ]]

L'errore sui parametri è dato da sqrt degli elementi sulla diagonale della matrice di covarianza

In [7]:
errp = np.sqrt(np.diag(cov))
for n,pp in enumerate(p):
    print('parametro ',n,pp,'+-',errp[n])

parametro  0 0.9529791291199698 +- 0.12424459172776187
parametro  1 1.0964370943363224 +- 0.7369829045938866

Scrivo sul plot

In [8]:
plt.plot(a,b,'o', label='data')
plt.plot(a,p[1]+p[0]*a,'r--', label='y=%.2fx+%.2f'%(p[0],p[1]))
plt.text(0,6,'Fit errors: %.2f, %.2f'%(errp[0], errp[1]))
plt.legend()
plt.show()

Gaus

In [9]:
a = np.random.normal(10,3,10000) 
plt.hist(a,100, histtype='step', ec='b')
plt.show()
In [10]:
def gaus(x, A, mu, sigma):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

import scipy.optimize as opt

a = np.random.normal(10,3,10000) 
h, bins, _ = plt.hist(a,100, histtype='step', ec='b', label='data') ## cacthcing the return values of hist -> bin and values
p, cov = opt.curve_fit(gaus,bins[:-1],h, p0=[10,10,2] )
plt.plot(bins,gaus(bins,*p),'--r', label='fit')
plt.legend()
plt.show()

Fittare subrage (un po' tricky)

In [11]:
a = np.random.normal(10,3,10000) 
h, bins, _ = plt.hist(a,100, histtype='step', ec='b', label='data')

# tuple defining the range
subrange=(10,15)

# slice h and bins array to fit the range (use comprehension)
subh = [e[0] for e in zip(h,bins) if subrange[0] < e[1] < subrange[1]]
subb = [e[1] for e in zip(h,bins) if subrange[0] < e[1] < subrange[1]]
p, cov = opt.curve_fit(gaus,subb,subh, p0=[10,10,2] )
errp = np.sqrt(np.diag(cov))
plt.plot(subb,gaus(subb,*p),'--r', label='fit')
plt.legend()
plt.show()
print('Parametri fit:')
for n,el in enumerate(p):
    print('Par %d\t %6.1f\t\u00B1 %.1f'%(n,p[n],errp[n])) ## \00B1 is the unicode string for +- symb

Parametri fit:
Par 0	  274.2	± 8.8
Par 1	   10.3	± 0.3
Par 2	    2.8	± 0.2

Errori nei dati e chisq test

Test manuale del chi quadro $$\chi = \sum\frac{\left(observed-expected\right)^2}{expected}$$ L'array degli observed (i dati) è in questo caso subb (il sottoinsieme dei bin nel range ricavato prima) mentre gli expected sono i valori della funzione di fit

In [12]:
chi = ((subh-gaus(subb,*p))**2 / gaus(subb,*p) ).sum()
print('Chi = ',chi)

Chi =  23.414285866493522

Si può usare però il modulo stats per fare lo stesso (dà anche il pvalue)

In [13]:
import scipy.stats as stats

kk = stats.chisquare(subh,gaus(subb,*p))

print('Chi = ',kk[0])
print('pvalue = ', kk[1])

Chi =  23.414285866493522
pvalue =  0.3786613947434653