我最初发布下面的基准是为了推荐numpy.corrcoef,愚蠢地没有意识到原来的问题已经使用了corrcoef,实际上是在询问高阶多项式拟合。我已经使用statsmodels为多项式r-squared问题添加了一个实际的解决方案,并且我已经离开了原始的基准测试,虽然偏离主题,但对某些人来说可能是有用的。
statsmodels能够直接计算多项式拟合的r^2,这里有2种方法......
import as sm
import 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 (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 (formula, data).fit().rsquared
为了进一步利用statsmodels,还应该查看拟合的模型摘要,可以在Jupyter / IPython笔记本中打印或显示为丰富的HTML表。除了rsquared之外,结果对象还提供对许多有用的统计指标的访问。
model = (y, xpoly)
results = ()
()
以下是我原来的答案,我对各种线性回归r ^ 2方法进行了基准测试...
问题中使用的corrcoef函数仅计算单个线性回归的相关系数r,因此它不能解决高阶多项式拟合的问题r^2。然而,对于它的价值,我发现对于线性回归,它确实是计算r的最快和最直接的方法。
def get_r2_numpy_corrcoef(x, y):
return (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倍,比快12倍。只是为了加强numpy为你做的事情,它比纯蟒蛇快28倍。我并不精通像numba和pypy这样的东西,所以其他人不得不填补这些空白,但我认为这很有说服力,corrcoef是计算简单线性回归的最佳工具。
这是我的基准测试代码。我从Jupyter笔记本上复制粘贴(很难不称它为IPython笔记本......),所以如果有什么事情发生,我道歉。 %timeit magic命令需要IPython。
import numpy as np
from scipy import stats
import as sm
import math
n=1000
x = (1000)*10
()
y = 10 * x + (5+(1000)*10-5)
x_list = list(x)
y_list = list(y)
def get_r2_numpy(x, y):
slope, intercept = (x, y, 1)
r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * (y, ddof=1)))
return r_squared
def get_r2_scipy(x, y):
_, _, r_value, _, _ = (x, y)
return r_value**2
def get_r2_statsmodels(x, y):
return (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 = (sum([(xi-x_bar)**2 for xi in x_list])/(n-1))
y_std = (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))/(x, ddof=1)
zy = ((y))/(y, ddof=1)
r = (zx*zy)/(len(x)-1)
return r**2
def get_r2_numpy_corrcoef(x, y):
return (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)