在Python中将曲线拟合到直方图

Question

我试图使曲线适合matplotlib生成的直方图中的值:

n, bins, patches = plt.hist(myData)

其中"plt"代表matplotlib.pyplot,myData是一个数组,每个索引都出现次数,如[9,3,3,....]

我希望箱子是我的x数据,n是我的y数据.也就是说,我想提取有关x号码与号码x的频率的信息.但是,我不能让箱子和n具有相同的尺寸.

所以基本上,我希望能够将曲线拟合到n(bins,params).

怎么会这样做？

Answer 1

来自以下文件matplotlib.pyplot.hist:

返回

n:数组或数组列表

直方图箱的值.请参阅normed并weights查看可能的语义.如果input x是一个数组,那么这是一个长度数组nbins.如果输入是序列数组[data1, data2,..],那么这是一个数组列表,其中每个数组的直方图值的顺序相同.

箱子:阵列

箱子的边缘.长度nbins + 1(nbins左边缘和最后一个bin的右边缘).即使传入多个数据集,也始终是一个数组.

补丁:列表或列表列表

如果有多个输入数据集,则用于创建直方图或此列表列表的各个修补程序的静默列表.

正如您所看到的那样,第二个返回实际上是垃圾箱的边缘,因此它包含的物品比垃圾箱多一个.

获取垃圾箱中心的最简单方法是:

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

这只是将np.diff两个边框(箱子的宽度)之间的宽度(有)的一半加到左边框边框.排除最后一个bin边框,因为它是最右边bin的右边框.

所以这实际上将返回bin中心 - 一个长度相同的数组n.

请注意,如果你有numba,你可以加速边界到中心的计算:

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

它快得多:

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)

即使它是更快的numba是一个非常重的依赖,你不要轻易添加.然而,玩起来很快很有趣,但在下面我将使用NumPy版本,因为它对大多数未来的访问者会更有帮助.

至于将函数拟合到直方图的一般任务:您需要定义一个适合数据的函数,然后才能使用scipy.optimize.curve_fit.例如,如果您想要拟合高斯曲线:

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

然后定义要拟合的函数和一些样本数据集.样本数据集仅用于此问题,您应该使用数据集并定义要拟合的函数:

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

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

拟合曲线并绘制曲线:

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()

请注意,您也可以使用NumPys histogram和Matplotlibs bar-plot.区别在于np.histogram不返回"补丁"数组,并且您需要Matplotlibs条形图的bin宽度:

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()

当然,您也可以在直方图中使用其他功能.我一般都喜欢Astropy的拟合模型,因为你不需要自己创建函数,它也支持复合模型和不同的拟合.

例如,使用Astropy将高斯曲线拟合到数据集:

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()

只需更换以下内容即可为数据拟合不同的模型:

t_init = models.Gaussian1D()

与不同的模型.例如a Lorentz1D(像高斯但有更宽的尾巴):

t_init = models.Lorentz1D()

根据我的样本数据,这不是一个好的模型,但如果已经有一个符合需求的Astropy模型,那么它真的很容易使用.