如何优化python中高斯-厄米特函数的拟合？

Question

到目前为止我做了什么：

我正在尝试将噪声数据（我通过向我的函数添加随机噪声自己生成）拟合到我定义的 Gauss-Hermite 函数。在某些情况下，它适用于较低的 h3 和 h4 值，但每隔一段时间，即使对于较低的 h3、h4 值，它也会产生非常糟糕的拟合，而对于较高的 h3、h4 值，它总是会产生不好的拟合。

我的代码：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as mpl
import matplotlib.pyplot as plt

# Let's define the Gauss-Hermite function
#a=amplitude, x0=location of peak, sig=std dev, h3, h4
def gh_func(x, a, x0, sig, h3, h4):
   
    return a*np.exp(-.5*((x-x0)/sig)**2)*(1+h3*((-np.sqrt(3))*((x-x0)/sig)+(2/np.sqrt(3))*((x-x0)/sig)**3)+ h4*((np.sqrt(6)/4)+(-np.sqrt(6))*((x-x0)/sig)**2+(np.sqrt(6)/3)*(((x-x0)/sig)**4)))

#generate clean data
x = np.linspace(-10, 20, 100)
y = gh_func(x, 10, 5, np.sqrt(3), -0.10,-0.03) #it gives okay fit for h3=-0.10, h4=-0.03 but bad fits for higher values like h3=-0.4 and h4=-0.3.

#add noise to data

noise=np.random.normal(0,np.sqrt(0.5),size=len(x))
yn = y + noise

fig = mpl.figure(1)
ax = fig.add_subplot(111)
ax.plot(x, y, c='k', label='analytic function')
ax.scatter(x, yn, s=5, label='fake noisy data')
fig.savefig('model_and_noise_h3h4.png')
# Executing curve_fit on noisy data
popt, pcov = curve_fit(gh_func, x, yn)

#popt returns the best fit values for parameters of the given model (func)


print('Fitted Parameters (Gaus_Hermite):\na = %.10f , x0 = %.10f , sig = %.10f\nh3 = %.10f , h4 = %.10f' \
%(popt[0],popt[1],popt[2],popt[3],popt[4]))

ym = gh_func(x, popt[0], popt[1], popt[2], popt[3], popt[4])
ax.plot(x, ym, c='r', label='Best fit')
ax.legend()
fig.savefig('model_fit_h3h4.png')
plt.legend(loc='upper left')
plt.xlabel("v")
plt.ylabel("f(v)")

我想做什么：

我想找到比 scipy.optimize 中的 curve_fit 更好的拟合方法，但我不确定我可以使用什么。即使我们最终使用curve_fit，我也需要一种方法来通过为自动生成的参数提供初始猜测来产生更好的拟合，例如这篇文章的公认答案（Jean Jacquelin 的方法）描述了一种仅用于单峰高斯的方法：gaussian fitting inaccurate for lower peak width using Python。但这仅适用于 mu、sigma 和振幅而不是 h3、h4。

除了来自 scipy.optimize 的 curve_fit，我认为还有一个叫做 lmfit:https://lmfit.github.io/lmfit-py/ 但我不确定我将如何在我的代码中实现它。我不想对参数使用手动初始猜测。我希望能够自动找到配件。

非常感谢！

Answer 1

使用lmfit 解决这个拟合问题会很简单，通过从您的gh_func 创建一个lmfit.Model，类似

from lmfit import Model
gh_model = Model(gh_func)
params = gh_model.make_params(x0=?, a=?, sig=?, h3=?, h4=?)

这些问号必须是相关参数的初始值。

您对scipy.optimize.curve_fit 的使用不提供变量的初始值。不幸的是，这不会引发错误或指出问题，因为scipy.optimize.curve_fit 的作者通过将初始值设为可选来欺骗您。所有非线性最小二乘分析总是需要所有参数的初始值。它不是实现的选择，而是数学算法的一个特征。 curve_fit 对你隐藏的是，离开 p0=None 会使所有初始值变为 1.0。这是否合适取决于模型和数据是否合适——它并不总是合理的。例如，如果您的数据的 x 值从 9000 扩展到 9500，并且峰值函数以 9200 为中心，那么从 1 开始的 x0 几乎肯定不会找到合适的拟合。以任何方式暗示初始值是可选的总是具有欺骗性的。他们只是不是。有时不是，不是针对明确定义的问题。初始值永远是可选的。

初始值不一定是完美的，但它们必须具有正确的比例。许多人会警告您“错误最小值”——陷入一个不错但不是“最佳”的解决方案。我的经验是，人们遇到的更常见的问题是初始值相差太远，以至于模型对其值的微小变化不敏感，因此永远无法优化（x0=1,sig=1 用于 x on [ 9000, 9500] 并且以 9200 为中心正好属于该类别）。

如果不知道模型函数的属性，很难回答“为自动生成的参数提供初始猜测”的基本问题。您也许可以使用一些启发式方法来猜测数据集中的值。 Lmfit 对于许多类似峰值的函数具有这样的启发式“从数据中猜测参数值”函数。您可能对h3 和h4 的物理（或其他领域）有所了解，并且知道哪些类型的值是合理的，并且可能会给出比h3=h4=1 更好的初始值。您可能想从猜测参数开始，就好像数据是简单的高斯分布（例如，使用 lmfit.Model.guess_from_peak()），然后使用它与您的数据的差异来获得 h3 和 h4 的比例。