11

Could someone explain how to get Chi^2/doF using numpy.polyfit?

casper
  • 181
  • 1
  • 2
  • 4

2 Answers2

18

Assume you have some data points

x = numpy.array([0.0, 1.0, 2.0, 3.0])
y = numpy.array([3.6, 1.3, 0.2, 0.9])

To fit a parabola to those points, use numpy.polyfit():

p = numpy.polyfit(x, y, 2)

To get the chi-squared value for this fit, evaluate the polynomial at the x values of your data points, subtract the y values, square and sum:

chi_squared = numpy.sum((numpy.polyval(p, x) - y) ** 2)

You can divide this number by the number of degrees of freedom if you like.

Sven Marnach
  • 574,206
  • 118
  • 941
  • 841
  • Thank you very much, Sven Marnach. Your answer completely solves my question. – casper Mar 30 '11 at 15:24
  • 5
    @casper: Based on your comment above, please accept this answer :) – SabreWolfy Feb 17 '12 at 09:54
  • 2
    For reference: here unitary uncertainty is assumed. The formula for the chi_square having an array s with the uncertainty on the measure is chi_squared = numpy.sum(((numpy.polyval(p, x) - y)/s) ** 2) – Daniele Jun 11 '15 at 13:15
  • without the additional weighting, isn't this just `sum of squared error` (not chi-squared)? or am i wrong? – Trevor Boyd Smith Nov 21 '19 at 17:40
  • @TrevorBoydSmith It's a χ² value assuming all "measurement errors" are the same. If you just have a set of points and look for the polynome that best approximates these points, this is what you'd generally use. If, on the other end, your data comes from some measurements, then you will also have error values associated with the data points, and you would scale the terms in the sum accordingly. – Sven Marnach Nov 22 '19 at 09:44
  • i get the whole weighting thing. my question was more about understanding `what is chi squared` vs `what is sum of squared errors`. – Trevor Boyd Smith Nov 22 '19 at 14:49
  • 1
    @TrevorBoydSmith If all errors are set to one, then they are both the same. :) – Sven Marnach Nov 22 '19 at 14:59
5

Numpy's polyfit has, at least since release 1.3, supported a full parameter. If that is set to True, polyfit will return a few more values, including the square of the residuals. Which is chi-squared (unnormalized by the degrees of freedom).

So a simple example would be

p, residuals, _, _, _ = numpy.polyfit(x, y, 2, full=True)
chisq_dof = residuals / (len(x) - 3)

I have not tried this myself with weights, but I assume polyfit does the right thing here (since numpy 1.7, polyfit accepts a parameter w to provide weights for the fit).