ADR Suggestion Fitting weights

[damskii9992]

General

When fitting data to a model it is possible to add weights to different data-points to prioritize the fit to certain parts of the data which might be more important or more accurate than other parts. This is especially useful when working with data with uncertainties, as it is possible to weigh more accurate data-points higher than more inaccurate points.

Current Implementation

Currently, if the user does not provide a list of weights to the fitter, a list of default weights will be used, according to:

if weights is None:
    weights = np.sqrt(np.abs(y))

i.e. the weights are set to the square-root of the data values y. This is an issue for 2 reasons:

1. When the data y-values are 0, the weights become 0 and since we divide by weights we get Inf/ NaN errors.
2. If the data y-values varies greatly, such as in a diffraction spectrum, the peaks would be valued infinitely higher than any other points in the spectrum. I.e. this is in general an issue when the data y-values are not neutron counts.

Proposed implementation

Based on the discussion of this issue in #133, we propose to completely remove default weights from the core EasyScience library and require users to provide a list of weights explicitly whenever making a fit. We additionally propose to add a check on this input to ensure that the user doesn't accidentally add weights with 0-values or NaNs/Infs to increase robustness.

Sub-projects can then implement their own defaults depending on the technique of the relevant project.

[seventil]

This needs to be converted to ADR

[henrikjacobsenfys]

I would furthermore like to discuss a better naming scheme than "weights", or at least making it very clear what we mean by it. The bug mentioned in #204 is caused by me assuming "weight" meant the same in lmfit and Bumps.

I am personally a fan of passing either sigma, (could also be called dy or y_error) or variance.

[henrikjacobsenfys]

Bumps and LMfit right now don't agree. This is because Bumps assumes that the weights it receives (called dy) are uncertainties (sigma) on the points. LMFit, on the other hand, has two settings: scale_covar=False, which is the same as Bumps, or scale_covar=True, (which is the default and therefore what we do), which assumes that the weights are not uncertainties but just some arbitrary scaling. In such cases, the result then gets scaled by sqrt(chi_red^2). This can be correct if chi_red^2 is larger than 1, since that implies that the weights are not sigma (or that the model is wrong).

We should fix this so LMFit gives the same result as Bumps. We should also make it more clear what we are doing. I therefore again propose that instead of asking for x,y,weights, we ask for x,y,sigma or x,y,dy

Here's a script that illustrates the issue (Generated by Claude, comments by me)

import numpy as np
from easyscience import AvailableMinimizers, Fitter
from easyscience.base_classes import ModelBase
from easyscience.variable import Parameter

# Make a simple model to fit data
class Line(ModelBase):
    def __init__(self):
        super().__init__()
        self._m = Parameter(name='m', value=1.0)
        self._c = Parameter(name='c', value=0.0)

    @property
    def m(self): return self._m
    @property
    def c(self): return self._c

    def __call__(self, x):
        return self._m.value * x + self._c.value

# Generate random data with a seed for reproducibility
rng = np.random.default_rng(1)
x = np.linspace(0, 10, 50)
sigma = np.ones_like(x)
y = 2.5 * x + 1.5 + rng.normal(0, 2.0, x.size)   # random noise is larger than sigma to make chi2_red larger than 1

# Fit using LMFit and Bumps and print the results
for minimizer in (AvailableMinimizers.LMFit, AvailableMinimizers.Bumps):
    line = Line()
    fitter = Fitter(line, line)
    fitter.switch_minimizer(minimizer)
    result = fitter.fit(x=x, y=y, weights=1 / sigma)
    print(f'{minimizer.name:6s} m = {line.m.value:.4f} +- {line.m.error:.4f}   '
          f'c = {line.c.value:.4f} +- {line.c.error:.4f}   reduced chi2 = {result.reduced_chi2:.3f}')

# Calculate the result analytically. 
# The derivate of a straight line mx+c with respect to m is x, and with respect to c is 1. We divide this by sigma to get the Jacobian of the weighted residuals

J = np.column_stack([x, np.ones_like(x)]) / sigma[:, None]

# For a straight line, the covariance is the inverse of $J^T J$, and the error is the square root of that

m_err, c_err = np.sqrt(np.diag(np.linalg.inv(J.T @ J))) 
print(f'exact  m error = {m_err:.4f}   c error = {c_err:.4f}')

Output:

LMFit  m = 2.4513 +- 0.0861   c = 1.6711 +- 0.4995   reduced chi2 = 3.214
Bumps  m = 2.4513 +- 0.0480   c = 1.6711 +- 0.2787   reduced chi2 = 3.214
exact  m error = 0.0480   c error = 0.2787

The analytical result matches Bumps but not LMFit.

We should also consider if we want to add chi_red^2 rescaling as an option for users if their chi_red^2 is large.

[damskii9992]

I'd like to propose to call it uncertaintes rather than sigma or dy. Just to be clear about what it is. Or variances, if that is what we are using.

I also think we should have both a weight argument, and this uncertainty argument. Because we definitely should still be able to scale our fit to prioritize certain parts more than others, regardless of their uncertainties. This weight argument can default to 1, and might be rarely used, but its good to have. By having it we also ensure less confusion, as users will see both a weight argument, and an uncertainty argument, thereby knowing that they are different and represent different things.

[henrikjacobsenfys]

Well, sigma or perhaps sigma_y is more descriptive. It tells users exactly what it is.

How would sigma and weights interact if a user supplies both?

[damskii9992]

Well, sigma or perhaps sigma_y is more descriptive. It tells users exactly what it is.

It tells a physicist exactly what it is. For anyone else it is confusing. Not everyone knows that uncertainties are commonly symbolized with a sigma 😅