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

[Daniel Mader]

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()