python - 完全矢量化 numpy polyfit

概述

我在使用 polyfit 时遇到了性能问题，因为它似乎无法接受广播数组。我知道 from this post 如果您使用 y ，则相关数据 numpy.polynomial.polynomial.polyfit 可以是多维的。但是，x 维度不能是多维的。有没有办法解决？

动机

我需要计算一些数据的变化率。为了与实验相匹配，我想使用以下方法:获取数据 y 和 x ，对于一小段数据拟合多项式，然后使用拟合系数作为变化率的估计值。

插图

import numpy as np
import matplotlib.pyplot as plt

n = 100
x = np.linspace(0, 10, n)
y = np.sin(x)

window_length = 10
ydot = [np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0]
                                  for j in range(n - window_length)]
x_mids = [x[j+window_length/2] for j in range(n - window_length)]

plt.plot(x, y)
plt.plot(x_mids, ydot)

plt.show()

蓝线是原始数据(正弦曲线)，而绿色是一阶微分(余弦曲线)。

问题

为了对此进行矢量化，我执行了以下操作:

window_length = 10
vert_idx_list = np.arange(0, len(x) - window_length, 1)
hori_idx_list = np.arange(window_length)
A, B = np.meshgrid(hori_idx_list, vert_idx_list)
idx_array = A + B

x_array = x[idx_array]
y_array = y[idx_array]

这将两个一维向量广播到形状为 (n-window_length, window_length) 的二维向量。现在我希望 polyfit 会有一个 axis 参数，这样我就可以并行计算，但没有这样的运气。

有没有人对如何做到这一点有任何建议？我愿意

最佳答案

polyfit 的工作方式是解决以下形式的最小二乘问题:

y = [X].a

其中 y 是您的依赖坐标，[X] 是相应独立坐标的 Vandermonde matrix，而 a 是拟合系数的向量。

在您的情况下，您总是在计算一阶多项式近似值，并且实际上只对一阶项的系数感兴趣。这有一个 well known closed form solution，你可以在任何统计书中找到，或者通过创建一个 2x2 线性方程系统，通过 [X] 的转置预乘上述方程的两边来产生你自己。这一切加起来就是您要计算的值:

>>> n = 10
>>> x = np.random.random(n)
>>> y = np.random.random(n)
>>> np.polyfit(x, y, 1)[0]
-0.29207474654700277
>>> (n*(x*y).sum() - x.sum()*y.sum()) / (n*(x*x).sum() - x.sum()*x.sum())
-0.29207474654700216

最重要的是，您有一个在数据上运行的滑动窗口，因此您可以使用类似于 1D summed area table 的内容，如下所示:

def sliding_fitted_slope(x, y, win):
    x = np.concatenate(([0], x))
    y = np.concatenate(([0], y))
    Sx = np.cumsum(x)
    Sy = np.cumsum(y)
    Sx2 = np.cumsum(x*x)
    Sxy = np.cumsum(x*y)

    Sx = Sx[win:] - Sx[:-win]
    Sy = Sy[win:] - Sy[:-win]
    Sx2 = Sx2[win:] - Sx2[:-win]
    Sxy = Sxy[win:] - Sxy[:-win]

    return (win*Sxy - Sx*Sy) / (win*Sx2 - Sx*Sx)

使用此代码，您可以轻松检查(注意我将范围扩大了 1):

>>> np.allclose(sliding_fitted_slope(x, y, window_length),
                [np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0]
                 for j in range(n - window_length + 1)])
True

和:

%timeit sliding_fitted_slope(x, y, window_length)
10000 loops, best of 3: 34.5 us per loop

%%timeit
[np.polyfit(x[j:j+window_length], y[j:j+window_length], 1)[0]
 for j in range(n - window_length + 1)]
100 loops, best of 3: 10.1 ms per loop

因此，对于您的样本数据，它的速度大约快了 300 倍。

关于python - 完全矢量化 numpy polyfit，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/28237428/

Polyfit

python - 完全矢量化 numpy polyfit