Python: ajuste gaussiano

¿Qué es la distribución normal o gaussiana?

Cuando graficamos un conjunto de datos como un histograma, la forma de ese gráfico graficado es lo que llamamos su distribución. La forma más comúnmente observada de valores continuos es la curva de campana, también llamada distribución gaussiana o normal.

Lleva el nombre del matemático alemán Carl Friedrich Gauss. Algunos conjuntos de datos de ejemplo comunes que siguen la distribución gaussiana son la temperatura corporal, la altura de las personas, el kilometraje del automóvil y las puntuaciones de coeficiente intelectual. 

Intentemos generar la distribución normal ideal y trazarla usando Python.

Cómo trazar la distribución gaussiana en Python

Tenemos bibliotecas como Numpy, scipy y matplotlib para ayudarnos a trazar una curva normal ideal.

import numpy as np
import scipy as sp
from scipy import stats
import matplotlib.pyplot as plt 
  
## generate the data and plot it for an ideal normal curve
  
## x-axis for the plot
x_data = np.arange(-5, 5, 0.001)
  
## y-axis as the gaussian
y_data = stats.norm.pdf(x_data, 0, 1)
  
## plot data
plt.plot(x_data, y_data)

Producción:

Los puntos en el eje x son las observaciones, y el eje y es la probabilidad de cada observación.

Generamos observaciones espaciadas regularmente en el rango (-5, 5) usando np.arange() . Luego lo ejecutamos a través de la función norm.pdf() con una media de 0,0 y una desviación estándar de 1, que arrojó la probabilidad de esa observación. Las observaciones alrededor de 0 son las más comunes, y las alrededor de -5.0 y 5.0 son raras. El término técnico para la función pdf() es función de densidad de probabilidad.

La función de Gauss:

Primero, ajustemos los datos a la función gaussiana. Nuestro objetivo es encontrar los valores de A y B que mejor se ajusten a nuestros datos. Primero, necesitamos escribir una función de Python para la ecuación de la función Gaussiana. La función debe aceptar la variable independiente (los valores de x) y todos los parámetros que la componen.

#Define the Gaussian function
def gauss(x, H, A, x0, sigma):
    return H + A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

Usaremos la función curve_fit del módulo de Python scipy.optimize para ajustar nuestros datos. Utiliza mínimos cuadrados no lineales para ajustar los datos a una forma funcional. Puede obtener más información sobre curve_fit utilizando la función de ayuda dentro del cuaderno de Jupyter Notebook o la documentación en línea de Scipy .

La función curve_fit tiene tres entradas requeridas: la función que desea ajustar, los datos x y los datos y que ajusta. Hay dos salidas. El primero es una array de los valores óptimos de los parámetros. La segunda es una array de la covarianza estimada de los parámetros a partir de la cual puede calcular el error estándar de los parámetros.

Ejemplo 1:

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
xdata = [ -10.0, -9.0, -8.0, -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
ydata = [1.2, 4.2, 6.7, 8.3, 10.6, 11.7, 13.5, 14.5, 15.7, 16.1, 16.6, 16.0, 15.4, 14.4, 14.2, 12.7, 10.3, 8.6, 6.1, 3.9, 2.1]
  
# Recast xdata and ydata into numpy arrays so we can use their handy features
xdata = np.asarray(xdata)
ydata = np.asarray(ydata)
plt.plot(xdata, ydata, 'o')
  
# Define the Gaussian function
def Gauss(x, A, B):
    y = A*np.exp(-1*B*x**2)
    return y
parameters, covariance = curve_fit(Gauss, xdata, ydata)
  
fit_A = parameters[0]
fit_B = parameters[1]
  
fit_y = Gauss(xdata, fit_A, fit_B)
plt.plot(xdata, ydata, 'o', label='data')
plt.plot(xdata, fit_y, '-', label='fit')
plt.legend()

Ejemplo 2:

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as mpl
  
# Let's create a function to model and create data
def func(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
  
# Generating clean data
x = np.linspace(0, 10, 100)
y = func(x, 1, 5, 2)
  
# Adding noise to the data
yn = y + 0.2 * np.random.normal(size=len(x))
  
# Plot out the current state of the data and model
fig = mpl.figure()
ax = fig.add_subplot(111)
ax.plot(x, y, c='k', label='Function')
ax.scatter(x, yn)
  
# Executing curve_fit on noisy data
popt, pcov = curve_fit(func, x, yn)
  
#popt returns the best fit values for parameters of the given model (func)
print (popt)
  
ym = func(x, popt[0], popt[1], popt[2])
ax.plot(x, ym, c='r', label='Best fit')
ax.legend()
fig.savefig('model_fit.png')

Producción: