Iteratively calculating weights for least-squares fitting

[Felix Boschetty]

Dear lmfit-community,

I am a new user and have a question about iteratively updating weights to use during least-squares fitting.

I want to fit a line using both x and y errors. To tackle this, I would like to incorporate the x error into the weights, producing an additional 'effective' error in y.

So my weights become:

where df/dx is the partial of the fit function w.r.t. x. and sigma_x is the error estimate in x. As the partials will contain best-fit parameters, the fitting process then becomes iterative.

Initially I am attempting to do this with simple polynomials, for which lmfit appears to be overkill, but I may need non-linear or more complex fits in the future.

Is there a straightforward way to incorporate this iterative weighting into an lmfit-based workflow? Furthermore is there a simple way to save the weights for each iteration?

Finally, is it possible to calculate the partial derivatives automatically for a given fit function? i've looked at using sympy but I don't know if I can incorporate sympy functions into lmfit.

Many thanks,

Felix Boschetty

[Matt Newville]

Hi Felix,

On Sat, Dec 16, 2023 at 12:18 PM Felix Boschetty <flix.bo...@gmail.com> wrote:

Dear lmfit-community,

I am a new user and have a question about iteratively updating weights to use during least-squares fitting.

I want to fit a line using both x and y errors. To tackle this, I would like to incorporate the x error into the weights, producing an additional 'effective' error in y.

So my weights become:

where df/dx is the partial of the fit function w.r.t. x. and sigma_x is the error estimate in x. As the partials will contain best-fit parameters, the fitting process then becomes iterative.

Should we be seeing more there?  For sure, if you are asking for help with code, not including *any* code is sort of odd. 

Initially I am attempting to do this with simple polynomials, for which lmfit appears to be overkill, but I may need non-linear or more complex fits in the future. 

Is there a straightforward way to incorporate this iterative weighting into an lmfit-based workflow?

 

Well, straightforward may be in the eye of the beholder ;).  But, I am not sure I would ever describe changing the weighting at each iteration that way.

Furthermore is there a simple way to save the weights for each iteration?

Can you save them into "global data" and/or write to disk.  I am not sure I would recommend that for each function evaluation.

Finally, is it possible to calculate the partial derivatives automatically for a given fit function? i've looked at using sympy but I don't know if I can incorporate sympy functions into lmfit.

Do you mean to do `np.gradient(f)/np.gradient(x)`?    Maybe I am not understanding what "partial" you are taking.

And, just to be clear, do you mean `x` to be an array of independent data for your `y` data?   I think you probably do. 

--Matt

[Felix Boschetty]

Hi Matt,

Thanks for the quick reply.

Here's a quick example, that I hope better illustrates what I am trying to do.

from lmfit import Parameters, report_fit, minimize

import numpy as np

# Define x and y data, and associated errors

x_data = np.array([0.130638, 0.001421, 0.117021, 0.033787, 0.034766, 0.033745])

y_data = np.array([3.07e+00, 1.00e-07, 2.82e+00, 6.40e-01, 7.50e-01, 8.20e-01])

x_err = np.array([3.12668e-04, 1.62169e-04, 2.16837e-04, 9.59135e-05, 0.000119049, 7.59220e-05])

y_err = np.array([0.12, 0.10, 0.11, 0.13, 0.11, 0.12])

# Define objective function to be minimised

def obj_fun(params, x_data, y_data, x_err, y_err):

p = params.valuesdict()

a = p['a']

b = p['b']

# Automatically calculate partials from fit function, using sympy, and parameters???

partial_x = b # partial derivative wrt of x of: y = a + bx.

partial_y = 1.

var_eff = ((partial_y**2)*y_err**2 + (partial_x**2)*x_err**2) # effective variance accounting for error in x

delta_calc_obs = (a + b*x_data) - y_data

delta_var = delta_calc_obs / var_eff

return delta_var

# Initialise Parameters

params = Parameters()

params.add_many(('a', 1.), ('b', 1.))

# Perform Fit

out = minimize(fcn=obj_fun, params=params, method="lstsq", args=(x_data, y_data, x_err, y_err))

report_fit(out)

This produces a fit that to my untrained eye seems reasonable.

I'm not sure however, that the fit statistics are then being calculated using the calculated weights (1 / var_eff).

To save the intermediate weights I suspect I need to seperate the calculating the weights step into its own function that also appends them to another variabel, from my objective function.

In the long run, I would like to be able to define a fit function with a string e.g., "a + b*x" or "a + b*x + c*x**2". Then calculate the partial derivative with sympy and perform the fit as above. I have managed to produce a partial derivative that can be evaluated a multiple values of x if needed (e.g., df/dx of a + b*x + c*x**2 is b + 2*c*x):

from sympy.abc import x, y

FIT_FUNCTION = 'a + b*x + c*x**2'

PARAMS = lmfit.create_params(a=1.0, b=1.0, c=1.0)

ff_s = sympy.parsing.sympy_parser.parse_expr(FIT_FUNCTION)

ff_pd_x_s = sympy.diff(ff_s, x) # Symbolic Partial Derivative wrt x

ff_pd_x = sympy.lambdify(list(ff_pd_x_s.free_symbols), ff_pd_x_s, 'numpy') # Partial Derivative capable of accepting vectors of values.

ff_s = sympy.parsing.sympy_parser.parse_expr(FIT_FUNCTION)

But I dont know if its possible to remove the need to explicitly define the parameters within the objective function.

Currently these have to be updated manually depending on the fit equation.

Thanks,

Felix

[Felix Boschetty]

Ah, just realised that the equation I attached to my original post as a screenshot didn't make it through, hence the gap and the resulting confusion.

[Matt Newville]

On Sat, Dec 16, 2023 at 4:17 PM Felix Boschetty <flix.bo...@gmail.com> wrote:

Hi Matt,

Thanks for the quick reply.

Here's a quick example, that I hope better illustrates what I am trying to do.

from lmfit import Parameters, report_fit, minimize

import numpy as np

# Define x and y data, and associated errors

x_data = np.array([0.130638, 0.001421, 0.117021, 0.033787, 0.034766, 0.033745])

y_data = np.array([3.07e+00, 1.00e-07, 2.82e+00, 6.40e-01, 7.50e-01, 8.20e-01])

x_err = np.array([3.12668e-04, 1.62169e-04, 2.16837e-04, 9.59135e-05, 0.000119049, 7.59220e-05])

y_err = np.array([0.12, 0.10, 0.11, 0.13, 0.11, 0.12])

# Define objective function to be minimised

def obj_fun(params, x_data, y_data, x_err, y_err):

p = params.valuesdict()

a = p['a']

b = p['b']

# Automatically calculate partials from fit function, using sympy, and parameters???

partial_x = b # partial derivative wrt of x of: y = a + bx.

partial_y = 1.

var_eff = ((partial_y**2)*y_err**2 + (partial_x**2)*x_err**2) # effective variance accounting for error in x

delta_calc_obs = (a + b*x_data) - y_data

delta_var = delta_calc_obs / var_eff

return delta_var

I think one normally scales the misfit by the square-root of the variance.  Like, if x_err -> 0,  "model-data" should be scaled by y_err, not y_err**2 (just from a dimensional analysis POV).   So I think you need a square-root of your `var_eff`.

I guess I have not really thought very much about how to handle uncertainties in `x` very much - it is not something I normally have much information about.  But, I do understand that it can be important in some cases.    There is ODRPACK, and the scipy wrapping of that, but it is a pretty different interface from the `scipy.optimize` functions and we don't have a wrapping for it.

I think your approach of adding  y_err and x_err * df/dx in quadrature makes sense.  From your example, I think you might be able to this much more generally (that is, for any model) by calculating df/dx numerically:

    

    model = a + b * x_data  # or  even something more complicated

    x_err_scaled = x_err * np.gradient(model) / np.gradient(x_data)

    return (model  - y_data) / np.sqrt(y_err**2 + x_err_scaled**2)

Maybe I am missing something, but that seems like a general solution to me.  Well, it ignores correlations between uncertainties in x and y, but that could be handled if one had a covariance matrix .   

If you have a linear model, this is sort of overkill.  We're not really aiming for the most efficient way to do linear regression ;)

As a micro-optimization, you could calculate `np.gradient(x_data)` ahead of time, and not in every function evaluation.

I wonder how the results from doing this would compare to the results from ODRPACK.

# Initialise Parameters

params = Parameters()

params.add_many(('a', 1.), ('b', 1.))

# Perform Fit

out = minimize(fcn=obj_fun, params=params, method="lstsq", args=(x_data, y_data, x_err, y_err))

report_fit(out)

This produces a fit that to my untrained eye seems reasonable.

I'm not sure however, that the fit statistics are then being calculated using the calculated weights (1 / var_eff).

To save the intermediate weights I suspect I need to seperate the calculating the weights step into its own function that also appends them to another variabel, from my objective function.

In the long run, I would like to be able to define a fit function with a string e.g., "a + b*x" or "a + b*x + c*x**2". Then calculate the partial derivative with sympy and perform the fit as above. I have managed to produce a partial derivative that can be evaluated a multiple values of x if needed (e.g., df/dx of a + b*x + c*x**2 is b + 2*c*x): 

from sympy.abc import x, y

FIT_FUNCTION = 'a + b*x + c*x**2'

PARAMS = lmfit.create_params(a=1.0, b=1.0, c=1.0)

ff_s = sympy.parsing.sympy_parser.parse_expr(FIT_FUNCTION)

ff_pd_x_s = sympy.diff(ff_s, x) # Symbolic Partial Derivative wrt x

ff_pd_x = sympy.lambdify(list(ff_pd_x_s.free_symbols), ff_pd_x_s, 'numpy') # Partial Derivative capable of accepting vectors of values.

ff_s = sympy.parsing.sympy_parser.parse_expr(FIT_FUNCTION)

But I dont know if its possible to remove the need to explicitly define the parameters within the objective function.

Using sympy would limit you to models that can be expressed as a simple expression. Why not do it numerically?  

--Matt

[Felix Boschetty]

Hi Matt,

Thanks! you're very correct about both the square root, and that it is easier to calculate the partial numerically.

I'm trying to implement the EV2 method discussed in: https://pubs.acs.org/doi/10.1021/acs.analchem.0c02178. The paper compares a lot of different proposed methods for dealing with errors in x and y. Orthogonal Distance Regression (used in ODRPack) is a special case which only really works when the ratio of the variances in x and y is 1. As far as my minimal understanding gets me, proper treatment of the uncertainties in both x and y doesn't affect the parameter estimation very much, but does affect the calculated parameter uncertainties. I'm not particularly math-literate, let alone stats-literate, so take that with a pinch of salt.

Although lmfit is likely overkill for what the functions i'm trying to fit right now, there may be a need to use more complex non-linear functions in the future.

I also wanted to use lmfit because I liked the idea of writing a function to calculate EV2 that could accept a variety of different fit functions (linear or otherwise), and I hoped it would save me the overhead of writing a variety of goodness of fit functions, confidence limit calculations and potentially exploring error structure further with emcee. And on top of that the documentation is far more friendly than statsmodels :)

One more question, the way I am writing my objective function means that the model has to be defined within it, which sort of defeats the point of using lmfit to access many different models, and means I have to change the parameters I unpack depending on the model. So I need to restructure.

I'm wondering if there's a way to wrap a built in lmfit model in a function that would allow me to calculate the effective weight at each iteration.

I've seen that there is a iter_cb flag in the Minimizer options, could that be used?

Felix