使用 lmfit 拟合时如何指定 Y 变量的误差？

Question

我几乎是 Python 新手，我正在尝试使用 lmfit 拟合大学的数据。 Y 变量的变量误差为 3%。如何将该错误添加到拟合过程中？我正在改变 scipy 的曲线拟合，而在 scipy 中这样做真的很容易，只需创建一个包含错误值的数组并在拟合时通过添加文本“sigma = [yourarray]”来指定错误 这是我当前的代码：

from lmfit import Minimizer, Parameters, report_fit
import matplotlib.pyplot as plt
w1, V1, phi1, scal1 = np.loadtxt("./FiltroPasaBajo_1.txt", delimiter = "\t", unpack = True)

t = w1
eV= V1*0.03 + 0.01

def funcion(parametros, x, y):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return modelo - y
parametros = Parameters()
parametros.add('R', value = 1000, min= 900, max = 1100)
parametros.add('C', value = 1E-6, min = 1E-7, max = 1E-5)

fit = Minimizer(funcion, parametros, fcn_args=(t,V1))
resultado = fit.minimize()


final = V1 + resultado.residual
report_fit(resultado)

try:
    plt.plot(t, V1, 'k+')
    plt.plot(t, final, 'r')
    plt.show()
except ImportError:
    pass

V1 是我测量的值，eV 是错误数组。 t 是 x 坐标。 感谢您的宝贵时间

Answer 1

minimize() 函数在最小二乘意义上最小化数组，调整变量参数以最小化目标函数返回的resid 数组的(resid**2).sum()。它真的对您的数据中的不确定性甚至您的数据一无所知。要使用适合的不确定性，您需要传入数组eV，就像传入t 和V1 一样，然后在计算要最小化的数组时使用它。

人们通常希望最小化Sum[ (data-model)^2/epsilon^2 ]，其中epsilon 是数据中的不确定性（您的eV），因此应将残差数组从data-model 更改为(data-model)/epsilon。为了你的健康，你会想要

def funcion(parametros, x, y, eps):
    R = parametros['R'].value
    C = parametros['C'].value

    modelo = 4/((1+(x**2)*(R**2)*(C**2))**1/2)
    return (modelo - y)/eps

然后将其与

一起使用

fit = Minimizer(funcion, parametros, fcn_args=(t, V1, eV))
resultado = fit.minimize()
...

如果您使用lmfit.Model 接口（专为曲线拟合而设计），那么您可以传入与data -model 相乘的weights 数组，1.0 / eV 也可以用来表示不确定性的权重（如上minimize)。使用lmfit.Model 接口并提供不确定性将如下所示：

from lmfit import Model
# model function, to model the data
def func(t, r, c):
    return 4/((1+(t**2)*(r**2)*(c**2))**1/2)

model  = Model(func)
parametros = model.make_params(r=1000, c=1.e-6)

parametros['r'].set(min=900, max=1100)
parametros['c'].set(min=1.e-7, max=1.e-5)

resultado = model.fit(V1, parametros, t=t, weights=1.0/eV)

print(resultado.fit_report())

plt.errorbar(t, V1, eV, 'k+', label='data')
plt.plot(t, resultado.best_fit, 'r', label='fit')
plt.legend()
plt.show()

希望对您有所帮助....

Answer 2

我认为您不能直接在 fit.minimize() 中提供 sigma。

但是我看到 fit.minimize() 使用 scipy 的 minimumsq 方法（默认情况下），这与 scipy 的 curve_fit 使用的方法相同。

如果您查看scipy's curve_fit source，它会跟随 sigma（对于一维情况）。

 transform = 1.0 / sigma
 jac = _wrap_jac(jac, xdata, transform)
 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)

由于 fit.minimize() 允许您为 minimumsq 传递 kwargs (Dfun)，因此您可以按照 scipy curve_fit 中的方式传递 jac。