# [Scipy-svn] r5548 - trunk/scipy/optimize

scipy-svn@scip... scipy-svn@scip...
Thu Feb 12 22:46:34 CST 2009

Author: oliphant
Date: 2009-02-12 22:46:29 -0600 (Thu, 12 Feb 2009)
New Revision: 5548

Modified:
trunk/scipy/optimize/minpack.py
Log:
Add the residual factor to convert the leastsq 'cov' matrix into a true estimate of the covariance matrix.

Modified: trunk/scipy/optimize/minpack.py
===================================================================
--- trunk/scipy/optimize/minpack.py	2009-02-13 04:15:48 UTC (rev 5547)
+++ trunk/scipy/optimize/minpack.py	2009-02-13 04:46:29 UTC (rev 5548)
@@ -191,9 +191,12 @@
unsuccessful call.

cov_x -- uses the fjac and ipvt optional outputs to construct an
-             estimate of the covariance matrix of the solution.
+             estimate of the jacobian around the solution.
None if a singular matrix encountered (indicates
-             infinite covariance in some direction).
+             very flat curvature in some direction).  This
+             matrix must be multiplied by the residual standard
+             deviation to get the covariance of the parameter
+             estimates --- see curve_fit.
infodict -- a dictionary of optional outputs with the keys:
'nfev' : the number of function calls
'fvec' : the function evaluated at the output
@@ -259,6 +262,8 @@

fixed_point -- scalar and vector fixed-point finder

+      curve_fit -- find parameters for a 1-d curve-fitting problem.
+
"""
x0 = array(x0,ndmin=1)
n = len(x0)
@@ -323,7 +328,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):
+def curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw):
"""Use non-linear least squares to fit a function, f, to data.

Parameters
@@ -350,13 +355,14 @@
Optimal values for the parameters so that the sum of the squared error of
f(xdata, *popt) - ydata is minimized
pcov : 2d array
-        A covariance matrix shouwing the curvature of the sum-of-squares
-        residual near the returned solution.  Returned directly from the call
-        to scipy.optimize.leastsq.
+        A covariance matrix providing an estimate of the covariance matrix of
+        the parameter estimates.   This may not be accurate when sigma
+        is given.

Notes
-----
The algorithm uses the Levenburg-Marquardt algorithm: scipy.optimize.leastsq
+    Additional keyword arguments are passed directly to that algorithm.

Example
-------
@@ -388,10 +394,18 @@
else:
func = _weighted_general_function
args += (1.0/asarray(sigma),)
-    popt, pcov, infodict, mesg, ier = leastsq(func, p0, args=args, full_output=1)
+    popt, pcov, infodict, mesg, ier = leastsq(func, p0, args=args, full_output=1, **kw)

if ier != 1:
+
+    if (len(xdata) > len(p0)):
+        s_sq = (func(popt, *args)**2).sum()/(len(xdata)-len(p0))
+        if sigma is not None:
+            s_sq /= args[-1].sum()
+        pcov = pcov * s_sq
+    else:
+        pcov = None

return popt, pcov