Ajustar una curva a un histograma en Python

Contamos con el resultado a esta interrogante, o por lo menos eso pensamos. Si continuas con preguntas dínoslo, que sin dudarlo te ayudaremos

Solución:

De la documentación de matplotlib.pyplot.hist:

Devoluciones

n : array o lista de arreglos

Los valores de los intervalos de histograma. Ver normed y weights para una descripción de la posible semántica. Si entrada x es un arrayentonces este es un array de longitud nbins. Si la entrada es una matriz de secuencias [data1, data2,..]entonces esta es una lista de arreglos con los valores de los histogramas para cada uno de los arreglos en el mismo orden.

contenedores: array

Los bordes de los contenedores. Longitud nbins + 1 (nbins borde izquierdo y borde derecho del último bin). siempre un solo array incluso cuando se pasan varios conjuntos de datos.

parches: lista o lista de listas

Lista silenciosa de parches individuales utilizados para crear el histograma o lista de dicha lista si hay múltiples conjuntos de datos de entrada.

Como puede ver, el segundo retorno son en realidad los bordes de los contenedores, por lo que contiene un artículo más que contenedores.

La forma más fácil de obtener los centros de contenedores es:

import numpy as np
bin_center = bin_borders[:-1] + np.diff(bin_borders) / 2

Lo que solo agrega la mitad del ancho (con np.diff) entre dos bordes (ancho de los contenedores) hasta el borde izquierdo del contenedor. Excluyendo el borde del último contenedor porque es el borde derecho del contenedor más a la derecha.

Así que esto devolverá los centros de los bins – un array con la misma longitud que n.

Tenga en cuenta que si tiene numba, podría acelerar el cálculo de bordes a centros:

import numba as nb

@nb.njit
def centers_from_borders_numba(b):
    centers = np.empty(b.size - 1, np.float64)
    for idx in range(b.size - 1):
        centers[idx] = b[idx] + (b[idx+1] - b[idx]) / 2
    return centers

def centers_from_borders(borders):
    return borders[:-1] + np.diff(borders) / 2

Es bastante más rápido:

bins = np.random.random(100000)
bins.sort()

# Make sure they are identical
np.testing.assert_array_equal(centers_from_borders_numba(bins), centers_from_borders(bins))

# Compare the timings
%timeit centers_from_borders_numba(bins)
# 36.9 µs ± 275 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit centers_from_borders(bins)
# 150 µs ± 704 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Incluso si es más rápido, numba es una dependencia bastante fuerte que no se agrega a la ligera. Sin embargo, es divertido jugar con él y es realmente rápido, pero a continuación usaré la versión NumPy porque será más útil para la mayoría de los futuros visitantes.

---

En cuanto a la tarea general de ajustar una función al histograma: debe definir una función que se ajuste a los datos y luego puede usar scipy.optimize.curve_fit. Por ejemplo, si desea ajustar una curva de Gauss:

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

Luego defina la función para ajustar y algún conjunto de datos de muestra. El conjunto de datos de muestra es solo para el propósito de esta pregunta, debe usar su conjunto de datos y definir la función que desea ajustar:

def gaussian(x, mean, amplitude, standard_deviation):
    return amplitude * np.exp( - ((x - mean) / standard_deviation) ** 2)

x = np.random.normal(10, 5, size=10000)

Ajustar la curva y trazarla:

bin_heights, bin_borders, _ = plt.hist(x, bins='auto', label='histogram')
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2
popt, _ = curve_fit(gaussian, bin_centers, bin_heights, p0=[1., 0., 1.])

x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 10000)
plt.plot(x_interval_for_fit, gaussian(x_interval_for_fit, *popt), label='fit')
plt.legend()

Tenga en cuenta que también puede usar NumPys histogram y Matplotlibs bar-parcela en su lugar. la diferencia es que np.histogram no devuelve los “parches” array y que necesita los anchos de bin para el diagrama de barras de Matplotlibs:

bin_heights, bin_borders = np.histogram(x, bins='auto')
bin_widths = np.diff(bin_borders)
bin_centers = bin_borders[:-1] + bin_widths / 2
popt, _ = curve_fit(gaussian, bin_centers, bin_heights, p0=[1., 0., 1.])

x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 10000)

plt.bar(bin_centers, bin_heights, width=bin_widths, label='histogram')
plt.plot(x_interval_for_fit, gaussian(x_interval_for_fit, *popt), label='fit', c='red')
plt.legend()

Por supuesto, también puede adaptar otras funciones a sus histogramas. En general, me gustan los modelos de Astropys para el ajuste, porque no necesita crear las funciones usted mismo y también admite modelos compuestos y diferentes ajustadores.

Por ejemplo, para ajustar una curva de Gauss usando Astropy al conjunto de datos:

from astropy.modeling import models, fitting

bin_heights, bin_borders = np.histogram(x, bins='auto')
bin_widths = np.diff(bin_borders)
bin_centers = bin_borders[:-1] + bin_widths / 2

t_init = models.Gaussian1D()
fit_t = fitting.LevMarLSQFitter()
t = fit_t(t_init, bin_centers, bin_heights)

x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 10000)
plt.figure()
plt.bar(bin_centers, bin_heights, width=bin_widths, label='histogram')
plt.plot(x_interval_for_fit, t(x_interval_for_fit), label='fit', c='red')
plt.legend()

Es posible ajustar un modelo diferente a los datos simplemente reemplazando:

t_init = models.Gaussian1D()

con otro modelo. por ejemplo un Lorentz1D (como un gaussiano pero con colas más anchas):

t_init = models.Lorentz1D()

No es exactamente un buen modelo dados mis datos de muestra, pero es realmente fácil de usar si ya hay un modelo de Astropy que se ajusta a las necesidades.

---

---

Utiliza Nuestro Buscador