Python 多重简单线性回归

Posted 2023-03-12

【中文标题】Python 多重简单线性回归

【英文标题】：Python Multiple Simple Linear Regression

【发布时间】：2014-05-04 00:07:58

【问题描述】：

注意，这不是关于多元回归的问题，而是关于在 Python/NumPy (2.7) 中多次进行简单（单变量）回归的问题。

我有两个 m x n 数组 x 和 y。这些行相互对应，每一对都是用于测量的 (x,y) 点的集合。也就是说，plt.plot(x.T, y.T, '.') 将绘制每个 m 个数据集/测量值。

我想知道执行 m 线性回归的最佳方法是什么。目前我遍历行并使用scipy.stats.linregress()。 （假设我不想要基于对矩阵进行线性代数的解决方案，而是想要使用此函数或等效的黑盒函数。）我可以尝试np.vectorize，但文档表明它也会循环。

通过一些实验，我还找到了一种将列表推导与map() 一起使用并获得正确结果的方法。我已经把这两种解决方案都放在下面了。在 IPython 中，`%%timeit` 返回，使用一个小数据集（注释掉）：

(loop) 1000 loops, best of 3: 642 µs per loop
(map) 1000 loops, best of 3: 634 µs per loop

为了尝试放大这一点，我制作了一个更大的随机数据集（维度 trials x trials）：

(loop, trials = 1000)  1 loops, best of 3: 299 ms per loop
(loop, trials = 10000) 1 loops, best of 3: 5.64 s per loop
(map, trials = 1000)   1 loops, best of 3: 256 ms per loop
(map, trials = 10000)  1 loops, best of 3: 2.37 s per loop

对于一个非常大的系列来说，这是一个不错的加速，但我期待更多。有没有更好的办法？

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
np.random.seed(42)
#y = np.array(((0,1,2,3),(1,2,3,4),(2,4,6,8)))
#x = np.tile(np.arange(4), (3,1))
trials = 1000
y = np.random.rand(trials,trials)
x = np.tile(np.arange(trials), (trials,1))
num_rows = shape(y)[0]
slope = np.zeros(num_rows)
inter = np.zeros(num_rows)
for k, xrow in enumerate(x):
    yrow = y[k,:]
    slope[k], inter[k], t1, t2, t3 = stats.linregress(xrow, yrow)
#plt.plot(x.T, y.T, '.')
#plt.hold = True
#plt.plot(x.T, x.T*slope + intercept)
# Can the loop be removed?
tempx = [x[k,:] for k in range(num_rows)]
tempy = [y[k,:] for k in range(num_rows)]
results = np.array(map(stats.linregress, tempx, tempy))
slope_vec = results[:,0]
inter_vec = results[:,1]
#plt.plot(x.T, y.T, '.')
#plt.hold = True
#plt.plot(x.T, x.T*slope_vec + inter_vec)
print "Slopes equal by both methods?: ", np.allclose(slope, slope_vec)
print "Inters equal by both methods?: ", np.allclose(inter, inter_vec)

【问题讨论】：

如果你有一台多核机器，一个简单的事情就是使用 IPythons 的并行map()。结帐ipython.org/ipython-doc/dev/parallel/…

感谢您的指点，没有意识到存在。不幸的是，我使用 IPython 只是为了隔离这段代码。在我的应用程序中，数据集足够小，任何这些技术都可以，我只是想找出最好的方法，以应对不可避免地遇到性能很重要的数据集。

【参考方案1】：

单变量线性回归很简单，可以手动对其进行矢量化：

def multiple_linregress(x, y):
    x_mean = np.mean(x, axis=1, keepdims=True)
    x_norm = x - x_mean
    y_mean = np.mean(y, axis=1, keepdims=True)
    y_norm = y - y_mean

    slope = (np.einsum('ij,ij->i', x_norm, y_norm) /
             np.einsum('ij,ij->i', x_norm, x_norm))
    intercept = y_mean[:, 0] - slope * x_mean[:, 0]

    return np.column_stack((slope, intercept))

一些虚构的数据：

m = 1000
n = 1000
x = np.random.rand(m, n)
y = np.random.rand(m, n)

它在一定程度上优于您的循环选项：

%timeit multiple_linregress(x, y)
100 loops, best of 3: 14.1 ms per loop

【讨论】：

谢谢，确实，手动矢量化相当简单——如果我想花一些时间玩线性代数，我认为可以使用numpy.linalg.lstsq 找到一个有效的解决方案。我的希望是找到一种 Pythonic 方法来矢量化函数 - 我认为 map() 可以解决问题，但从性能来看，它似乎也是通过循环实现的。

循环并不是 Python 中唯一昂贵的东西。如果您对循环进行 cythonized 处理，但在每次迭代时保留一个 Python 函数调用，我认为您不会看到很大的改进。

迟来的接受，感谢您在回答中提供性能数据。