Search code examples
pythonnumpypolynomials

Vectorizing numpy polyfit for moving window


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.]

Solution

  • 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.)