有效扩展线性回归

Question

我有一个包含两列（X 和 Y 坐标）的数据框。我需要从 df 开始的扩展线性回归。例如，在第 2 点，我需要前 2 点的回归；在第 3 点，前 3 点我需要它，依此类推。根据文档，对于具有一个输入和一个输出的函数 expanding_apply 可以使用，但 linregress 提供 5 个输出。

目前我正在对所有行进行 for 循环，它工作正常，但速度非常慢，几乎无法使用。

我尝试了几件事，但遭到拒绝。尝试将输入作为元组发送：

pd.expanding_apply((df.x, df.y), linregress)
> AttributeError: 'tuple' object has no attribute 'dtype'

尝试将输入作为 df 发送：

pd.expanding_apply(df[['x','y']], linregress)
> IndexError: tuple index out of range

每个回归几乎与之前的回归相同（仅添加了一个数据点），因此大概有很大的加速空间。有什么方法可以在 pandas 或 numpy/scipy 域中更有效地实现这一点？

编辑：linregress 可以选择接受一个二维数组（而不是 2 个单独的一维数组），因此 linregress(df[['x','y']]) 本身可以正常工作。然而expanding_apply 可能期待 arg 的系列，而不是 df。

Answer 1

要计算y = a*x + b 的线性回归参数，您必须求解一个超定方程组X*a = y，其中：

X = [[1, x0], [1, x1], ..., [1, x(n-1)]]
a = [b, a]
y = [ y0, y1, ..., y(n-1)]

如果你只是在a和b的值之后，你可以将你的系统两边都预乘X.T，然后求解得到的2x2系统。仔细观察，得到的数组可以写成：

np.dot(X.T, X) = [[n, np.sum(x)],
                  [np.sum(x), np.sum(x*x)]]
np.dot(X.T, y) = [np.sum(y), np.sum(x*y)]

将所有这些放在一起，给定两个长度相等的一维数组 x 和 y，您可以使用 numpy >= 1.8 执行以下操作：

n = 10
x, y = np.random.rand(2, n)

lhs = np.empty((n-1, 2, 2))
rhs = np.empty((n-1, 2))

lhs[:, 0, 0] = np.arange(2, n+1)
lhs[:, 0, 1] = np.cumsum(x)[1:]
lhs[:, 1, 0] = lhs[:, 0, 1]
lhs[:, 1, 1] = np.cumsum(x*x)[1:]

rhs[:, 0] = np.cumsum(y)[1:]
rhs[:, 1] = np.cumsum(x*y)[1:]

a = np.linalg.solve(lhs, rhs)

您可以检查a 与polyfit 的结果相比是否持有正确的参数：

In [49]: a
Out[49]:
array([[ 0.64778976, -0.39918768],
       [ 0.76225593, -0.41054035],
       [ 0.72598372, -0.35430181],
       [ 0.70608159, -0.33873589],
       [ 0.6899674 , -0.34941498],
       [ 0.68270772, -0.34834723],
       [ 0.71031366, -0.59487271],
       [ 0.7422803 , -0.74757567],
       [ 0.65982282, -0.48593478]])

In [50]: for j in range(2, n+1):
   ....:     print np.polynomial.polynomial.polyfit(x[:j], y[:j], 1)
   ....:
[ 0.64778976 -0.39918768]
[ 0.76225593 -0.41054035]
[ 0.72598372 -0.35430181]
[ 0.70608159 -0.33873589]
[ 0.6899674  -0.34941498]
[ 0.68270772 -0.34834723]
[ 0.71031366 -0.59487271]
[ 0.7422803  -0.74757567]
[ 0.65982282 -0.48593478]