Python多個簡單線性迴歸

Question

請注意這是而不是關於多重回歸的問題，這是一個關於在Python/NumPy（2.7）中多次執行簡單（單變量）迴歸的問題。

我有兩個米 X Ñ陣列x和y。行相互對應，每對是測量的（x，y）點集合。也就是說，plt.plot(x.T, y.T, '.')會繪製數據集/測量結果中的每一個。

我想知道執行線性迴歸的最佳方法是什麼。目前我遍歷行並使用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) 

Answer 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