python - 为移动窗口向量化 numpy polyfit

Question

我想向量化以下函数

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

到目前为止做了什么:

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

但看起来我遗漏了什么或做错了什么。你能纠正我吗？也许还有另一种更有效的解决方案？谢谢。

测试用例:

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)

结果是:

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

Answer 1

polyfit方法返回多项式系数，最高幂优先。

polyval方法期望系数具有最低幂优先。将一种方法的输出提供给另一种方法时，请考虑到这一点。

此外，polyval 的 x 参数不合逻辑:np.zeros(len(pnp_array))。为什么要多次要求 polyval 计算同一点 0 处的多项式？特别是因为您的非矢量化函数计算了 win_size 处的多项式。为了与非矢量化方法保持一致，将行替换为

res[win_size-1:] = np.polynomial.polynomial.polyval(np.zeros(len(pnp_array)), fit).T[len(pnp_array) - 1,:]

与

res[win_size-1:] = np.polynomial.polynomial.polyval(win_size, fit[::-1])

那么测试用例的两个输出是相同的。

(也就是说，我也不知道你为什么要在窗口的右边缘计算多项式；中间的值会更具代表性吗？但这由你决定。)