# [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
```