Ошибка в коэффициентах для полинома с использованием Numpy

Коэффициенты полинома предназначены для масштабированных полиномов и полиномов смещения для повышения числовой стабильности. Вы можете либо преобразовать в «нормальный» многочлен, либо использовать ряд напрямую, если вы замените off + scl*x на x, где off и scl возвращаются pol.mapparms. Чтобы преобразовать в стандартную форму (не рекомендуется), выполните pol.convert(domain=[-1, 1]).

Руипы.

В вашем примере есть три проблемы:

1. Вы подбираете полином пятого порядка только с четырьмя точками данных. Это случай недоопределенности, и он, вероятно, приведет к ранжированию. Однако это случайность, а не основная часть вашей проблемы.

2. Вы ожидаете, что pol(0) будет работать как numpy.polyval, но это не так. Я вообще-то не уверен, что он делает. Класс предоставляет __call__, который заставляет pol(0) работать, но, насколько я могу судить, документации для вызываемого объекта нет (см. Polynomial docs). numpy.polynomial.polynomial содержит собственную версию polyval. Я протестирую его np.polyval и самодельную версию test_polyval вместе.

3. Наиболее важно то, что коэффициенты порядка Polynomial класса отличаются от numpy.polyfit и numpy.polyval. В Polynomial, как вы описали, коэффициент наивысшего порядка стоит последним в списке / массиве. Однако в функциях numpy коэффициент наивысшего порядка идет первым (см. polyval docs).

В приведенном ниже фрагменте кода показано, как оценить полином, представленный вашим объектом Polynomial, при произвольном наборе значений x, а также показано, что для получения такого же поведения из numpy.polyval вы должны изменить порядок коэффициентов с помощью coef[::-1] . Я мог бы эквивалентным образом использовать numpy.fliplr, чтобы изменить порядок коэффициентов.

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!'