我最初发布下面的基准是为了推荐numpy.corrcoef，愚蠢地没有意识到原始问题已经使用corrcoef并且实际上是在询问高阶多项式拟合。我已经使用 statsmodels 为多项式 r-squared 问题添加了一个实际的解决方案，并且我留下了原始基准，这些基准虽然离题，但可能对某人有用。
statsmodels具有r^2直接计算多项式拟合的能力，这里有两种方法......
import statsmodels.api as sm
import statsmodels.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 # or rsquared_adj
要进一步利用statsmodels，还应该查看拟合模型摘要，它可以在 Jupyter/IPython 笔记本中打印或显示为丰富的 HTML 表格。除了 .results 对象之外，结果对象还提供对许多有用统计指标的访问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) 点的一堆方法的 timeit 结果：
- 纯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 为您所做的事情，它比纯 python 快 28 倍。我不精通 numba 和 pypy 之类的东西，所以其他人必须填补这些空白，但我认为这对我来说很有说服力，它是计算简单线性回归corrcoef的最佳工具。r
这是我的基准测试代码。我从 Jupyter Notebook 复制粘贴（很难不称它为 IPython Notebook...），所以如果途中出现任何问题，我深表歉意。%timeit 魔术命令需要 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_list)
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)
21 年 7 月 28 日基准测试结果。（Python 3.7、numpy 1.19、scipy 1.6、statsmodels 0.12）
Python
2.41 ms ± 180 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numpy polyfit
318 µs ± 44.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Numpy Manual
79.3 µs ± 4.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Numpy corrcoef
83.8 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Scipy
221 µs ± 7.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Statsmodels
375 µs ± 3.63 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
