Scipy - How do you do a 'weighted' least squares fit to data?

[fordie1000]

Hi,

I have performed a fit to data using scipy's 'leastsq' function. However,
I wanted to be able to weight the fit by individual errors on the data points.

Here is an example of what I'm doing at the moment :

# fitting functions for the data
fitfunc = lambda p,x: p[0] + p[1]*x
errfunc = lambda p,x,y : y-fitfunc(p,x)
 
 
x = [8365 7584 6956 6471 6149 6471 6956 7586]
y = [-0.70728186 -0.48169299  0.          0.42631948  0.9592345   0.66865681 0.24757267 -0.28228264]
errors = [ 0.10181198  0.24058398  0.10275219  0.2954378   0.09397324 0.21624879 0.10871897  0.34591412]
 
p0 = [5.0,-6.30]
out = optimize.leastsq(errfunc, p0,args=(xfit_rot_data,yfit_rot_data,yfit_err),full_output=1)
 
bestparams= out[0]
 
 

Thanks ...

[bvdet]

You assignments to x, y and errors are not valid Python lists. I do not understand what "weight the fit" is.

BV

[fordie1000]

yes ... sorry the x,y, and error lists should have commas between them ... I just copied and pasted them from a text file as an example.

What I mean by weighting the fit is that there is an error associated with
each x,y data point and when scipy does the fitting I would like it to take into
account that error so that points with huge errors don't have the same influence on the fit as a point with a very small error.

Thanks

[Glenton]

Hi

It seems that the problem is mathematical, rather than coding! I won't help with the python, cos you clearly know what you're doing, but I'm sure people would appreciate posting back your solution.

I can think of two methods, one crude and quick, the other more rigorous.

Method 1:
- Create an integer weighting, but inverting the errors (1/error), multiplying by some suitable constant, and rounding to the nearest integer.
- Create a new data set by adding multiple copies of each data point, corresponding to the above integer.
- Do a least square fit on this new data set.

Obviously by picking the constant suitably large you can get the weighting quite accurate. The big advantage is that it's a small tweak on your code.

METHOD 2:
- Create the weighted least square function yourself (Sum ((data-f(x))^2)/error).
- Use optimise to minimise this error based on your parameters (ie don't use least squares, but the other scipy optimize function).

Let me know if you don't have any joy or need some clarification and I'll make a bit more effort.

[smearedink]

In your error function, rather than simply having y-fitfunc(p,x), have (y-fitfunc(p,x))/errors, where errors is your array of 1-sigma error bars. When leastsq squares these terms and sums them, it's effectively including the 1/sigma^2 weights you want.

[eggmilk]

In your err fun:

P = 1/ERRORS
errfun = VTPV ( errfun = min )