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

[josef.pktd at gmail.com]

On Mon, Feb 21, 2011 at 5:12 PM, Jonathan Rocher <jrocher at enthought.com> wrote:
> 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.

I'm not sure this is clear.

(If I remember correctly) if you want to weight with 1/ys, then the
argument should be sigma=ys
ys**2 would be (proportional to the) the error variance sigma**2, and
the errors are reweighted by 1/ys.

Josef

>
> 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
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>