scipy 修复curve_fit中的拟合参数

irlmq6kh  于 2022-11-09  发布在  其他
关注(0)|答案(3)|浏览(243)

我有一个函数Imaginary,它描述了一个物理过程,我想将其拟合到数据集x_interpolate, y_interpolate。该函数是洛伦兹峰函数的一种形式,我有一些用户给定的初始值,除了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])

由于很多原因,这不是一个最佳的方式来做这件事,所以我想知道是否有一个更Python的方式来做这件事?谢谢你的帮助

o7jaxewo

o7jaxewo1#

如果某个参数是固定的,那么它就不是真正的参数,因此应该从参数列表中删除它。定义一个模型,该模型将该参数替换为固定值,并拟合该参数。下面的示例是为了简洁而简化的,并且是独立的:

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)

fit的第二个调用返回[-0.02350478, 0.35048631],这是a和c的最佳值。B的值固定为0.5。
当然,参数应该从初始向量pI和边界中移除。

gj3fmq9x

gj3fmq9x2#

您可能会发现lmfit(https://lmfit.github.io/lmfit-py/)很有帮助。这个库为scipy优化例程添加了一个更高级的接口,旨在提供一种更像Python的优化和曲线拟合方法。例如,它使用Parameter对象来允许设置边界和固定参数,而不必修改目标或模型函数。对于曲线拟合,它定义了可以使用的高级模型函数。
例如,您可以使用Imaginary函数,就像您用编写它一样

from lmfit import Model
lmodel = Model(Imaginary)

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

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作为自变量当你想要改变峰或背景的函数形式时,这种方法会变得非常方便。或者当将多个峰值拟合到光谱时,可以使用该方法。

xriantvc

xriantvc3#

我能够解决这个关于任意数量的参数和固定参数的任意定位的问题:

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...

这样,你就可以改变根函数和参数的数量,而不需要改变fit例程。而且,在任何时候,你都可以通过改变“listparam”来修复或让fit任何参数。
使用方法如下:

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

“popt”和“pcov”是“listparam”中数字“1”大小的1D数组,用于得出拟合参数的结果(最佳值和错误矩阵)
“param”将创建一个与原始(输入)“param”大小相同的1D数组,但它将自动更新为拟合值(与“popt”相同),并根据“listparam”保持固定值
希望能派上用场!
观察结果1:x = 1D-阵列独立值,y = 1D-阵列相关值
这是我的第一篇文章。如果我能改进的话,请告诉我!

相关问题