如何编写函数将数据拟合到 N 个类高斯峰值的总和,而无需明确定义每个可能的 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)

其中 Ae0hfSfwhm通过一些良好的初始猜测从拟合中确定。 重要的是,参数 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 的同时计算表达式并返回结果。

我尝试了各种各样的事情。我的主要障碍是我不知道如何告诉程序使用特定的 N 和相应的 i 值。最后,经过一番搜索,我想我找到了一种使用 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,n),x,y,p)
A,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,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,fwhm = fittedParameters

print(fittedParameters)
plt.plot(x,"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()

示例输出:

enter image description here

P.S.:从表面上看,这些数据点已经很好地拟合了 n=1