Specifying absolute errors in curve_fit

scipy/optimize/minpack.py

@@ -448,7 +448,7 @@ def _weighted_general_function(params, xdata, ydata, function, weights):
     return weights * (function(xdata, *params) - ydata)
 
 
-def curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw):
+def curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, **kw):
     """
     Use non-linear least squares to fit a function, f, to data.
 
@@ -470,18 +470,32 @@ def curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw):
         values will all be 1 (if the number of parameters for the function
         can be determined using introspection, otherwise a ValueError
         is raised).
-    sigma : None or N-length sequence
-        If not None, this vector will be used as relative weights in the
+    sigma : None or N-length sequence, optional
+        If not None, these values are used as weights in the
         least-squares problem.
+    absolute_sigma : bool, optional
+        If False, `sigma` denotes relative weights of the data points.
+        The returned covariance matrix `pcov` is based on *estimated*
+        errors in the data, and is not affected by the overall
+        magnitude of the values in `sigma`. Only the relative
+        magnitudes of the `sigma` values matter.
+
+        If True, `sigma` describes one standard deviation errors of
+        the input data points. The estimated covariance in `pcov` is
+        based on these values.
 
     Returns
     -------
     popt : array
         Optimal values for the parameters so that the sum of the squared error
         of ``f(xdata, *popt) - ydata`` is minimized
     pcov : 2d array
-        The estimated covariance of popt.  The diagonals provide the variance
-        of the parameter estimate.
+        The estimated covariance of popt. The diagonals provide the variance
+        of the parameter estimate. To compute one standard deviation errors
+        on the parameters use ``perr = np.sqrt(np.diag(pcov))``.
+
+        How the `sigma` parameter affects the estimated covariance
+        depends on `absolute_sigma` argument, as described above.
 
     See Also
     --------
@@ -497,13 +511,13 @@ def curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw):
     >>> import numpy as np
     >>> from scipy.optimize import curve_fit
     >>> def func(x, a, b, c):
-    ...     return a*np.exp(-b*x) + c
+    ...     return a * np.exp(-b * x) + c
 
-    >>> x = np.linspace(0,4,50)
-    >>> y = func(x, 2.5, 1.3, 0.5)
-    >>> yn = y + 0.2*np.random.normal(size=len(x))
+    >>> xdata = np.linspace(0, 4, 50)
+    >>> y = func(xdata, 2.5, 1.3, 0.5)
+    >>> ydata = y + 0.2 * np.random.normal(size=len(xdata))
 
-    >>> popt, pcov = curve_fit(func, x, yn)
+    >>> popt, pcov = curve_fit(func, xdata, ydata)
 
     """
     if p0 is None:
@@ -537,11 +551,16 @@ def curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw):
         msg = "Optimal parameters not found: " + errmsg
         raise RuntimeError(msg)
 
-    if (len(ydata) > len(p0)) and pcov is not None:
-        s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
-        pcov = pcov * s_sq
-    else:
-        pcov = inf
+    if pcov is None:
+        # indeterminate covariance
+        pcov = zeros((len(popt), len(popt)), dtype=float)
+        pcov.fill(inf)
+    elif not absolute_sigma:
+        if len(ydata) > len(p0):
+            s_sq = (asarray(func(popt, *args))**2).sum() / (len(ydata) - len(p0))
+            pcov = pcov * s_sq
+        else:
+            pcov.fill(inf)
 
     if return_full:
         return popt, pcov, infodict, errmsg, ier

scipy/optimize/tests/test_minpack.py

@@ -336,6 +336,45 @@ def f_double_gauss(x, x0, x1, A0, A1, sigma, c):
         popt, pcov = curve_fit(f_double_gauss, x, y, guess, maxfev=10000)
         assert_allclose(popt, good, rtol=1e-5)
 
+    def test_pcov(self):
+        xdata = np.array([0, 1, 2, 3, 4, 5])
+        ydata = np.array([1, 1, 5, 7, 8, 12])
+        sigma = np.array([1, 2, 1, 2, 1, 2])
+
+        def f(x, a, b):
+            return a*x + b
+
+        popt, pcov = curve_fit(f, xdata, ydata, p0=[2, 0], sigma=sigma)
+        perr_scaled = np.sqrt(np.diag(pcov))
+        assert_allclose(perr_scaled, [ 0.20659803, 0.57204404], rtol=1e-3)
+
+        popt, pcov = curve_fit(f, xdata, ydata, p0=[2, 0], sigma=3*sigma)
+        perr_scaled = np.sqrt(np.diag(pcov))
+        assert_allclose(perr_scaled, [ 0.20659803, 0.57204404], rtol=1e-3)
+
+        popt, pcov = curve_fit(f, xdata, ydata, p0=[2, 0], sigma=sigma,
+                               absolute_sigma=True)
+        perr = np.sqrt(np.diag(pcov))
+        assert_allclose(perr, [0.30714756, 0.85045308], rtol=1e-3)
+
+        popt, pcov = curve_fit(f, xdata, ydata, p0=[2, 0], sigma=3*sigma,
+                               absolute_sigma=True)
+        perr = np.sqrt(np.diag(pcov))
+        assert_allclose(perr, [3*0.30714756, 3*0.85045308], rtol=1e-3)
+
+        # infinite variances
+
+        def f_flat(x, a, b):
+            return a*x
+
+        popt, pcov = curve_fit(f_flat, xdata, ydata, p0=[2, 0], sigma=sigma)
+        assert_(pcov.shape == (2, 2))
+        assert_array_equal(pcov, np.inf)
+
+        popt, pcov = curve_fit(f, xdata[:2], ydata[:2], p0=[2, 0])
+        assert_(pcov.shape == (2, 2))
+        assert_array_equal(pcov, np.inf)
+
 
 class TestFixedPoint(TestCase):