[SciPy-User] fitting sigmoidal data with weighted least squares

[Jonathan Rocher]

Yes, you should give to sigma an array of errors of the same length as xs or
ys. Note that I tried you code and if you try to use 1/ys**2 without
tweaking, a "good enough" fit will not be found with the default max value
on the number of calls to the function.  One of the values of y is 0, so
that creates an infinite error bar. You might therefore tweak certain values
of your error array and potentially also play with the "maxfev" keyword as
well.

Hope this helps.
Jonathan

On Mon, Feb 21, 2011 at 9:50 AM, Daniel Mader <
danielstefanmader at googlemail.com> wrote:

> Dear Jonathan,
>
> thanks for such a quick answer! However, I am not 100% sure if I
> understood correctly:
>
> Currently, I am passing raw values for x and y to curve_fit:
> popt,pcov = scipy.optimize.curve_fit(fitfunc,xs,ys,p0=guess)
>
> Now, when I want weighting, I just need to pass an additional array to
> the function, with one entry per x,y pair?
>
> I could use the standard deviation alright but I am not sure is this
> would yield the same results as the 1/Y² weighting used by GraphPad.
>
> What would you suggest to use instead?
> weights = yin**2
> so that
> popt,pcov = scipy.optimize.curve_fit(fitfunc,xs,ys,p0=guess,sigma=weights)
>
> Thanks again in advance for any hint,
> Daniel
>
>
> 2011/2/21 Jonathan Rocher wrote:
> > Hi Daniel,
> >
> > I used this recently. You need to pass the optional argument 'sigma' to
> > curve_fit:
> > In [36]: from scipy import optimize
> >
> > In [37]: optimize.curve_fit?
> >
> > (...)
> > sigma : None or N-length sequence
> >         If not None, it represents the standard-deviation of ydata.
> >         This vector, if given, will be used as weights in the
> >         least-squares problem.
> >
> > Hope this helps,
> > Jonathan
> >
> > On Mon, Feb 21, 2011 at 8:11 AM, Daniel Mader
> > <danielstefanmader at googlemail.com> wrote:
> >>
> >> Hi,
> >>
> >> currently, I am using scipy.optimize.fit_curve in order to do least
> >> squares fitting to sigmoidal data. Basically, this works OK but now
> >> I'd like to introduce some kind of weighting to the fit.
> >>
> >> >From the help of GraphPad Prism:
> >> "Regression is most often done by minimizing the sum-of-squares of the
> >> vertical distances of the data from the line or curve. Points further
> >> from the curve contribute more to the sum-of-squares. Points close to
> >> the curve contribute little. This makes sense, when you expect
> >> experimental scatter to be the same, on average, in all parts of the
> >> curve.
> >>
> >> In many experimental situations, you expect the average distance (or
> >> rather the average absolute value of the distance) of the points from
> >> the curve to be higher when Y is higher. The points with the larger
> >> scatter will have much larger sum-of-squares and thus dominate the
> >> calculations. If you expect the relative distance (residual divided by
> >> the height of the curve) to be consistent, then you should weight by
> >> 1/Y2."
> >>
> >> This is exactly the case for my data, so I'd like to give this a try
> >> but I have no clue how.
> >>
> >> Attached is a basic script which works besides weighting. Maybe
> >> someone could point out how to pass this to the underlying
> >> scipy.optimize.leastsq function?
> >>
> >> Thanks a lot in advance,
> >> Daniel
> >>
> >> # -*- coding: utf-8 -*-
> >>
> >> import scipy, pylab, scipy.optimize
> >>
> >> def findNearest(array,value):
> >>  return abs(scipy.asarray(array)-value).argmin()
> >>
> >> def sigmoid(x,EC50,k,base,amp):
> >>  return amp / (1 + scipy.exp(-k*(x-EC50))) + base
> >>
> >> xs = [-5.80914299,
> >>      -4.60517019,
> >>      -3.5065579,
> >>      -2.30258509,
> >>      -1.2039728,
> >>      0.,
> >>      1.09861229,
> >>      2.30258509,
> >>      3.40119738
> >>      ]
> >> ys = [5.15459766e-04,
> >>      0.00000000e+00,
> >>      8.57757267e-04,
> >>      6.35666594e-03,
> >>      1.23643898e-01,
> >>      5.36029832e-01,
> >>      7.95598054e-01,
> >>      8.96318087e-01,
> >>      1.00000000e+00
> >>      ]
> >>
> >> print "guessing parameters ..."
> >> xmin = min(xs)
> >> xmax = max(xs)
> >> ymin = min(ys)
> >> ymax = max(ys)
> >> y50 = (ymax - ymin) / 2 + ymin
> >> idx = findNearest(ys,y50)
> >> EC50 = xs[idx]
> >> k = 1
> >> baseline = ymin
> >> amplitude = ymax - ymin
> >>
> >> guess = [EC50,k,baseline,amplitude]
> >> print "guess: ", guess
> >>
> >> print "fitting data ..."
> >> fitfunc = sigmoid
> >> popt,pcov = scipy.optimize.curve_fit(fitfunc,xs,ys)
> >> print "popt: ", popt
> >>
> >> x = scipy.linspace(min(xs),max(xs),100)
> >> pylab.figure()
> >> pylab.plot(xs,ys,'x', label='raw')
> >> pylab.plot(x,fitfunc(x,*guess), label='guess')
> >> pylab.plot(x,fitfunc(x,*popt), label='fit')
> >> pylab.legend()
> >> pylab.grid()
> >> pylab.show()
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>



-- 
Jonathan Rocher,
Enthought, Inc.
jrocher at enthought.com
1-512-536-1057
http://www.enthought.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110221/42975bc9/attachment.html>