我有屏蔽数组的 3D 堆栈。我想对每一行的值执行线性回归，col(空间索引)沿轴 0(时间)。这些堆栈的尺寸各不相同，但典型的形状可能是 (50, 2000, 2000)。我的空间有限但时间密集的测试用例具有以下维度:
stack.ma_stack.shape
(1461, 390, 327)
我对每一行进行了快速循环测试，列:
from scipy.stats.mstats import linregress
#Ordinal dates
x = stack.date_list_o
#Note: idx should be row, col
def sample_lstsq(idx):
b = stack.ma_stack[:, idx[0], idx[1]]
#Note, this is masked stats version
slope, intercept, r_value, p_value, std_err = linregress(x, b)
return slope
out = np.zeros_like(stack.ma_stack[0])
for row in np.arange(stack.ma_stack.shape[1]):
for col in np.arange(stack.ma_stack.shape[2]):
out[row, col] = sample_lstsq((row, col))
这有效(缓慢)。我知道必须有更有效的方法。
我开始尝试使用索引数组和 np.vectorize，但我认为这实际上不会提供任何真正的改进。我考虑过将所有东西转储到 Pandas 或尝试移植到 Cython，但我希望我能坚持使用 NumPy/SciPy。或者也许并行解决方案是我提高性能的最佳选择？
此外，是否有人对 NumPy/SciPy 线性回归选项进行了基准测试？我遇到过以下选项，但我自己没有测试过:
- scipy.stats.linregress
- numpy.linalg.leastsq
- numpy.polyfit(deg=1)
我希望有一种现有方法可以显着提高性能，而无需进行大量实现工作。谢谢。
编辑 12/3/13 @02:29
@HYRY 建议的方法非常适用于上述样本数据集(运行时间约 15 秒)，该数据集在所有维度(空间和时间)上都是连续的(未屏蔽)。但是，当将包含缺失数据的掩码数组传递给 np.linalg.leastsq 时，所有掩码值都将填充 fill_value(默认为 1E20)，这会导致虚假的线性拟合。
幸运的是，numpy masked array 模块有 np.ma.polyfit(deg=1)，它可以像 np.linalg.leastsq 一样处理 2D y 数组。查看源代码(https://github.com/numpy/numpy/blob/v1.8.0/numpy/ma/extras.py#L1852)，ma polyfit 只是 np.polyfit 的包装器，它使用 x 和 y 掩码的组合掩码。当 y 中缺失数据的位置不变时，这对 2D y 非常有效。
不幸的是，我的数据在空间和时间上有可变的缺失数据位置。这是来自另一个堆栈的示例:
In [146]: stack.ma_stack.shape
Out [146]: (57, 1889, 1566)
对单个索引进行采样会返回具有 6 个未屏蔽值的时间序列:
In [147]: stack.ma_stack[:,0,0]
Out [147]:
masked_array(data = [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
519.7779541015625 -- -- -- 518.9047241210938 -- -- -- -- -- -- --
516.6539306640625 516.0836181640625 515.9403686523438 -- -- -- --
514.85205078125 -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --],
mask = [ True True True True True True True True True True True True
True True True True True True True True True False True True
True False True True True True True True True False False False
True True True True False True True True True True True True
True True True True True True True True True],
fill_value = 1e+20)
对不同的位置进行采样会从不同的时间片返回不同数量的未屏蔽值:
In [148]: stack.ma_stack[:,1888,1565]
Out[148]:
masked_array(data = [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
729.0936889648438 -- -- -- 724.7155151367188 -- -- -- -- -- -- --
722.076171875 720.9276733398438 721.9603881835938 -- 720.3294067382812 --
-- 713.9591064453125 709.8037719726562 707.756103515625 -- -- --
703.662353515625 -- -- -- -- 708.6276245117188 -- -- -- -- --],
mask = [ True True True True True True True True True True True True
True True True True True True True True True False True True
True False True True True True True True True False False False
True False True True False False False True True True False True
True True True False True True True True True],
fill_value = 1e+20)
每个索引的最小未屏蔽值数量为 6，最大数量为 45。因此每个位置至少有一些屏蔽值。
作为引用，我的 x(时间序数)值都是未屏蔽的:
In [150]: stack.date_list_o
Out[150]:
masked_array(data = [ 733197.64375 733962.64861111 733964.65694444 733996.62361111
733999.64236111 734001.63541667 734033.64305556 734071.64722222
734214.675 734215.65694444 734216.625 734226.64722222
734229.63819444 734232.65694444 734233.67847222 734238.63055556
734238.63055556 734245.65277778 734245.65277778 734255.63125
734255.63125 734307.85 734326.65138889 734348.63888889
734348.63958333 734351.85 734363.70763889 734364.65486111
734390.64722222 734391.63194444 734394.65138889 734407.64652778
734407.64722222 734494.85 734527.85 734582.85
734602.65486111 734664.85555556 734692.64027778 734741.63541667
734747.85 734807.85555556 734884.85555556 734911.65763889
734913.64375 734917.64236111 734928.85555556 734944.71388889
734961.62777778 735016.04583333 735016.62777778 735016.85555556
735036.65347222 735054.04583333 735102.63125 735119.61180556
735140.63263889],
mask = False,
fill_value = 1e+20)
所以我 reshape stack.ma_stack 并运行 polyfit:
newshape = (stack.ma_stack.shape[0], stack.ma_stack.shape[1]*stack.ma_stack.shape[2])
print newshape
#(57, 2958174)
y = stack.ma_stack.reshape(newshape)
p = np.ma.polyfit(x, y, deg=1)
但是到了 ~1500 列，y 中的每一行都被“累积”屏蔽了，我得到了一些提示和空输出:
RankWarning: Polyfit may be poorly conditioned
** On entry to DLASCL, parameter number 4 had an illegal value
...
因此，看起来在不同位置使用具有缺失数据的 2D y 是行不通的。我需要一个 leastsq 拟合，它在每个 y 列中使用所有可用的未屏蔽数据。可能有一种方法可以通过仔细压缩 x 和 y 并跟踪未屏蔽的索引来做到这一点。
还有其他想法吗？ pandas 开始看起来像是一个很好的解决方案。
13 年 12 月 3 日 @20:29 编辑
@HYRY 提供了一种解决方案，适用于时间(轴=0)维度中的缺失值。我不得不稍微修改以处理空间 (axes=1,2) 维度中的缺失值。如果一个特定的空间索引在时间上只有一个未屏蔽的条目，我们当然不想尝试线性回归。这是我的实现:
def linreg(self):
#Only compute where we have n_min unmasked values in time
n_min = 3
valid_idx = self.ma_stack.count(axis=0).filled(0) >= n_min
#Returns 2D array of unmasked columns
y = self.ma_stack[:, valid_idx]
#Extract mask for axis 0 - invert, True where data is available
mask = ~y.mask
#Remove masks, fills with fill_value
y = y.data
#Independent variable is time ordinal
x = self.date_list_o
x = x.data
#Prepare matrices and solve
X = np.c_[x, np.ones_like(x)]
a = np.swapaxes(np.dot(X.T, (X[None, :, :] * mask.T[:, :, None])), 0, 1)
b = np.dot(X.T, (mask*y))
r = np.linalg.solve(a, b.T)
#Create output grid with original dimensions
out = np.ma.masked_all_like(self.ma_stack[0])
#Fill in the valid indices
out[valid_idx] = r[:,0]
运行时非常快 - 此处讨论的数组维度仅需 5-10 秒。
最佳答案
如果我没理解错的话，你想做很多线性回归y = k * x + b，只有一个x，但是很多y，对于您要计算的每个 y 和 b。
如果x的形状是50，y是(50, 1000)，你可以使用numpy.linalg.lstsq，这里是一些演示:
import numpy as np
x = np.random.rand(50)
k = np.random.rand(1000)
b = np.random.rand(1000)
y = np.outer(x, k) + b + np.random.normal(size=(50, 1000), scale=1e-10)
r = np.linalg.lstsq(np.c_[x, np.ones_like(x)], y)[0]
print np.allclose(r[0], k)
print np.allclose(r[1], b)
对于形状为 (50, 2000, 2000) 的 y，您可以将其 reshape 为 (50, 2000*2000)。
编辑
这是我的屏蔽数组解决方案。公式为:
将 Y 准备为形状为 (1889*1566, 57) 的二维数组，将 X 准备为形状为 (57, 2) 的二维数组。 mask 为与 Y 形状相同的 bool 数组，True 表示 Y 中的值可用。
计算形状为(1889*1566, 2, 2)的数组a，形状为(1889*1566, 2)的b，然后调用numpy。 linalg.solve(a, b)，这里是一些演示代码:
import numpy as np
N = 50
M = 1000
x = np.random.rand(N)
X = np.c_[x, np.ones_like(x)]
beta = np.random.rand(M, 2)
Y = np.dot(beta, X.T)
Y += np.random.normal(scale=0.1, size=Y.shape)
mask = np.random.randint(0, 2, size=Y.shape).astype(np.bool)
a = np.swapaxes(np.dot(X.T, (X[None, :, :] * mask[:, :, None])), 0, 1)
b = np.dot(X.T, (mask*Y).T)
beta2 = np.linalg.solve(a, b.T)
i = 123
print "real:", beta[i]
print "by solve:", beta2[i]
m = mask[i]
x2 = X[m]
y2 = Y[i, m]
print "by lstsq:", np.linalg.lstsq(x2, y2)[0]
输出第123个系数:
real: [ 0.35813131 0.29736779]
by solve: [ 0.38088499 0.30382547]
by lstsq: [ 0.38088499 0.30382547]
你也可以通过下面的代码计算a数组，我认为它会比上面的方法占用更少的内存:
a2 = np.empty((M, 2, 2))
xm = mask * x
a2[:, 0, 0] = (xm*xm).sum(1)
a2[:, 1, 0] = (xm*mask).sum(1)
a2[:, 0, 1] = a2[:, 1, 0]
a2[:, 1, 1] = (mask).sum(1)
print np.allclose(a2, a)
关于python - 3D numpy 数组的每个元素的高效一维线性回归，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/20343500/