[SciPy-user] optimize.leastsq fit errors?

Andrew Straw strawman at astraw.com
Thu Apr 15 17:11:45 CDT 2004


I couldn't find this in scipy.  Here's what I use.  I hereby place it 
in the public domain.

import math
import Numeric as na

def leastsq(x,y, full_output=False):
     """find least-squares fit to y = a + bx

     inputs:
       x            sequence of independent variable data
       y            sequence of dependent variable data
       full_output  (optional) return dictionary of all results and stats

     outputs (not using full_output):
       a      y intercept of regression line
       b      coeffecient (slope) of regression line

     full output (using full_output):
       stats  dictionary containing statistics on fit
         key   value
         ---   -----------------------
         a     a in: y = a + bx
         b     a in: y = a + bx
         ap    a' in: x = a' + b'y
         bp    b' in: x = a' + b'y
         r2    correlation coeffecient
         var_x variance of x (sigma**2)
         var_y variance of y (sigma**2)
         cov   covariance
         SEa   standard error for a
         SEb   standard error for b
     """
     # notation from 
http://mathworld.wolfram.com/CorrelationCoefficient.html
     # and http://mathworld.wolfram.com/LeastSquaresFitting.html

     x=na.asarray(x)
     y=na.asarray(y)

     n = len(x)
     assert n == len(y)
     mean_x = na.sum(x.astype(na.Float))/n
     mean_y = na.sum(y.astype(na.Float))/n

     SSxx = na.sum( (x-mean_x)**2 )
     SSyy = na.sum( (y-mean_y)**2 )
     SSxy = na.sum( x*y ) - n*mean_x*mean_y

     # y = a + b x
     # x = ap + bp y

     b = SSxy/SSxx
     bp = SSxy/SSyy

     a = mean_y - b*mean_x

     if not full_output:
         return a, b

     ap = mean_x - bp*mean_y

     s2 = (SSyy - b*SSxy)/(n-2)
     s = math.sqrt(s2)

     SEa = s*math.sqrt( 1/n + mean_x**2/SSxx )
     SEb = s/math.sqrt(SSxx)

     stats = dict(
         r2 = b*bp,      # correlation coefficient
         var_x = SSxx/n, # variance of x (sigma**2)
         var_y = SSyy/n, # variance of y (sigma**2)
         cov = SSxy/n,   # covariance
         SEa = SEa,      # standard error for a
         SEb = SEb,      # standard error for b
         a = a,          # a in: y = a + bx
         b = b,          # b in: y = a + bx
         ap = ap,        # a' in: x = a' + b'y
         bp = bp,        # b' in: x = a' + b'y
         )
     return stats

On Apr 15, 2004, at 12:24 PM, Scott Ransom wrote:

> Hi all,
>
> Is it possible (or does anyone have any add-on functions available) to 
> get
> error estimates on the fit parameters using optimize.leastq?  (or
> alternatively, the resulting covariance matrix from the fit?)
>
> Thanks,
>
> Scott
>
> -- 
> Scott M. Ransom              Address:  McGill Univ. Physics Dept.
> Phone:  (514) 398-6492                 3600 University St., Rm 338
> email:  ransom at physics.mcgill.ca       Montreal, QC  Canada H3A 2T8
> GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 155 bytes
Desc: This is a digitally signed message part
Url : http://www.scipy.net/pipermail/scipy-user/attachments/20040415/93b40f12/PGP.bin


More information about the SciPy-user mailing list