How do I specify error of Y-variable when fitting with lmfit?

Question

I'm almost new to Python and I'm trying to fit data from college using lmfit. The Y variable has a variable error of 3%. How do I add that error to the fitting process? I am changing from scipy's curve fit and in scipy it was really easy to do so, just creating an array with the error values and specifying the error when fitting by adding the text "sigma = [yourarray]" This is my current code:

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 are the values I measured, and eV would be the array of errors. t is the x coordinate. Thank you for your time

Answer 1

The minimize() function minimizes an array in the least-square sense, adjusting the variable parameters in order to minimize (resid**2).sum() for the resid array returned by your objective function. It really does not know anything about the uncertainties in your data or even about your data. To use the uncertainties in your fit, you need to pass in your array eV just as you pass in t and V1 and then use that in your calculation of the array to be minimized.

One typically wants to minimize Sum[ (data-model)^2/epsilon^2 ], where epsilon is the uncertainty in the data (your eV), so the residual array should be altered from data-model to (data-model)/epsilon. For your fit, you would want

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

and then use this with

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

If you use the lmfit.Model interface (designed for curve-fitting), then you could pass in weights array that multiplies data -model, and so would be 1.0 / eV to represent weighting for uncertainties (as above with minimize). Using the lmfit.Model interface and providing uncertainties would then look like this:

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()

hope that helps....