[SciPy-User] fitting sigmoidal data with weighted least squares
Daniel Mader
danielstefanmader at googlemail.com
Tue Feb 22 03:58:15 EST 2011
Dear Joseph and Jonathan,
thank you both very very much, you helped me a lot! It works as I'd
expect. However, given my data, a square weight sigma=ydata**2 makes
no sense, plain ydata or yStd is fully enough for a perfect fit in the
lower concentrations.
Once again, thank you very much and have a great day!
Daniel
2011/2/22 <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
>>
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list