I'd like to vectorize following function
def nppolyfit(pnp_array, **kwargs):
""" Moving polyfit
"""
win_size = kwargs['length']
degree = kwargs['degree']
xdata = range(win_size)
res = np.zeros(pnp_array.shape)
for i in range(win_size, len(pnp_array) + 1):
res[i-1] = np.poly1d(np.polyfit(xdata , pnp_array[i - win_size : i], deg = degree))(win_size)
return res
What was done so far:
def rolling_window(a, window):
shp = a.shape[:-1] + (a.shape[-1] - window + 1, window)
strides = a.strides + (a.strides[-1],)
return np.lib.stride_tricks.as_strided(a, shape=shp, strides=strides)
def nppolyfitv(pnp_array, **kwargs):
""" Moving polyfit
"""
win_size = kwargs['length']
degree = kwargs['degree']
xdata = np.arange(win_size)
ydata = rolling_window(pnp_array, win_size).T
fit = np.polyfit(xdata, ydata, deg = degree)
res = np.zeros(len(pnp_array))
res[win_size-1:] = np.polynomial.polynomial.polyval(np.zeros(len(pnp_array)), fit).T[len(pnp_array) - 1,:]
return res
But it looks like I am missing something or doing wrong. Could you please correct me? Maybe there is another more effective solution? Thanks.
Test case:
import numpy as np
npd = np.arange(30)
win_size1 = 11
degree = 1
c1 = nppolyfit(npd, length=win_size1, degree=degree)
c1v = nppolyfitv(npd, length=win_size1, degree=degree)
print(c1)
print(c1v)
And results are:
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 11. 12. 13. 14. 15.
16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 30.
59. 88. 117. 146. 175. 204. 233. 262. 291. 320. 349. 378.
407. 436. 465. 494. 523. 552.]
The polyfit
method returns polynomial coefficients, highest powers first.
The polyval
method expects the coefficients to have lowest powers first. Take this into account when feeding the output of one method to another.
Also, the x-argument of polyval
is illogical: np.zeros(len(pnp_array))
. Why do you ask polyval
to evaluate the polynomial at the same point 0, many times? Especially since your non-vectorized function evaluated the polynomial at win_size
. To be consistent with the non-vectorized method, replace the line
res[win_size-1:] = np.polynomial.polynomial.polyval(np.zeros(len(pnp_array)), fit).T[len(pnp_array) - 1,:]
with
res[win_size-1:] = np.polynomial.polynomial.polyval(win_size, fit[::-1])
then both outputs for the test case are the same.
(That said, I also don't know why you evaluate the polynomial at the right edge of the window; would the middle be the more representative value? But that is something for you to decide.)