11

I am currently using numpy.polyfit(x,y,deg) to fit a polynomial to experimental data. I would however like to fit a polynomial that uses weighting based on the errors of the points.

I have found scipy.curve_fit which makes use of weights and I suppose I could just set the function, 'f', to the form a polynomial of my desired order, and put my weights in 'sigma', which should achieve my goal.

I was wondering is there another, better way of doing this?

Many Thanks.

Jdog
  • 10,071
  • 4
  • 25
  • 42

3 Answers3

16

For weighted polynomial fitting you can use: 

numpy.polynomial.polynomial.polyfit(x, y, deg, rcond=None, full=False, w=weights)

see http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polyfit.html

Important to note that in this function the weights should not be supplied as 1/variance (which is the usual form in many weighted applications), but as 1/sigma

Although curve_fit and leastsq are much more general and powerful optimization tools than polyfit (in that they can fit just any function), polyfit has the advantage that it yields an (exact) analytical solution and is therefore probably much faster than iterative approximation methods like curve_fit and leastsq - especially in the case of fitting polynomials to multiple sets of y-data (obtained at the same x-vector)

Update: As of numpy version 1.7, numpy.polyfit also takes weights as an input (which ideally should be supplied as 1/sigma, not 1/variance)

Rolf Bartstra
  • 1,643
  • 1
  • 16
  • 19
  • Why is it `1/sigma` instead of `1/sigma**2`? (The documentation clearly and explicitly agrees with you that it is `1/sigma` and not `1/sigma**2`). – EL_DON Sep 19 '16 at 15:16
8

Take a look at http://scipy-cookbook.readthedocs.io/items/FittingData.html in particular the section 'Fitting a power-law to data with errors'. It shows how to use scipy.optimize.leastsq with a function that includes error weighting.

Tobias Kienzler
  • 25,759
  • 22
  • 127
  • 221
Bernhard
  • 8,583
  • 4
  • 41
  • 42
  • My data doesn't follow a power law however, so I cannot fit a straight line. – Jdog Jul 12 '11 at 13:22
  • 1
    That does not matter, the leastsq fitting is suitable for linear functions (linear in the paramters) which is true for polynominals. You just have to define your function to be minimized as f(x) = (data - model(x))/error(data). You weigh by a the inverse of the error on your data. Check out http://en.wikipedia.org/wiki/Least_squares#Weighted_least_squares – Bernhard Jul 12 '11 at 15:25
  • 1
    I just checked the docs of `scipy.curve_fit` it does the same as I suggest. It uses the least squares method and probably weighs the sum of the squares by 1/sigma**2. – Bernhard Jul 13 '11 at 08:38
  • Ok, just going to implement it using curve_fit. Thanks. – Jdog Jul 15 '11 at 10:35
  • Note that the method described in the link actually propagates the error into the fit values. I had a difficult time finding another example that actually does this. – Steven C. Howell Dec 30 '16 at 16:44
-2

Here is how I did it, with lots of comments!

Note: I did it with qth and nth order polynomial fits.

from numpy import *
import pylab

# get data
fn = 'cooltemp.dat'
x, y, xerr, yerr = loadtxt(fn,unpack=True, usecols=[0,1,2,3])

# create nth degree polynomial fit
n = 1
zn = polyfit(x,y,n) 
pn = poly1d(zn) # construct polynomial 

# create qth degree polynomial fit
q = 5
zq = polyfit(x,y,q) 
pq = poly1d(zq)

# plot data and fit
xx = linspace(0, max(x), 500)
pylab.plot(xx, pn(xx),'-g', xx, pq(xx),'-b')
pylab.errorbar(x, y, xerr, yerr, fmt='r.')

# customise graph
pylab.legend(['degree '+str(n),'degree '+str(q),'data'])
pylab.axis([0,max(x),0,max(y)])
pylab.xlabel('x label (unit)')
pylab.ylabel('y label (unit)')

pylab.show()
8765674
  • 1,234
  • 4
  • 17
  • 32
  • 8
    This plots a polynomial fit with error bars on the points but doesn't actually include the errors (`xerr` and `yerr`) in the construction of the polynomial fit. – Jdog Aug 17 '12 at 11:42