曲线与自定义功能的直方图不符合

Question

我试图将曲线拟合到直方图，但是即使直方图不是直方图，所得曲线也是平坦的。如何正确安装曲线？

我当前的代码：

import numpy as np

import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

import pandas as pd

import scipy.optimize as optimization


x = np.random.normal(1e-10, 1e-7, size=10000)

def func(x, a, b, c):
    return a * (np.exp(-b*(x-c)**2))

bins=25

logbins = np.logspace(np.log10(1.0E-10),np.log10(1E-07),bins)


bin_heights, bin_borders, _ = plt.hist(x, bins=logbins, edgecolor='black', color='b')
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

x0    = np.array([0.0, 0.0, 0.0])

popt,cov = optimization.curve_fit(func, bin_centers, bin_heights,x0)
a,b,c=popt

popt, pcov = curve_fit(func, bin_centers, bin_heights, p0=[a,b,c])

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

plt.plot(x_interval_for_fit, func(x_interval_for_fit, *popt), label='Fit',color='r')


plt.xscale('log')

plt.show()

结果：

Answer 1

您会获得不良结果，因为您使用的函数适合直方图根本不像直方图的形状。通过使用简单的二阶插值函数，结果要好得多（尽管您可能会说不理想）：

def func(x, a, b, c):
    return a * x**2 + b * x + c  # Simple 2nd-order polynomial

与代码一起使用（您可以删除两个优化步骤，并且只能执行一次），我得到了以下结果：

您代码中可能无意识的一件事是，尽管您是事实，但在直方图上创建了一个正态分布，您决定以令人惊讶的方式汇总它们：仅考虑您在分布的一侧的东西（因为您从1e-10开始，并且您的分布是中心在1E-10中），然后将垃圾桶大小对数增加，最终会在较大的垃圾箱上获得更多点。另外，您忽略了一半以上的积分（这些要小于1E-10，检查hist's 文档）。您可以通过：np.sum（bin_heights）检查它，您会发现计数不到len（x）的一半。

如果我上面所说的是无意的，您真的想找出随机数的原始生成高斯（但使用计数而不是密度），那么您应该直接在未修改的直方图上做到这一点。但是，一个问题，因为您使用的是非常小的标准偏差，所以问题的条件非常糟糕（从某种意义上说，雅各布的价值差异很大，具体取决于您走到哪个维度）。您可以感觉到优化者通过手动计算错误并在改变参数时查看其更改（或不改变）来理解每个参数的影响是多么困难。由于您正在进行最小二乘优化，因此您的成本函数看起来像这样：

f = lambda a, b, c: sum((bin_heights - func(bin_centers, a, b, c))**2)

如果您使用此功能稍微播放一些初始条件，请查看优化者所看到的，您会注意到为什么会融合很难：

f(0, 0, 0)   # 8_407_566
f(0, 0, 1)   # 8_407_566, i.e. no change
f(0, 1, 0)   # 8_407_566, i.e. no change
f(1, 0, 0)   # 8_387_591, clearly an improvement

这不是数学上严格的论点（a和b接近0时，雅各布的不良条件均接近0），但它很好地直觉了。

因此，您需要一个很好的“第一个猜测”。当然，在这个问题中，由于您知道自己正在尝试适合高斯，因此可以将平均值，最大值和标准偏差用作初始参数（在您的情况下，是标准偏差正方形的倒数）。下面的代码显示了一组“较少调谐”的初始参数，这些参数也有效：

def func(x, a, b, c):  # Original function
    return a * np.exp(- b * (x - c)**2)

bins=25  # Let hist do its magic
bin_heights, bin_borders, _ = plt.hist(x, bins=bins, edgecolor='black', color='b')
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

x0 = np.array([600.0, 1.0e7, 0.0])

popt, cov = curve_fit(func, bin_centers, bin_heights, x0)

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

plt.plot(x_interval_for_fit, func(x_interval_for_fit, *x0), label='Initial guess',color='y')
plt.plot(x_interval_for_fit, func(x_interval_for_fit, *popt), label='Fit',color='r')
plt.show()

结果，结果是由初始猜测（以黄色）创建的曲线以及插值的结果（红色）：