[SciPy-user] linear (polynomial) fit with error bars

[massimo sandal]

Hi,

Fabrice Silva ha scritto:
> Using the diag input argument of leastsq :
> 
>         from scipy import optimize
>         def errfunc(a,X,Y):
>             return Y-(a[0]*X+a[1])
>         #b may be the vector containing the error bars sizes.
>         weigths = 1./b
>         a, success = optimize.leastsq(errfunc, [0,0],args=(X,Y), diag=weigths)
> 
> You here give more importance to points having small error bars.

Thanks for your advice. I am trying however to use your code, but I am
stuck upon an error.
Here is the script, that reads a very raw data file (see below):

#!/usr/bin/env python
from scipy import optimize

def errfunc(a,X,Y):
	return Y-(a[0]*X+a[1])
#b may be the vector containing the error bars sizes.

f=open("data.txt", "r")
datf=f.readlines()


numofdatapoints=len(datf)/3
xval=[0.0]*numofdatapoints
yval=[0.0]*numofdatapoints
err=[0.0]*numofdatapoints

# print len(datf)
# print datf

for i in range(numofdatapoints):
    xval[i]=float(datf[i])
    yval[i]=float(datf[i+numofdatapoints])
    err[i]=float(datf[i+2*numofdatapoints])

weigths=[0.0]*numofdatapoints
for i in range(numofdatapoints):
	weigths[i] = 1./err[i]

w=[0.0]*numofdatapoints
success=[0.0]*numofdatapoints

w, success = optimize.leastsq(errfunc, [0,0], args=(xval,yval),
diag=weigths)

print valA
print success

-------
data.txt:
118.877580092022
110.450590941286
108.684062758621
109.314800167624
103.090778781767
98.5714370869397
29.1
31.42
33.74
36.06
38.38
40.7
2.76010170015786
3.52143474842509
2.45059986418858
3.21254530326032
2.11363073382134
2.14664809861522

-------
and the script dies with the following error:

massimo at calliope:~/Python/linfit$ python linfit.py
Traceback (most recent call last):
    File "linfit.py", line 31, in <module>
      w, success = optimize.leastsq(errfunc, [0,0], args=(xval,yval),
diag=weigths)
    File "/usr/lib/python2.5/site-packages/scipy/optimize/minpack.py",
line 262, in leastsq
      m = check_func(func,x0,args,n)[0]
    File "/usr/lib/python2.5/site-packages/scipy/optimize/minpack.py",
line 12, in check_func
      res = atleast_1d(apply(thefunc,args))
    File "linfit.py", line 4, in errfunc
      return Y-(a[0]*X+a[1])
ValueError: shape mismatch: objects cannot be broadcast to a single shape

which baffles me. What should I look for to understand what I am doing
wrong?

m.
-- 
Massimo Sandal
University of Bologna
Department of Biochemistry "G.Moruzzi"

snail mail:
Via Irnerio 48, 40126 Bologna, Italy

email:
massimo.sandal at unibo.it

tel: +39-051-2094388
fax: +39-051-2094387

-------------- next part --------------
A non-text attachment was scrubbed...
Name: massimo.sandal.vcf
Type: text/x-vcard
Size: 274 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20080410/62d3445c/attachment.vcf>