Ajustando una distribución normal
15 Dec 2022Como el teorema del límite central es tan común en la naturaleza y aparece con tanta frecuencia cuando se levantan grandes cantidades de datos disponibles en nuestra era, que tener una función para ajustar una distribución gaussiana o normal no es una mala idea.
Me explico: supongamos que tenemos un set de datos que es de nuestro interés, y al visualizarlo notamos que podría ser descrito por un perfíl gaussiano. ¿Cómo obtenemos o extraemos los valores de dicha distribución, con el fin de estudiar sus parámetros?
Creemos una distribución artificial y veamos cómo entender dicha distribución, usando un puñado de puntos (100 mediciones) para poder caracterizar la muestra:
from numpy.random import seed
from numpy.random import normal
import matplotlib.pyplot as plt
import numpy as np
#celda magica
%matplotlib inline
#semilla para reproduccion
seed(40)
#generamos datos con media cero, dispersion 1 y de tamaño 100
centro = 0
desviacion = 1
plt.figure(figsize=(6,3))
data = normal(loc=centro, scale=desviacion, size=100)
_,_,_ = plt.hist(data, 18, label='distribucion generada')
plt.legend(loc='best')
plt.show()
Estos datos podrían corresponder a variados fenómenos físicos y/o describir el funcionamiento de algún proceso conocido de antemano. Supongamos que son reales y es imperante entenderlos y caracterizarlo adecuadamente.
Entonces, usaremos la librería scipy, particularmente la función optimize para crear una función que ajuste un perfil gaussiano descrito explicitamente, y los restos (la comparación entre la distribución y el modelo gaussiano) serán minimizados en el sentido de los mínimos cuadrados. Esto será logrado (ajustado o fiteado) usando una función de bondad de ajuste llamada chi cuadrado $\chi ^2$.
Vamos a ello:
def fit_gauss(X,guess,plot=False,large=[6,4],normed=True,xlabel='',Nbins=None,lw=3,fontsize=15,points=False):
"""
Ajusta una función gaussiana a una distribucion dada.
Entrada:
X: Distribucion.
guess: Valores iniciales para realizar el ajuste (guess=[mu,sigma,amplitud])
plot: mostrar el ajuste.
large: tamaño del plot.
normed: Si se normaliza el histograma y la funcion.
Salida:
mu: Valor medio de la gaussiana.
sigma: desviacion estandar del ajuste.
amplitud: amplitud de la funcion.
"""
from scipy import optimize
import numpy as np
import matplotlib.pyplot as plt
#Entendemos la distribución de entrada:
(altura, bins) = np.histogram(X, bins=Nbins, density=normed)
hist = np.array([ (bins[i]+bins[i+1])/2. for i in range(len(bins)-1)])
#describimos explicitamente la función gaussiana:
def gaussian(val, x):
return (val[2])*np.exp(-0.5*((x - val[0])/val[1])**2)
#describimos explicitamente la métrica de bondad de ajuste, chi cuadrado (chi2)
def chi2(val,x,y):
return (y-gaussian(val, x))**2 #Asumiendo que los errores son 1
x = np.linspace(np.min(X),np.max(X),Nbins)
xi = np.linspace(np.min(X),np.max(X),1000)
#realizamos el ajuste usando los minimos cuadrados.
fit = optimize.leastsq(chi2,guess,args=(hist,altura))
#obtenemos los parámetros calculados.
mu = fit[0][0]
sigma = abs(fit[0][1])
ampl = abs(fit[0][2])
#agregamos algunas funciones de visualización:
if plot==True:
plt.figure(figsize=large)
if normed == True:
plt.ylabel(r'${\rm Normalized\ Amplitude }$',fontsize=fontsize)
else:
plt.ylabel(r'${\rm Amplitude }$',fontsize=fontsize)
plt.xlabel(r'${\rm %s }$' %xlabel,fontsize=fontsize)
a,b,c=plt.hist(X,Nbins,density=normed, label='distribucion generada')
if points==True:
plt.plot(hist,altura,'o--',color='black', label='datos ajustados')
plt.plot(xi,gaussian(fit[0],xi),'-',color='red',linewidth=lw, label=' modelo gaussiano')
#listoco
return mu,sigma,ampl
Para ayudar a la convergencia de nuestro fit, "adivinaremos" los parámetros que podrían llegar a describir la distribución gaussiana, esto es normalmente llamado "guess". Usaremos los promedios de la distribución para obtener el $\mu$, la desviación estandar para obtener el $\sigma$, y la altura es el parámetro más laxo, usualmente lo fijo en 10 sin ningun análisis previo.
#Obtenemos el guess
guess = [np.mean(data), np.std(data),10]
#Aplicamos la función y obtenemos los parámetros de la distribución:
mu,sigma,ampl = fit_gauss(data,guess=guess,plot=True,large=[8,4],normed=True,xlabel='distribucion\ generada',Nbins=18, points=True)
#Definimos el factor de threshold:
factor = 2
#preparamos el gráfico y ploteamos:
plt.axvline(x= mu, color='red',)
plt.axvline(x= mu+factor*sigma)
plt.text( mu+factor*sigma,0.33,'$\mu+2\sigma$',rotation=90,fontsize=20)
plt.axvline(x= mu-factor*sigma)
plt.text( mu-factor*sigma - 0.27 ,0.33,'$\mu-2\sigma$',rotation=90,fontsize=20)
plt.legend(loc='best', fontsize=6)
print('Valor mu original: ',centro)
print('Valor sigma original: ',desviacion)
print('Valor mu extraido: ',mu)
print('Valor sigma extraido: ',sigma)
print('Amplitud gaussiana: ',ampl)
plt.show()
Entonces, desde la función hemos obtenido parámetros similares desde los que fueron creados los datos sintéticos desde una media $\mu=0$ ($\mu_{fit} = 0.06$), con $\sigma=1$ ($\sigma_{fit} = 0.95$), y una altura de 0.39.
Se ve que los puntos negros que describen cada bin por separado, difieren bastante del modelo gaussiado utilizado (en rojo). Aún así, en promedio, los restos se minimizan y el ajuste se logra con cierto éxito. De hecho, es comprobable (¿tarea para el/la revisor/a?) que la distribución de los restos entre los datos y el modelo generan una nueva distribución gaussiana, y es por esto que el fit funciona. De hecho si los restos no fueran gaussianos, me preocuparía mucho el fit, independiente del aspecto que tenga.
Estas funciones se podrían hacer recursivas para caracterizar distribuciones aún más complejas dependiendo de los requerimientos, y por supuesto que, al generar datos con más puntos, la estimación será aun más precisa.
Veamos esto y finalicemos el caso de uso:
#usamos diferentes parámetros de entrada sólo para comprobar la función:
centro = 43245.567
desviacion = 200
Ndatos = 10000
data = normal(loc=centro, scale=desviacion, size=Ndatos)
#adivinamos dónde están, aproximadamente, los valores de ajuste.
guess = [np.mean(data), np.std(data),10]
#Aplicamos la función y obtenemos los parámetros de la distribución:
mu,sigma,ampl = fit_gauss(data,guess=guess,plot=True,large=[8,3],normed=True,xlabel='distribucion\ generada',Nbins=30, points=True)
plt.legend(loc='best')
print('Valor mu original: ',centro)
print('Valor sigma original: ',desviacion)
print('Valor mu extraido: ',mu)
print('Valor sigma extraido: ',sigma)
print('Amplitud gaussiana: ',ampl)
plt.show()
Bueno, y si hasta aquí ya no resolviste tu duda, quizá lo que buscabas era una forma de tomar una distribución dada, y luego probar diferentes modelos, cosa que lo haga rápido y sin tantas preguntas, que me lo haga ahora yaaaaaa.
Se pueden usar librerias especializadas, como Fitter, veamos un ejemplo.
from fitter import Fitter
from fitter import get_common_distributions
plt.figure(figsize = (10,3))
modelos = Fitter(data, distributions=get_common_distributions())
#Ejecutamos el analisis
modelos.fit()
modelos.summary()
Vemos que existe un sumario visual de las estadísticas y modelos posibles de generación de una distribución. Los mejores indicadores son una distribución normal, para sorpresa de nadie, pero no se ve porque está debajo del siguiente modelo que ajusta a la perfección: la distribución gamma. Interesante, podríamos indagar en esto, pero lástima que está fuera de mis intereses justo ahora. Pero antes de irnos, asegurémoos de que los parámetros son correctos:
#Usando el método de minimización de la diferencia cuadrada de los errores y obtenemos el mejor modelo
f.get_best(method = 'sumsquare_error')
Como vemos, el mejor modelo es la distribución normal (norm), usando una media (loc) de $\mu=43243.15$ [alguna unidad], y un sigma (scale) de $\sigma=201.27$. Me parece suficiente para mis propósitos.
Gracias por leer!.