# Solving a polynomial with NumPy

I'm trying to do something that I thought would be very straightforward but somehow I'm struggling.

I have a time series and I want to extrapolate it, assuming a linear trend, to forecast when will it reach a certain point.

For that, I tried using NumPy to fit the polynomial

$$y = c_0 + c_1 x$$

with this code:

from numpy.polynomial import Polynomial

p = Polynomial.fit(arr[:, 0], arr[:, 1], 1)

I can evaluate it forward by doing p(x_val).

The last step is to solve $$x$$ for a certain $$\hat{y}$$, therefore

$$x = \frac{\hat{y} - c_0}{c_1}$$

but it turns out I cannot use p.coef directly, because there's some sort of smart mapping happening. Indeed, naïvely multiplying the coefficients no longer suffices:

In [85]: p(0)
Out[85]: 10.173221334851732

In [86]: p.coef[0] + p.coef[1] * 0
Out[86]: 93.18158032460339

In [87]: p(1e6)
Out[87]: 33.488650125538086

In [88]: p.coef[0] + p.coef[1] * 1e6
Out[88]: 83008452.17133197

What is the simplest way of solving for $$x$$ given the NumPy polynomial then?

• It seems to be working for me. Check this out. Jan 23 at 14:41
• Thanks @nicoguaro! So the key is the .convert() method then. Did you get that from the docs? Anyway, do you want to expand your comment to an answer? And I'll mark it
– user782
Jan 23 at 14:48
• As mentioned in the docs, as well as this answer on StackOverflow, the coefficients stored are scaled and shifted. Jan 23 at 14:49
• > Note that the coefficients are given in the scaled domain defined by the linear mapping between the window and domain. convert can be used to get the coefficients in the unscaled data domain. numpy.org/doc/stable/reference/… I had to squint a bit to see it, but ok. Thanks!
– user782
Jan 23 at 14:53

It seems to work correctly.

Keep in mind that the domain is scaled and shifted so you need to use the .convert() method to get the coefficients in the original domain.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng()
x = np.arange(10)
y = np.arange(10) + 0.01*rng.standard_normal(10)

p_fitted = np.polynomial.Polynomial.fit(x, y, deg=1)

print(p_fitted.convert().coef)

This snippet gives the result

[-0.00559272  1.00221273]