linear (polynomial) fit with error bars

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

linear (polynomial) fit with error bars

massimo sandal
Hi,

Is there a scipy implementation of linear fit that keeps into account
the existence and width of error bars for each point? I googled but I
can't find that.

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

snail mail:
Via Irnerio 48, 40126 Bologna, Italy

email:
[hidden email]

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

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user

massimo.sandal.vcf (286 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: linear (polynomial) fit with error bars

Fabrice Silva-2
Le mardi 08 avril 2008 à 13:00 +0200, massimo sandal a écrit :
> Hi,
>
> Is there a scipy implementation of linear fit that keeps into account
> the existence and width of error bars for each point? I googled but I
> can't find that.

What about computing least squares approximation using weightening
coefficients inversely proportional to value incertainty (i.e. length of
each error bars) ?
--
Fabrice Silva
LMA UPR CNRS 7051

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: linear (polynomial) fit with error bars

massimo sandal
Fabrice Silva ha scritto:

> Le mardi 08 avril 2008 à 13:00 +0200, massimo sandal a écrit :
>> Hi,
>>
>> Is there a scipy implementation of linear fit that keeps into account
>> the existence and width of error bars for each point? I googled but I
>> can't find that.
>
> What about computing least squares approximation using weightening
> coefficients inversely proportional to value incertainty (i.e. length of
> each error bars) ?
Yes, that's what I was trying to ask.
But I can't find the function that does it (if it exists).

m.

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

snail mail:
Via Irnerio 48, 40126 Bologna, Italy

email:
[hidden email]

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

_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user

massimo.sandal.vcf (286 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: linear (polynomial) fit with error bars

Gerard Vermeulen-2
On Tue, 08 Apr 2008 18:13:04 +0200
massimo sandal <[hidden email]> wrote:

> Fabrice Silva ha scritto:
> > Le mardi 08 avril 2008 à 13:00 +0200, massimo sandal a écrit :
> >> Hi,
> >>
> >> Is there a scipy implementation of linear fit that keeps into
> >> account the existence and width of error bars for each point? I
> >> googled but I can't find that.
> >
> > What about computing least squares approximation using weightening
> > coefficients inversely proportional to value incertainty (i.e.
> > length of each error bars) ?
>
> Yes, that's what I was trying to ask.
> But I can't find the function that does it (if it exists).
>
> m.
>
You are looking for scipy.odr (requires scipy-0.6 or later).
It is much more useful to fit experimental data than leastsq.

Gerard
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: linear (polynomial) fit with error bars

massimo sandal
In reply to this post by 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@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:
[hidden email]

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


_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user

massimo.sandal.vcf (286 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: linear (polynomial) fit with error bars

Robert Kern-2
On Thu, Apr 10, 2008 at 7:09 AM, massimo sandal <[hidden email]> wrote:

>  massimo@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?

Print out the various individual parts of that expression to make sure
they are compatible (or use a debugger to do the same interactively).
E.g.

  print Y
  print X
  print a

In this case, the problem is that you are passing in X and Y as lists
instead of arrays. Since a[0] is a numpy scalar type that inherits
from the builtin int type, a[0]*X uses the list type's multiplication
which leaves you with an empty list.

Instead of constructing xval and yval as lists, use numpy.empty().

  xval=numpy.empty([numofdatapoints])
  yval=numpy.empty([numofdatapoints])
  err=numpy.empty([numofdatapoints])

Also, you don't need to "declare" the variables "w" and "success".

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco
_______________________________________________
SciPy-user mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/scipy-user
Loading...