如何编写一个函数来将数据拟合到 N 个类高斯峰的总和,而无需为每个可能的 N 显式定义表达式?

How to write a function to fit data to a sum of N Gaussian-like peaks without explicitly defining the expression for every possible N?

我正在尝试将高斯峰的级数拟合到光谱线形。 级数是 N 均匀分布的高斯峰的总和。当编码为函数时,N=1 的公式如下所示:

A * ((e0-i*hf)/e0)**3 * ((S**i)/np.math.factorial(i)) * np.exp(-4*np.log(2)*((x-e0+i*hf)/fwhm)**2)

where A, e0, hf, Sfwhm 将通过一些良好的初始猜测从拟合中确定。 重要的是,参数 i 从 0 开始,每增加一个组件就递增 1。

因此,对于 N = 3,表达式将采用以下形式:

A * ((e0-0*hf)/e0)**3 * ((S**0)/np.math.factorial(0)) * np.exp(-4*np.log(2)*((x-e0+0*hf)/fwhm)**2) +
A * ((e0-1*hf)/e0)**3 * ((S**1)/np.math.factorial(1)) * np.exp(-4*np.log(2)*((x-e0+1*hf)/fwhm)**2) +
A * ((e0-2*hf)/e0)**3 * ((S**2)/np.math.factorial(2)) * np.exp(-4*np.log(2)*((x-e0+2*hf)/fwhm)**2)

除了i之外的所有参数对于求和中的每个分量都是常数,这是有意的。 i 正在根据参数的数量以受控方式变化。

我正在使用 curve_fit。编写拟合例程的一种方法是为任何合理的 N 显式定义表达式,并只使用适当的表达式。就像,这里是 5 或 6,取决于间距,间距由 hf 决定。我可以只定义一个包含 N 个组件的长函数,将适当的 i 值写入每个组件。我知道如何做到这一点(并且做到了)。但我想更智能地编写代码。我的目标是编写一个接受任何 N 值的函数,如上所述添加适当数量的组件,计算表达式,同时适当地递增 i 和 return 结果。

我尝试了很多事情。我的主要障碍是我不知道如何告诉程序使用特定的 Ni 的相应值。最后,经过一番搜索,我认为我找到了一种使用 lambda 函数对其进行编码的好方法。

from scipy.optimize import curve_fit
import numpy as np

def fullfunc(x,p,n):
    def func(x,A,e0,hf,S,fwhm,i):
        return A * ((e0-i*hf)/e0)**3 * ((S**i)/np.math.factorial(i)) * np.exp(-4*np.log(2)*((x-e0+i*hf)/fwhm)**2)
    y_fit = np.zeros_like(x)
    for i in range(n):
        y_fit += func(x,p[0],p[1],p[2],p[3],p[4],i)
    return y_fit
p = [1,26000,1400,1,1000]
x = [27027,25062,23364,21881,20576,19417,18382,17452,16611,15847,15151]
y = [0.01,0.42,0.93,0.97,0.65,0.33,0.14,0.06,0.02,0.01,0.004]
n = 7
fittedParameters, pcov = curve_fit(lambda x,p: fullfunc(x,p,n), x, y, p)
A,e0,hf,S,fwhm = fittedParameters

这给出:

TypeError: <lambda>() takes 2 positional arguments but 7 were given

我不明白为什么。我觉得 lambda 函数无法处理初始参数列表。 我将不胜感激关于如何在不明确写出所有方程式的情况下如何进行这项工作的任何建议,因为我发现这有点太死板了。 提供的 x 和 y 范围是真实数据的样本,可以大致了解形状是什么。

由于您仅在 i=0, 1, ..., n-1 范围内使用求和,因此无需参考复杂的 lambda 结构,这些结构在曲线拟合的上下文中可能有效,也可能无效。只需将拟合函数定义为 n 个分量的总和:

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

def func(x, A, e0, hf, S, fwhm):
    return sum((A * ((e0-i*hf)/e0)**3 * ((S**i)/np.math.factorial(i)) * np.exp(-4*np.log(2)*((x-e0+i*hf)/fwhm)**2)) for i in range(n))

p = [1,26000,1400,1,1000]
x = [27027,25062,23364,21881,20576,19417,18382,17452,16611,15847,15151]
y = [0.01,0.42,0.93,0.97,0.65,0.33,0.14,0.06,0.02,0.01,0.004]
n = 7
fittedParameters, pcov = curve_fit(func, x, y, p0=p)
#A,e0,hf,S,fwhm = fittedParameters

print(fittedParameters)
plt.plot(x, y, "ro", label="data")
x_fit = np.linspace(min(x), max(x), 100)
y_fit = func(x_fit, *fittedParameters)
plt.plot(x_fit, y_fit, label="fit")
plt.legend()
plt.show()

示例输出:

P.S.: 看起来,这些数据点已经很好地拟合了n=1.