Bugfix: propagation of measurement uncertainties into parameter covariances

astropy/modeling/fitting.py

@@ -1240,15 +1240,26 @@
                     for _ in weights * np.array(model.fit_deriv(x, y, *params))
                 ]
 
-    def _compute_param_cov(self, model, y, init_values, cov_x, fitparams, farg):
+    def _compute_param_cov(
+        self, model, y, init_values, cov_x, fitparams, farg, weights=None
+    ):
         # now try to compute the true covariance matrix
         if (len(y) > len(init_values)) and cov_x is not None:
-            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
+            self.fit_info["param_cov"] = cov_x
+            if weights is None:
+                # if there are no measurement uncertainties given in `weights`,
+                # fall back on the default behavior in scipy.optimize.curve_fit
+                # when `absolute_sigma == False`. If there are uncertainties,
+                # assume they are "absolute" and not "relative".
+                # For details, see curve_fit:
+                #   https://github.com/scipy/scipy/blob/
+                #   c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/
+                #   optimize/_minpack_py.py#L591-L602
+                sum_sqrs = np.sum(self.objective_function(fitparams, *farg) ** 2)
+                dof = len(y) - len(init_values)
+                self.fit_info["param_cov"] *= sum_sqrs / dof
         else:
             self.fit_info["param_cov"] = None
-
         if self._calc_uncertainties is True:
             if self.fit_info["param_cov"] is not None:
                 self._add_fitting_uncertainties(model, self.fit_info["param_cov"])
@@ -1306,9 +1317,14 @@
         z : array, optional
            input coordinates
         weights : array, optional
-            Weights for fitting.
-            For data with Gaussian uncertainties, the weights should be
-            1/sigma.
+            Weights for fitting. For data with Gaussian uncertainties, the weights
+            should be 1/sigma.
+
+            .. versionchanged:: 5.3
+                Calculate parameter covariances while accounting for ``weights``
+                as "absolute" inverse uncertainties. To recover the old behavior,
+                choose ``weights=None``.
+
         maxiter : int
             maximum number of iterations
         acc : float
@@ -1349,7 +1365,9 @@
             model_copy, farg, maxiter, acc, epsilon, estimate_jacobian
         )
 
-        self._compute_param_cov(model_copy, y, init_values, cov_x, fitparams, farg)
+        self._compute_param_cov(
+            model_copy, y, init_values, cov_x, fitparams, farg, weights
+        )
 
         model.sync_constraints = True
         return model_copy

astropy/modeling/tests/test_fitters.py

@@ -619,7 +619,7 @@ def test_param_cov(self, fitter):
 
         with NumpyRNGContext(_RANDOM_SEED):
             x = np.linspace(0, 1, 100)
-            # y scatter is amplitude ~1 to make sure covarience is
+            # y scatter is amplitude ~1 to make sure covariance is
             # non-negligible
             y = x * a + b + np.random.randn(len(x))
 
@@ -638,6 +638,40 @@ def test_param_cov(self, fitter):
         assert_allclose(fmod.parameters, beta.ravel())
         assert_allclose(olscov, fitter.fit_info["param_cov"])
 
+    @pytest.mark.parametrize("fitter", non_linear_fitters)
+    def test_param_cov_with_uncertainties(self, fitter):
+        """
+        Tests that the 'param_cov' fit_info entry gets the right answer for
+        *linear* least squares, where the answer is exact
+        """
+        fitter = fitter()
+
+        a = 2
+        b = 100
+
+        with NumpyRNGContext(_RANDOM_SEED):
+            x = np.linspace(0, 1, 100)
+            # y scatter is amplitude ~1 to make sure covariance is
+            # non-negligible
+            y = x * a + b + np.random.normal(size=len(x))
+            sigma = np.random.normal(loc=1, scale=0.1, size=len(x))
+
+        # compute the ordinary least squares covariance matrix
+        # accounting for measurement uncertainties `sigma`
+        X = np.vstack([x, np.ones(len(x))]).T
+        inv_N = np.linalg.inv(np.diag(sigma) ** 2)
+        cov = np.linalg.inv(X.T @ inv_N @ X)
+        beta = cov @ X.T @ inv_N @ y.T
+
+        # now do the non-linear least squares fit
+        mod = models.Linear1D(a, b)
+
+        with pytest.warns(AstropyUserWarning, match=r"Model is linear in parameters"):
+            fmod = fitter(mod, x, y, weights=sigma**-1)
+
+        assert_allclose(fmod.parameters, beta.ravel())
+        assert_allclose(cov, fitter.fit_info["param_cov"])
+
 
 class TestEntryPoint:
     """Tests population of fitting with entry point fitters"""

docs/changes/modeling/14519.api.rst

@@ -0,0 +1,2 @@
+Propagate measurement uncertainties via the ``weights`` keyword argument into the
+parameter covariances.

docs/whatsnew/5.3.rst

@@ -19,6 +19,7 @@ In particular, this release includes:
 .. * :ref:`whatsnew-5.3-unit-formats-fraction`
 .. * :ref:`whatsnew-5.3-unit-formats-order`
 .. * :ref:`whatsnew-5.3-nddata-collapse-arbitrary-axes`
+.. * :ref:`whatsnew-5.3-modeling-measurement-uncertainties`
 .. * :ref:`whatsnew-5.3-coordinates-refresh-site-registry`
 
 In addition to these major changes, Astropy v5.3 includes a large number of
@@ -198,6 +199,7 @@ we can take the sum of ND masked quantities along the ``1`` axis like so::
      array([ True, False]),
      StdDevUncertainty([1.41421356, 1.73205081]))
 
+
 .. _whatsnew-5.3-coordinates-refresh-site-registry:
 
 Refresh cached observatory site registry for |EarthLocation| methods
@@ -214,6 +216,16 @@ however, is updated from time to time to include new locations.  The user
 may refresh their locally cached site registry by passing the new
 ``refresh_cache=True`` option to these two functions.
 
+.. _whatsnew-5.3-modeling-measurement-uncertainties:
+
+Support for collapse operations on arbitrary axes in ``nddata``
+===============================================================
+
+Propagate measurement uncertainties into the best-fit parameter covariances
+via the ``weights`` keyword argument in non-linear fitters. Decreasing
+the ``weights`` will now increase the uncertainties on the best-fit parameters.
+
+
 Full change log
 ===============