固定拟合参数在curve_fit中



我有一个函数 Imaginary,它描述了物理过程,我想将其适合到数据集x_interpolate, y_interpolate。该函数是Lorentzian峰值函数的一种形式,除了f_peak(峰位置)外,我有一些用户给出的初始值,我使用峰查找算法找到了这些初始值。除偏移量外,所有拟合参数均应为正,因此我已经相应地设置了bounds_I

def Imaginary(freq, alpha, res, Ms, off):
    numerator = (2*alpha*freq*res**2)
    denominator = (4*(alpha*res*freq)**2) + (res**2 - freq**2)**2
    Im = Ms*(numerator/denominator) + off
    return Im
pI = np.array([alpha_init, f_peak, Ms_init, 0])
bounds_I = ([0,0,0,0, -np.inf], [np.inf,np.inf,np.inf, np.inf])
poptI, pcovI = curve_fit(Imaginary, x_interpolate, y_interpolate, pI, bounds=bounds_I)

在某些情况下,我想在拟合过程中保持参数f_peak修复。我通过将bounds_I更改为

来尝试简单的解决方案。
bounds_I = ([0,f_peak+0.001,0,0, -np.inf], [np.inf,f_peak-0.001,np.inf, np.inf])

这是出于多种原因而不是一种最佳的方式,所以我想知道是否有一种更重要的方式来做到这一点?谢谢您的帮助

如果固定参数,则不是真正的参数,因此应将其从参数列表中删除。定义一个模型,该模型被固定值替换并拟合。下面的示例,简化为简洁,要自我包容:

x = np.arange(10)
y = np.sqrt(x)    
def parabola(x, a, b, c):
  return a*x**2 + b*x + c
fit1 = curve_fit(parabola, x, y)  #  [-0.02989396,  0.56204598,  0.25337086]
b_fixed = 0.5
fit2 = curve_fit(lambda x, a, c: parabola(x, a, b_fixed, c), x, y) 

拟合的第二个调用返回 [-0.02350478, 0.35048631],这是a和c的最佳值。B的值固定为0.5。

当然,也应从初始向量pi和边界中删除参数。

您可能会发现lmfit(https://lmfit.github.io/lmfit-py/)有用。该库将更高级别的接口添加到Scipy优化例程中,旨在采用更优化的优化和曲线拟合方法。例如,它使用参数对象允许设置边界和固定参数,而无需修改目标或模型函数。对于曲线拟合,它定义了可以使用的高级模型函数。

为您示例,您可以使用

编写的Imaginary函数
from lmfit import Model
lmodel = Model(Imaginary)

然后创建参数(LMFit将根据您的函数签名命名参数对象),提供初始值:

params = lmodel.make_params(alpha=alpha_init, res=f_peak, Ms=Ms_init, off=0)

默认情况下,所有参数都是未结合的,并且在拟合中会有所不同,但是您可以修改这些属性(而无需重写模型函数):

params['alpha'].min = 0
params['res'].min = 0
params['Ms'].min = 0

您可以将一个(或更多)参数设置为与以下内容一样在拟合中不变。

params['res'].vary = False

要清楚:这不需要更改模型功能,使其更容易使用固定,可以强加的范围等等。

然后,您将使用模型和这些参数执行拟合:

result = lmodel.fit(y_interpolate, params, freq=x_interpolate)

您可以获取

参数的拟合统计信息,最佳拟合值和不确定性的报告
print(result.fit_report())

最佳拟合参数将在result.params中保存。

fwiw,lmfit还为许多常见形式(包括洛伦兹(Lorentzian)和持续偏移)设有内置模型。因此,您可以将此模型构建为

from lmfit.models import LorentzianModel, ConstantModel
mymodel = LorentzianModel(prefix='l_') + ConstantModel()
params = mymodel.make_params()

将具有名为l_amplitudel_centerl_sigmac的参数(其中c为常数),并且该模型将使用 x用于自变量(您的freq)。当您可能想更改峰或背景的功能形式时,或将多个峰拟合到光谱时,这种方法可能会变得非常方便。

我能够就任意数量的参数和固定参数的任意定位解决此问题:

def d_fit(x, y, param, boundMi, boundMx, listparam):
    Sparam, SboundMi, SboundMx = asarray([]), asarray([]), asarray([])
    Nparam, NboundMi, NboundMx = asarray([]), asarray([]), asarray([])
    for i in range(len(param)):
        if(listparam[i] == 1):
            Sparam = append(Sparam,asarray(param[i]))
            SboundMi = append(SboundMi,asarray(boundMi[i]))
            SboundMx = append(SboundMx,asarray(boundMx[i]))
        else:
            Nparam = append(Nparam,asarray(param[i]))
    def funF(x, Sparam):
        j = 0
        for i in range(len(param)):
            if(listparam[i] == 1):
                param[i] = Sparam[i-j]
            else:
                param[i] = Nparam[j]
                j = j + 1
        return fun(x, param)
    return curve_fit(lambda x, *Sparam: funF(x, Sparam), x, y, p0 = Sparam, bounds = (SboundMi,SboundMx)) 

在这种情况下:

param = [a,b,c,...] # parameters array (any size)
boundMi = [min_a, min_b, min_c,...] # minimum allowable value of each parameter
boundMx = [max_a, max_b, max_c,...] # maximum allowable value of each parameter
listparam = [0,1,1,0,...] # 1 = fit and 0 = fix the corresponding parameter in the fit routine

,根函数被定义为

def fun(x, param):
    a,b,c,d.... = param
    return a*b/c... # any function of the params a,b,c,d...

这样,您可以在不更改拟合例程的情况下更改根函数和参数数。而且,您可以随时修复或允许更改" ListParam"。

来拟合任何参数。

这样使用:

popt, pcov = d_fit(x, y, param, boundMi, boundMx, listparam)

" Popt"one_answers" pcov"是" 1"数量的1D数组。在" ListParam"中带来拟合参数的结果(最佳值和错误矩阵)

"参数"将对拟合值的原始(输入)"参数"(输入)的大小相同(输入)的"参数"的1D数组(输入)"参数"(paramquord; listParam"

希望可以很有用!

obs1:x = 1D阵列独立值,y = 1d阵列依赖值

obs2:这是我的第一篇文章。请让我知道我是否可以授予它!

最新更新