La funzione "spigolosa"

(Con l'aggiunta della correzione sui centri dei bin da usare per il fit)

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

# 1000 numeri gaussiani (media 10, sigma 1)
d = np.random.normal(10, 1, 1000)

# histo con un certo binning
h, bins, _ = plt.hist(d, bins=12, color='gray')
plt.show()
In [2]:
# NOTA: in h e bins ho raccolto ciò che viene restituito dal comando hist, vale a dire, rispettivamente, l'array dei valori
# e l'array dei bin. La variabile "_" viene usata per intercettare il resto da ignorare
# Per capire basta stamparli...
print(h)
print(bins)
[  4.  21.  55. 116. 182. 210. 185. 136.  65.  17.   7.   2.]
[ 6.89165402  7.44380641  7.99595881  8.5481112   9.10026359  9.65241598
 10.20456837 10.75672076 11.30887316 11.86102555 12.41317794 12.96533033
 13.51748272]
In [3]:
# Attenzione che il size di "bins" è quello di h +  1 (contiene gli edge dei bin)
# bisogna tenerne conto quando si usano insieme
print('h -> ',len(h))
print('bins -> ',len(bins))
h ->  12
bins ->  13
In [4]:
# definisco gaussiana per poi fare fit
def gaus(x, A, mu, sigma):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

# creo un array di "centri" dei bin in modo che il fit usi questi valori
# e non quelli dei bordi

binw = bins[1]-bins[0]
binc = [b+binw/2 for b in bins[:-1]] 


# faccio fit usando "curve_fit" 
p, _ = opt.curve_fit(gaus,binc,h, p0=[10,10,1])

# p contiene il risultato del fit
print('const ->', p[0])
print('mean ->', p[1])
print('sigma ->', p[2])
const -> 212.32241901681223
mean -> 9.979413859626975
sigma -> 1.045408050340565
In [5]:
# Se faccio il plot usando lo stesso asse dei dati questo avrà la stessa discretizzazione
# in particolare la funzione è valutata ogni edge di ogni bin
plt.hist(d, bins=12)
plt.plot(binc, gaus(binc, *p), '--r')
plt.show()
In [6]:
# Per un risultato graficamente migliore creo un asse x più fine. Ad esempio usando linspace() di numpy
# che crea un array di numeri equidistanziati fra due estremi. 

xfini = np.linspace(7, 13, 1000) # 1000 numeri equispaziati fra 7 e 13, estremi inclusi

plt.hist(d, bins=12)
plt.plot(xfini, gaus(xfini, *p), '--r')
plt.show()