Error in coefficients for a polynomial using Numpy

2
votes

I am using the numpy.polynomial.polynomial.Polynomial class (Numpy library) in order to fit with the method fit() certain data to a polynomial function. The polynomial obtained is all right and I can plot it and substitute points in order to get the ‘y’ value, and I get correct responses. The problem is that the .coef attribute of the Polynomial class returns a set of coefficients that are somehow normalized or changed and I cant see how. What do I mean? The code follows:

x_vid = array([0.0, 50.0, 75.0, 100.0])
y_vid = array([0.0, 30.0, 55.0, 100.0])
pol = Polynomial.fit(x_vid, y_vid, 5) # The polynomial is OK!!
print  pol.coef

The .coef attribute returns the next array:

30   38.16   17.93   9.98    2.06   1.85

The coefficients are in increasing order so, so those coefficients represent the following polynomial function:

30 + 38.16x + 17.93x^2 + 9.98x^3 + 2.06x^4 + 1.85x^5

However and here comes the problem, if I substitute any value from my range of values [0-100] there, it will not return the proper value, despite that if I do for example:

pol(0) → I will get a 0 which is OK, but it is immediate to see that in the polynomial I have written, it will not return 0 at x=0.

I think the polynomial function might be normalized or displaced. I might be facing here a mathematical problem here instead of a programming one, but any help is really welcome, because I need to write down the polynomial and I am not sure of the correct form of it. Thanks.

More info: http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.Polynomial.html#numpy.polynomial.polynomial.Polynomial

2
pythonpython-2.7numpypolynomial-mathpolynomials
Did you try using numpy.polyfit? Also, I don't get the same coefficients as you, mine are two orders of magnitude larger. - darthbith
thanks @darthbith the two methods should be very similar, anyhow the polynomial is well fitted the problem is the coefficients to write the polynomial function. I have corrected the two orders of magnitude I divided by 100 for different reasons at the time. - Ruips
I found the coefficients returned by polyfit were reasonable and the y-intercept matched the value of polyval(pol, 0). I have no idea what the coefficients returned from the Polynomial.fit() method are - darthbith

2 Answers

3
votes

The Polynomial coefficients are for scaled and offset polynomials for improved numerical stability. You can either convert to a "normal" polynomial, or use the series directly if you substitute off + scl*x for x, where off and scl are returned by pol.mapparms. To convert to standard form (not recomended), do pol.convert(domain=[-1, 1]).

1
votes

Ruips.

There are three problems with your example:

  1. You're fitting a fifth order polynomial with only four data points. This is a case of underdetermination, and it's likely to produce RankWarnings. This is incidental, though, and not the main part of your problem.

  2. You expect pol(0) to work like numpy.polyval, but it doesn't. I'm actually not sure what it does. The class provides a __call__, which makes pol(0) work, but as far as I can tell there's no documentation for the callable (see Polynomial docs). numpy.polynomial.polynomial contains its own version of polyval. I will test it, np.polyval and a homebrewed version test_polyval together.

  3. Most importantly, the Polynomial class orders coefficients differently from numpy.polyfit and numpy.polyval. In Polynomial, as you described, the highest order coefficient is last in the list/array. In the numpy functions, though, the highest order coefficient is first (see polyval docs).

The code snippet below illustrates how to evaluate the polynomial represented by your Polynomial object at an arbitrary set of x-values, and also shows that in order to get the same behavior out of numpy.polyval, you have to reverse the order of the coefficients using coef[::-1]. I could have equivalently used numpy.fliplr to reverse the coefficient order.

import numpy as np
from numpy.polynomial.polynomial import Polynomial,polyval
from numpy import array
import sys


x_vid = array([0.0, 50.0, 75.0, 100.0])
y_vid = array([0.0, 30.0, 55.0, 100.0])
pol = Polynomial.fit(x_vid, y_vid, 5) # The polynomial is OK!!

# I've written this, which should do what numpy.polynomial.polynomial.polyval 
# does, as a sanity check:
def test_polyval(polynomialInstance,xArray):
    # check that xArray is a numpy.ndarray, using ndarray.shape
    try:
        y = np.zeros(xArray.shape)
    except Exception as e:
        sys.exit('polyval error: %s'%e)

    # manually sum the polynomial terms on xArray
    for exp,c in enumerate(polynomialInstance.coef):
        y = y + c*x**exp

    return y

# Define some random x values for testing, in the range of points used
# for fitting:
x = np.random.rand(100)*100

# Compute, using our own polyval function, then Polynomial.polyval,
# and finally using numpy.polyval, making sure to reverse the
# coefficient order for the last:
y_test_polyval = test_polyval(pol,x)
y_Polynomial_polyval = polyval(x,pol.coef)
y_numpy_polyval = np.polyval(pol.coef[::-1],x)

# Make sure the two results are within machine epsilon:
if np.allclose(y_test_polyval,y_numpy_polyval) and \
        np.allclose(y_test_polyval,y_Polynomial_polyval):
    print 'Hurray!'