向量化 numpy polyfit 以移动 window

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 步骤

In [304]: def nppolyfit(pnp_array, **kwargs):
     ...:   """ Moving polyfit 
     ...:   """
     ...:   win_size = kwargs['length']
     ...:   degree = kwargs['degree']
     ...: 
     ...:   xdata = np.arange(win_size)
     ...:   res = np.zeros(pnp_array.shape)
     ...:   fits = []
     ...:   for i in range(win_size, len(pnp_array) + 1):
     ...:       fit = np.polyfit(xdata , pnp_array[i - win_size : i], deg = degr
     ...: ee)
     ...:       res[i-1] = np.poly1d(fit)(win_size)
     ...:       fits.append(fit)
     ...:   return res, fits
     ...: 
     ...: 
In [305]: 
In [305]: nppolyfit(npd,length=win_size1, degree=degree)
Out[305]: 
(array([  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.]),
 [array([  1.00000000e+00,   1.60677522e-15]),
  array([ 1.,  1.]),
  array([ 1.,  2.]),
  array([ 1.,  3.]),
  array([ 1.,  4.]),
  array([ 1.,  5.]),
  array([ 1.,  6.]),
  array([ 1.,  7.]),
  array([ 1.,  8.]),
  array([ 1.,  9.]),
  array([ 1., 10.]),
  array([  1.,  11.]),
  array([  1.,  12.]),
  array([  1.,  13.]),
  array([  1.,  14.]),
  array([  1.,  15.]),
  array([  1.,  16.]),
  array([  1.,  17.]),
  array([  1.,  18.]),
  array([  1.,  19.])])

然后我应该将其与多维 polyfit 案例中的 fit 变量进行比较。

将您的最后一个函数更改为 return fit 以及 res:

In [308]: nppolyfitv(npd, length=win_size1, degree=degree)
Out[308]: 
(array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,
         1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
         1.,  1.,  1.,  1.]),
 array([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           1.00000000e+00,   1.00000000e+00],
        [  1.60677522e-15,   1.00000000e+00,   2.00000000e+00,
           3.00000000e+00,   4.00000000e+00,   5.00000000e+00,
           6.00000000e+00,   7.00000000e+00,   8.00000000e+00,
           9.00000000e+00,   1.00000000e+01,   1.10000000e+01,
           1.20000000e+01,   1.30000000e+01,   1.40000000e+01,
           1.50000000e+01,   1.60000000e+01,   1.70000000e+01,
           1.80000000e+01,   1.90000000e+01]]))

拟合似乎匹配。所以大概问题出在 poly1d 步骤。

Answer 2

polyfit方法returns多项式系数，最高次幂优先。

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

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

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

向量化 numpy polyfit 以移动 window

Vectorizing numpy polyfit for moving window

python

numpy

polynomials