Calculation of residuals: Opposite of common calculation method

[Timothy-J-Warner]

I noticed that the calculations of residuals have the opposite sign of what I expected based on the typical calculation method for residuals.

Commonly residuals are calculated:

• residual = data - model

but LMfit calculates:

• residual = model - data

Is there a specific reason the residuals are calculated this way? I know it doesn't affect squares and is an easy fix to reverse the sign but I wanted to know if there was reason for this definition. The reason I ask this is that I was confused about results I was receiving using LMfit (reversed sign for residuals) until I checked the source code and found the definition of residuals below:

def _residual(self, params, data, weights, **kwargs):
"""Return the residual.

    Default residual: ``(data-model)*weights``.

    If the model returns complex values, the residual is computed by
    treating the real and imaginary parts separately. In this case, if
    the weights provided are real, they are assumed to apply equally
    to the real and imaginary parts. If the weights are complex, the
    real part of the weights are applied to the real part of the
    residual and the imaginary part is treated correspondingly.

    Since the underlying `scipy.optimize` routines expect
    ``numpy.float`` arrays, the only complex type supported is
    ``complex``.

    The "ravels" throughout are necessary to support `pandas.Series`.

    """
    model = self.eval(params, **kwargs)
    if self.nan_policy == 'raise' and not np.all(np.isfinite(model)):
        msg = ('The model function generated NaN values and the fit '
               'aborted! Please check your model function and/or set '
               'boundaries on parameters where applicable. In cases like '
               'this, using "nan_policy=\'omit\'" will probably not work.')
        raise ValueError(msg)

    diff = model - data

    if diff.dtype is complex:
        # data/model are complex
        diff = diff.ravel().view(float)
        if weights is not None:
            if weights.dtype is complex:
                # weights are complex
                weights = weights.ravel().view(float)
            else:
                # real weights but complex data
                weights = (weights + 1j * weights).ravel().view(float)
    if weights is not None:
        diff *= weights
    return diff

[newville]

@Timothy-J-Warner I have no preference for residual = data - model or residual = model - data. Is one of these more "common"? I cannot think of any reason for choosing one sign over the other.

But, I would say that we should probably be internally consistent in the code, docstrings, and external documentation. And, as you say, lModel._residual says that it returns (data-model)*weights but it actually returns (model-data)*weights. A PR to fix that (and any other) inconsistency(ies) would be welcome.

I don't have a preference for whether code or doc is changed. Model._residual is really an implementation detail. OK, it's not __residual__, but why would anyone call Model._residual and try to add that back to the data to recover the model when it already exists in ModelResult.best_fit?

[Timothy-J-Warner]

The reason I say it is more "common" is when I do I quick search, (data - model) most commonly comes up but I didn't find a textbook reference etc. to confirm.

I can create a PR to change the documentation or code but I am unsure if it will cause issues in other locations. For example I noticed the difference in residuals using ModelResult.plot and ModelResult.plot_residuals. Taking a quick look I think ModelResult.plot_residuals also defines (model - data)

residuals = reduce_complex(self.eval()) - reduce_complex(self.data)

while ModelResult.plot appears to use ModelResult.plot_residuals rather than defining residuals again.

[newville]

@Timothy-J-Warner I think the goal should be consistency. And I think it would be fine to change the sign in the code. I have no preference.

[Timothy-J-Warner]

How do you suggest that I do a PR for this change (I would be my first pull request to a large open source project).
Should I create a fork, make the changes then initiate a pull request back to the master branch.
What tests would you suggest before initiating a pull request?

[newville]

@Timothy-J-Warner Make a fork into your own account, work there and push your work to your fork (probably doing pip install -e ".[all]" from there to get all of the developer-needed add-ons. If altering the docs, check that they build with sphinx (cd docs ; make), if altering the code, check that the tests pass (cd tests ; pytest).

When you are happy with the results, submit a Pull Request from there.

I think there are probably good GitHub tutorials.

[Timothy-J-Warner]

Thanks very much. I will create a fork and work on the consistency of the code with the documentation.