与 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接口(专为曲线拟合而设计)，那么你可以传入weights乘以data -model的数组，因此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 的curve_fit来源，它会遵循 sigma(对于 1-d 情况)。

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。