Bugfix: propagation of measurement uncertainties into parameter covariances

[bmorris3]

Description

astropy.modeling can fit nonlinear models to observations with uncertainties, and return best-fit parameters and their covariances. However, if you increase the measurement uncertainties and re-run the fit, you will find that the parameter covariances do not respond to the changing uncertainties. An example is collapsed below.

Example

Example

The code below currently runs without error on main, meaning the parameter covariance does not respond to increases in the uncertainties in the observations:

import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling import fitting
from astropy.modeling.models import Linear1D

linear_model = Linear1D()
fitter = fitting.LMLSQFitter(calc_uncertainties=True)

err_scales = np.geomspace(1, 100, 10)

np.random.seed(42)
x = np.linspace(0, 1)
yerrs = np.random.normal(
    loc=0.1, 
    scale=0.01, 
    size=x.size
)
true_slope, true_intercept = 2, 3

y = np.random.normal(
    loc=true_slope * x + true_intercept, 
    scale=yerrs,
    size=x.size
)

uncertainties = []

for err_scale in err_scales:
    fit_model = fitter(
        linear_model, x, y, weights=1/err_scale/yerrs
    )

    uncertainties.append(fit_model.intercept.std)
    
# this assert block should FAIL because we don't expect the
# best-fit parameter uncertainties to be identical as we 
# increase the measurement uncertainties:
for i in range(1, len(uncertainties)):
    np.testing.assert_almost_equal(uncertainties[0], uncertainties[i])

If the code were working correctly, the parameter uncertainties should increase.

---

This problem arises because astropy.modeling.fitting renormalizes the covariance matrix on these lines:

astropy/astropy/modeling/fitting.py

Lines 1246 to 1248 in eccdccf

	sum_sqrs = np.sum(self.objective_function(fitparams, *farg) ** 2)
	dof = len(y) - len(init_values)
	self.fit_info["param_cov"] = cov_x * sum_sqrs / dof

An analogous normalization is done in scipy's curve_fit but only in the case where absolute_sigma == False. We can see from the definition of this argument in the curve_fit docstring that this choice is equivalent to assuming that only the relative scale of the uncertainties matters, not their actual ("absolute") values.

This is not how I (or most astronomers?) might expect the weight = 1 / sigma parameter in astropy.modeling to behave. Usually when we bother to specify uncertainties, we are giving "absolute" ones, not "relative" ones.

This PR preserves the current behavior if no measurement uncertainties are provided via the weights kwarg. If weights is not None, the renormalization for "relative" sigmas is skipped, and the measurement uncertainties affect the covariance matrix.

I've added a test with an ordinary least squares example with heteroscedastic errors, and check that the resulting covariances are correct.

[github-actions]

Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see "When to rebase and squash commits".
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• Is a milestone set? Milestone must be set but we cannot check for it on Actions; do not let the green checkmark fool you.
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

[WilliamJamieson]

Should we also give the option to toggle the behavior like the absolute_sigma option that scipy.optimize.curve_fit gives users for the case when no weights are passed?

[bmorris3]

I understand giving people options is generally a good thing, and the generic name of the kwarg is weights. But we say in the docstring that weights = 1 / sigma, which is not true if absolute_sigma == False.

[WilliamJamieson]

Right, so then the question becomes "should we add an option to disable this change". That is add an option which allows users to return the behavior of the weights to how they previously worked? My concern here lies with the fact that this change will technically change this computation's output in some cases. This can have knock-on effects for downstream users, who may wish to return to the previous behavior.

@pllim what is your opinion here?

[pllim]

That is add an option which allows users to return the behavior of the weights to how they previously worked?

Unless the old way is completely wrong, it is generally good to provide backward compatibility.

[bmorris3]

I would argue the old way is completely wrong. We offered a keyword to adjust the weights, but did not compute the result based on the weights.

[bmorris3]

Also – the user can have "perfect" backwards compatibility if they set weights = None.

[pllim]

As I told Brett offline, API change is acceptable as long as it is not backported. And if an extra option lets people get the old results, then fine by me but ultimately it is up to @astropy/modeling .

[WilliamJamieson]

In principle, I am okay with this change so long as it is not backported. However, we should probably advertise the change with a "whats-new". So that it is clear when reading the changes for astropy 5.3 that this behavior will be slightly different.

[WilliamJamieson]

These changes look fine to me. Please add a "whats new" item detailing both how the results will change and the circumstance that the results will change under. I think this is necessary since this is a very subtle API change could be difficult to figure out.

[bmorris3]

Thanks @WilliamJamieson!

[saimn]

Too late to fix, but the whatsnew section title is wrong. (Please be more careful in future PRs)

[bmorris3]

Yikes! My apologies 😱

[pllim]

Hmm, this is super confusing... Is there a way to fix just the doc without doing another release?

https://docs.astropy.org/en/stable/whatsnew/5.3.html#whatsnew-5-3-modeling-measurement-uncertainties

Also, while we're at it, is it normal to have dev tag displayed on the "stable" URL?

[saimn]

It will be for the next bugfix release, 5.3 is done, but good point, we need to fix that in the v5.3.x branch.

About the dev tag, this is the goal of #14860. (Not a new issue, just managed to find a solution before forgetting about it ;)).

[pllim]

I fixed it directly on v5.3.x 31f3290

[saimn]

Thanks !