Analisi statistiche in python
Ad esempio, possiamo fare una analisi statistica per vedere come si distribuisce una serie di dati, che ipotizziamo possano assumere qualunque valore reale (variabile continua). In genere, tale analisi inizia con lo scrivere la Funzione densità di Probabilità (PDF) \(P(x)\) dei dati, che rappresenta la probabilità \(p(x)\) che un dato valore \(x\) della variabile continua sia contenuto nell'intervallo \(x \div x+\Delta x\), divisa per la larghezza dell'intervallo \(\Delta x\), quando questo tende a zero:
\[P(x) = \lim_{\Delta x \to 0} \frac{p(x+\Delta x)-p(x)}{\Delta x} = \frac{dp(x)}{dx}\]
Segue dalla definizione che la probabilità di trovare un dato valore \(x\) nell'intervallo \([a,b]\) è data da:
\[p(x) = \int_{a}^{b} P(x) dx.\]
Naturalmente deve valere la proprietà:
\[\int_{-\infty}^{+\infty} P(x) dx = 1\]
in quanto è certo che \(x\) assumerà un valore reale.
Nella pratica, uno NON ha a disposizione un insieme infinito di valori reali, ma bensì un insieme discreto \(N\) di numeri reali \(x_j\) e, suddiviso l'intervallo \([x_{min},x_{max}]\) in un certo numero \(N_c\) di sottointervalli (canali, o bins) di ampiezza \(\Delta x\), tenuto conto della definizione di probabilità come rapporto tra il numero di casi favorevoli sul numero di casi possibili, si ha che:
\[P(x_i) = \frac{n_i}{\Delta x N}\]
dove \(n_i\) è il numero di valori di \(x\) che cadono nell'\(i\)-esimo sottointervallo.
Naturalmente se \(N \to \infty\) e \(\Delta x \to 0\), si ritrova la definizione data sopra, ma nella pratica questo non succede mai.
Il calcolo della PDF si riconduce a suddividere l'intervallo \([x_{min},x_{max}]\) in \(N_c\) sottointervalli e a contare quanti valori delle \(x_j\) cadono in ciascuno di questi sottointervalli (cioé valutare \(n_i\)), quindi dividere ciascun valore per \(\Delta x N\).
La libreria numpy di python contiene una funzione che calcola la pdf di un vettore x di valori:
pdf, bins = histogram( x, nbins, density = True|False )
La variabile booleana density fa sì che la funzione restituisca il numero dei conteggi \(n_i\) se density=False o la PDF se density=True.
Esempio: utilizziamo la funzione random.random_sample() per ottenere un vettore di valori distribuito secondo una distribuzione uniforme e la funzione histogram per plottare la pdf dei valori stessi:
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as ran
N = 100
x = ran.random_sample(N)
print(x)
print( "x minimo = %le, x massimo = %le " % (x.min(), x.max()) )
Nc = 10
pdf, bins = np.histogram( x, Nc, density = True )
print( pdf )
print( bins )
Possiamo quindi fare il grafico della PDF:
%matplotlib inline
nextr = len( bins )
dx = bins[1] - bins[0]
plt.xlabel( r"$x$" )
plt.ylabel( "PDF" )
plt.bar( bins[:nextr-1], pdf, dx, color = "g" )
plt.show()
Possiamo quindi verificare che l'integrale vale effettivamente 1 moltiplicando il valore della pdf per la larghezza degli intervalli e sommando:
\[\sum_{i=1}^{N_c} P(x_i) \Delta x\]
print( dx * pdf.sum() )
Naturalmente la distribuzione diventa più piatta se si sceglie N più grande, in modo da aumentare la statistica in ogni singolo canale:
N = 1000
x = ran.random_sample( N )
Nc = 10
pdf, bins = np.histogram( x, Nc, density = True )
nextr = len( bins )
dx = bins[1] - bins[0]
plt.xlabel( r"$x$" )
plt.ylabel( "PDF" )
plt.bar( bins[:nextr-1], pdf, dx, color = "g" )
plt.show()
Ad esempio, la funzione:
numpy.random.normal( xc, sigma, N )
genera in insieme di N numeri random distribuiti secondo una distribuzione Gaussiana di centro xc e deviazione standard sigma:
\[p(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]
dove \(\mu\) è il centro della distribuzione e \(\sigma\) la deviazione standard.
Esempio: generiamo N = 1000 valori distribuiti in maniera Gaussiana e plottiamone l'istogramma
x = ran.normal( 0.0, 1.0, 1000 )
Nc = 50
pdf, bins = np.histogram( x, Nc, density = True )
nextr = len( bins )
dx = bins[1] - bins[0]
plt.xlabel( r"$x$" )
plt.ylabel( "PDF" )
plt.bar( bins[:nextr-1], pdf, dx, color = "g" )
plt.show()
Il valore medio di una \(P(x)\) è definito come:
\[\bar{x} = \int_{-\infty}^{+\infty} x P(x) dx = \sum_{i=1}^{N_c} x_i P(x_i) \Delta x\]
e vale \(\mu\) per la distribuzione Gaussiana, \(0.5\) per la distribuzione uniforme, ecc.
La deviazione standard o momento del secondo ordine rappresenta come si sparpagliano i valori attorno alla media:
\[s = \int_{-\infty}^{+\infty} (x-\bar{x})^2 P(x) dx = \sum_{i=1}^{N_c} (x_i - \bar{x})^2 P(x_i) \Delta x\]
e vale esattamente \(\sigma\) per la Gaussiana, \(\frac{1}{12}\) per la distribuzione uniforme, ecc.
Più in generale, i momenti statistici di una distribuzione aiutano ad identificare il tipo di distribuzione seguita da una serie di numeri e sono definiti come:
\[m(n;P) = \int_{-\infty}^{+\infty} (x - \bar{x})^n P(x) dx = \sum_{i=1}^{N_c} (x_j - \bar{x})^n P(x_i) \Delta x\]
In genere il momento del terzo ordine si chiama Skewness e tiene conto delle simmetrie della funzione di distribuzione e vale 0 per una Gaussiana; quello del quarto ordine si chiama flatness (o kurtosis) e tiene conto di quanto sono riempite le code della distribuzione (vale \(3\sigma^2\) per una Gaussiana), ecc.
Esempio: definiamo una funzione che calcoli i momenti di ordine n di una PDF di dati, ad esempio per una distribuzione uniforme:
def moment( n, x, pdf ):
if ( n < 2 ):
print( "Errore! il primo parametro DEVE essere maggiore o uguale a 2!" )
return Null
else:
dx = x[1] - x[0]
f_ave = x * pdf
xm = dx * f_ave.sum()
f_ave = ( ( x - xm )**n ) * pdf
return dx * f_ave.sum()
N = 1000
xu = ran.random_sample( N )
Nc = 10
pdf, bins = np.histogram( xu, Nc, density = True )
nextr = len( bins )
# Calcola i centri dei bins
xc = 0.5 * ( bins[0:nextr-1] + bins[1:nextr] )
print( "Deviazione standard: ", moment( 2, xc, pdf ) )
print( "Valore teorico: ", 1.0 / 12.0 )
print( "Skewness: ", moment( 3, xc, pdf ) )
print( "Valore teorico: ", 0.0 )
print( "Flatness: ", moment( 4, xc, pdf ) )
print( "Valore teorico: ", 1.0 / 80.0 )
o per una distribuzione Gaussiana:
N = 1000
mu = 0.0
sigma = 1.0
xg = ran.normal( mu, sigma, N )
Nc = 50
pdf, bins = np.histogram( xg, Nc, density = True )
nextr = len( bins )
# Calcola i centri dei bins
xc = 0.5 * ( bins[0:nextr-1] + bins[1:nextr] )
print( "Deviazione standard: ", moment( 2, xc, pdf ) )
print( "Valore teorico: ", sigma )
print( "Skewness: ", moment( 3, xc, pdf ) )
print( "Valore teorico: ", 0.0 )
print( "Flatness: ", moment( 4, xc, pdf ) )
print( "Valore teorico: ", 3.0 * sigma**2 )
Nota che la libreria Scipy ha tutta una collezione di funzioni statistiche che consentono di calcolare i momenti delle funzioni di distribuzione di una serie di dati.
\[Y = min_{\alpha_1,\ldots,\alpha_n} \sum_i [y_i - f(x_i)]^2\]
Imponendo la condizione di minimo rispetto ai parametri della funzione \(\alpha_1,\ldots,\alpha_n\) si perviene ad un sistema di equazioni algebriche, purtroppo in generale non lineare, del tipo:
\begin{eqnarray*} \frac{\partial Y}{\partial \alpha_1} &=& 2 \sum_i [y_i - f(x_i)]\frac{\partial f}{\partial \alpha_1} = 0 \\ &\ldots& \\ \frac{\partial Y}{\partial \alpha_n} &=& 2 \sum_i [y_i - f(x_i)]\frac{\partial f}{\partial \alpha_n} = 0 \\ \end{eqnarray*}
che, una volta risolto, consente di determinare i parametri \(\alpha_j\) che minimizzano le differenze tra la funzione e i dati.
Supposto che il vettore dei dati x sia distribuito in maniera Gaussiana, trovare i due parametri \(\mu\) e \(\sigma\) della distribuzione.
In questo caso, \(n=2\) e \(\alpha_1 = \mu\) e \(\alpha_2 = \sigma\),
\[f(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]
che conviene scrivere in termini del logaritmo della \(f(x)\) come:
\[F(x; \mu,\sigma) = \log f(x; \mu, \sigma) = -\frac{1}{2} \log(2\pi) - \log\sigma - \frac{(x-\mu)^2}{2\sigma^2}\]
e, derivando parzialmente rispetto a \(\mu\) e \(\sigma\):
\begin{eqnarray*} \frac{\partial F}{\partial \mu} &=& \frac{x-\mu}{\sigma^2} \\ \frac{\partial F}{\partial \sigma} &=& -\frac{1}{\sigma} + \frac{(x-\mu)^2}{\sigma^3} \\ \end{eqnarray*}
che, inserito nelle equazioni, dà:
\begin{eqnarray*} \sum_i [\log y_i - F(x_i)] \frac{x_i-\mu}{\sigma^2} &=& 0 \\ \sum_i [\log y_i - F(x_i)] \left[ -\frac{1}{\sigma} + \frac{(x-\mu)^2}{\sigma^3} \right] &=& 0 \\ \end{eqnarray*}
che ci dà \(\mu\) e \(\sigma\) in termini di coefficienti che dipendono dalle somme su \(i\) di prodotti di quantità note: \(y_i\), \(x_i\), \(F(x_i)\), ecc.
Il sistema, in questo caso, è non lineare, quindi va risolto numericamente!
Il package scipy contiene un sotto-modulo optimize, che contiene una funzione per risolvere il problema in maniera semplice:
par, corr = curve_fit( f, xdata, ydata )
dove: par è il risultato, contenente il valore cercato dei parametri, corr è la matrice di correlazione tra i parametri, f è la funzione che descrive l'andamento dei dati e xdata e ydata sono i due vettori contenenti le \(x_i\) e le \(y_i\).
Esempio: best fit dei della distribuzione Gaussiana calcolata in precedenza:
import scipy.optimize
def F( x, mu, sigma ):
ff = -0.5 * np.log( 2.0 * np.pi ) - np.log( sigma ) - ( ( x - mu )**2.0 ) / ( 2.0 * sigma * sigma )
return ff
# Limita l'intervallo del fit ai valori per cui y_i != 0 per evitare problemi con il log!
imin = 2
imax = -5
xdata = xc[imin:imax]
ydata = np.log( pdf[imin:imax] )
par, corr = scipy.optimize.curve_fit( F, xdata, ydata )
print( "Parametri del fit [mu,sigma]: ", par )
Infine, plottiamo insieme la pdf e la curva per verificare la qualità del best-fit:
nextr = len( bins )
dx = bins[1] - bins[0]
plt.xlabel( r"$x$" )
plt.ylabel( "PDF" )
plt.bar( bins[:nextr-1], pdf, dx, color = "g" )
yc = np.exp( F( xc, par[0], par[1] ) )
plt.plot( xc, yc, "r-o" )
plt.show()
