我最初發布下面的基準是為了推薦numpy.corrcoef，愚蠢地沒有意識到原來的問題已經使用了corrcoef，實際上是在詢問高階多項式擬合。我已經使用statsmodels為多項式r-squared問題添加了一個實際的解決方案，并且我已經離開了原始的基準測試，雖然偏離主題，但對某些人來說可能是有用的。

statsmodels能夠直接計算多項式擬合的r^2，這里有2種方法......

import statsmodels.api as sm

import stasmodels.formula.api as smf

# Construct the columns for the different powers of x

def get_r2_statsmodels(x, y, k=1):

xpoly = np.column_stack([x**i for i in range(k+1)])

return sm.OLS(y, xpoly).fit().rsquared

# Use the formula API and construct a formula describing the polynomial

def get_r2_statsmodels_formula(x, y, k=1):

formula = 'y ~ 1 + ' + ' + '.join('I(x**{})'.format(i) for i in range(1, k+1))

data = {'x': x, 'y': y}

return smf.ols(formula, data).fit().rsquared

為了進一步利用statsmodels，還應該查看擬合的模型摘要，可以在Jupyter / IPython筆記本中打印或顯示為豐富的HTML表。除了rsquared之外，結果對象還提供對許多有用的統計指標的訪問。

model = sm.OLS(y, xpoly)

results = model.fit()

results.summary()

以下是我原來的答案，我對各種線性回歸r ^ 2方法進行了基準測試...

問題中使用的corrcoef函數僅計算單個線性回歸的相關系數r，因此它不能解決高階多項式擬合的問題r^2。然而，對于它的價值，我發現對于線性回歸，它確實是計算r的最快和最直接的方法。

def get_r2_numpy_corrcoef(x, y):

return np.corrcoef(x, y)[0, 1]**2

這些是我通過比較1000個隨機（x，y）點的方法得到的時間結果：

純Python（直接r計算）

1000循環，最佳3：每循環1.59毫秒

Numpy polyfit（適用于n次多項式擬合）

1000個循環，每循環最佳3：326μs

Numpy手冊（直接r計算）

10000循環，最佳3：每循環62.1μs

Numpy corrcoef（直接r計算）

10000循環，最佳3：每循環56.6μs

Scipy（線性回歸以r為輸出）

1000個循環，最佳3：676μs/循環

Statsmodels（可以做n次多項式和許多其他擬合）

1000個循環，最佳3：每循環422μs

corrcoef方法使用numpy方法以“手動”方式勉強計算r ^ 2。它比polyfit方法快5倍，比scipy.linregress快12倍。只是為了加強numpy為你做的事情，它比純蟒蛇快28倍。我并不精通像numba和pypy這樣的東西，所以其他人不得不填補這些空白，但我認為這很有說服力，corrcoef是計算簡單線性回歸的最佳工具。

這是我的基準測試代碼。我從Jupyter筆記本上復制粘貼（很難不稱它為IPython筆記本......），所以如果有什么事情發生，我道歉。 ％timeit magic命令需要IPython。

import numpy as np

from scipy import stats

import statsmodels.api as sm

import math

n=1000

x = np.random.rand(1000)*10

x.sort()

y = 10 * x + (5+np.random.randn(1000)*10-5)

x_list = list(x)

y_list = list(y)

def get_r2_numpy(x, y):

slope, intercept = np.polyfit(x, y, 1)

r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * np.var(y, ddof=1)))

return r_squared

def get_r2_scipy(x, y):

_, _, r_value, _, _ = stats.linregress(x, y)

return r_value**2

def get_r2_statsmodels(x, y):

return sm.OLS(y, sm.add_constant(x)).fit().rsquared

def get_r2_python(x_list, y_list):

n = len(x)

x_bar = sum(x_list)/n

y_bar = sum(y_list)/n

x_std = math.sqrt(sum([(xi-x_bar)**2 for xi in x_list])/(n-1))

y_std = math.sqrt(sum([(yi-y_bar)**2 for yi in y_list])/(n-1))

zx = [(xi-x_bar)/x_std for xi in x_list]

zy = [(yi-y_bar)/y_std for yi in y_list]

r = sum(zxi*zyi for zxi, zyi in zip(zx, zy))/(n-1)

return r**2

def get_r2_numpy_manual(x, y):

zx = (x-np.mean(x))/np.std(x, ddof=1)

zy = (y-np.mean(y))/np.std(y, ddof=1)

r = np.sum(zx*zy)/(len(x)-1)

return r**2

def get_r2_numpy_corrcoef(x, y):

return np.corrcoef(x, y)[0, 1]**2

print('Python')

%timeit get_r2_python(x_list, y_list)

print('Numpy polyfit')

%timeit get_r2_numpy(x, y)

print('Numpy Manual')

%timeit get_r2_numpy_manual(x, y)

print('Numpy corrcoef')

%timeit get_r2_numpy_corrcoef(x, y)

print('Scipy')

%timeit get_r2_scipy(x, y)

print('Statsmodels')

%timeit get_r2_statsmodels(x, y)