使用 lmfit 拟合多个数据集而不编写目标函数

Question

本主题描述如何使用 lmfit 拟合多个数据集： Python and lmfit: How to fit multiple datasets with shared parameters?

但是它使用用户编写的拟合/目标函数。

我想知道是否可以使用 lmfit 拟合多个数据集，而无需编写目标函数并使用模型类的 model.fit() 方法。

例如：假设我们有多个 (x,y) 坐标数据集，我们希望使用相同的模型函数来拟合这些数据集，以便找到平均而言最适合所有数据的参数集。

import numpy as np 
from lmfit import Model, Parameters
from lmfit.models import GaussianModel

def gauss(x, amp, cen, sigma):
    return amp*np.exp(-(x-cen)**2/(2.*sigma**2))

x1= np.arange(0.,100.,0.1)
x2= np.arange(0.,100.,0.09)
y1= gauss(x1, 1.,50.,5.)+ np.random.normal(size=len(x1), scale=0.1)
y2= gauss(x2, 0.8,48.4.,4.5)+ np.random.normal(size=len(x2), scale=0.1)

mod= GaussianModel()
params= mod.make_params()

mod.fit([y1,y2], params, x= [x1, x2])

我猜如果可能的话，数据必须以正确的类型传递给 mod.fit。文档只说 mod.fit 需要一个类似数组的数据输入。

我试图给它列表和数组。如果我将不同的数据集作为列表传递，我会得到一个 ValueError: setting an array element with a sequence

如果我传递一个数组，我会得到一个 AttributeError: 'numpy.ndarray' has no atribute 'exp'

所以我只是想做一些不可能的事情还是我做错了什么？

Answer 1

嗯，我认为答案是“有点”。 lmfit.Model 类旨在表示数据数组的模型。因此，如果您可以将多个数据集映射到一个 numpy ndarray（例如，使用np.concatenate），您可能可以编写一个模型函数来表示这一点，方法是为不同的数据集构建子模型并以相同的方式连接它们。

我认为任何内置模型都无法做到这一点。我还认为，一旦你开始编写复杂的模型函数，编写目标函数并不是一个很大的飞跃。也就是说，会是什么

def model_function(x, a, b, c):
   ### do some calculation with x, a, b, c values
   result = a + x*b + x*x*c
   return result

可能变成

def objective_function(params, x, data):
     vals = params.valuesdict()
     return data - model_function(x, vals['a'], vals['b'], vals['c'])

如果 do_calc() 正在做任何复杂的事情，解包参数和减去数据的额外负担非常小。而且，特别是如果某些参数将用于多个数据集，而某些参数仅用于特定数据集，则您必须在模型函数或目标函数中进行管理。在您链接到的示例中，我的答案包括数据集的循环，按名称为每个数据集挑选参数。你可能想做这样的事情。您可以通过将模型函数视为对连接数据集进行建模来在模型函数中执行此操作，但我不确定您是否真的会通过这样做获得很多。

Answer 2

我发现了问题。实际上 model.fit() 将很好地处理多个数据集的数组并执行适当的拟合。对具有多个数据集的 model.fit() 的正确调用是：

import numpy as np 
from lmfit import Model, Parameters
from lmfit.models import GaussianModel
import matplotlib.pyplot as plt

def gauss(x, amp, cen, sigma):
   "basic gaussian"
    return amp*np.exp(-(x-cen)**2/(2.*sigma**2))

x1= np.arange(0.,100.,0.1)
x2= np.arange(0.,100.,0.1)
y1= gauss(x1, 1.,50.,5.)+ np.random.normal(size=len(x1), scale=0.01)
y2= gauss(x2, 0.8,48.4,4.5)+ np.random.normal(size=len(x2), scale=0.01)

mod= GaussianModel()
params= mod.make_params()

params['amplitude'].set(1.,min=0.01,max=100.)
params['center'].set(1.,min=0.01,max=100.)
params['sigma'].set(1.,min=0.01,max=100.)

result= mod.fit(np.array([y1,y2]), params,method='basinhopping',
x=np.array([x1,x2]))

print(result.fit_report(min_correl=0.5))

fig, ax = plt.subplots()

plt.plot(x1,y1, lw=2, color='red')
plt.plot(x2,y2, lw=2, color='orange')
plt.plot(x1,result.eval(x=x1), lw=2, color='black')

plt.show()

原始代码中的问题实际上在于我的数据集没有相同的长度。但是我完全不确定如何以最优雅的方式处理这个问题？