我最初发布以下基准测试是为了推荐numpy.corrcoef
,愚蠢地没有意识到原始问题已经使用了corrcoef
,实际上是在询问更高阶多项式拟合。 我添加了一个使用statsmodels解决多项式r平方问题的实际解决方案,并保留了原始基准测试,虽然与主题无关,但对某些人可能有用。
statsmodels
有直接计算多项式拟合的r^2
能力,以下是两种方法...
import statsmodels.api as sm
import statsmodels.formula.api as smf
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
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)点上的timeit结果:
- 纯Python (直接
r
计算)
- 1000次循环,3次取最优: 每次1.59毫秒
- Numpy polyfit (适用于n次多项式拟合)
- 1000次循环,3次取最优: 每次326微秒
- Numpy Manual (直接
r
计算)
- 10000次循环,3次取最优: 每次62.1微秒
- Numpy corrcoef (直接
r
计算)
- 10000次循环,3次取最优: 每次56.6微秒
- Scipy (使用线性回归,输出
r
)
- 1000次循环,3次取最优: 每次676微秒
- Statsmodels (可以进行n次多项式和其他许多拟合)
- 1000次循环,3次取最优: 每次422微秒
corrcoef方法比使用numpy方法手动计算r^2要快。 它比polyfit方法快5倍以上,比scipy.linregress快约12倍。 仅为了加强对numpy为您完成的工作的理解,它比纯Python快28倍。 我不太精通numba和pypy之类的东西,因此其他人必须填补这些空白,但我认为这已经足以让我相信corrcoef是计算简单线性回归的最佳工具。
这是我的基准测试代码。 我从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)
2021年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)