Первоначально я разместил тесты ниже с целью рекомендации numpy.corrcoef, глупо не понимая, что в исходном вопросе уже используется corrcoef и фактически спрашивал о полиномиальных подстановках более высокого порядка. Я добавил фактическое решение к полиномиальному вопросу r-squared, используя statsmodels, и оставил исходные тесты, которые, хотя и не по теме, потенциально полезны для кого-то.
statsmodels имеет возможность напрямую вычислять r^2 полиномиального подбора, вот 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, следует также ознакомиться с краткой сводкой модели, которая может быть напечатана или отображена в виде таблицы HTML в записной книжке Jupyter / IPython. Объект результатов обеспечивает доступ ко многим полезным статистическим метрикам в дополнение к 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) точек:
- Pure Python (прямой
r расчет)
- 1000 циклов, лучшее из 3: 1,59 мс на цикл
- Полифит Numpy (применим к полиномам n-й степени)
- 1000 петель, лучшее из 3: 326 мкс на петлю
- Numpy Manual (прямой
r расчет)
- 10000 петель, лучшее из 3: 62,1 мкс на петлю
- Numpy corrcoef (прямой
r расчет)
- 10000 петель, лучшее из 3: 56,6 мкс на петлю
- Сципи (линейная регрессия с
r в качестве выхода)
- 1000 петель, лучшее из 3: 676 мкс на петлю
- Statsmodels (может делать полином n-й степени и много других подгонок)
- 1000 петель, лучшее из 3: 422 мкс на петлю
Метод corrcoef едва ли превосходит вычисление r ^ 2 «вручную», используя простые методы. Это в 5 раз быстрее, чем метод polyfit, и в 12 раз быстрее, чем scipy.linregress. Просто чтобы подчеркнуть, что NumPy делает для вас, это в 28 раз быстрее, чем чистый Python. Я не очень разбираюсь в таких вещах, как numba и pypy, поэтому кто-то другой должен был бы заполнить эти пробелы, но я думаю, что это достаточно убедительно для меня, что corrcoef - лучший инструмент для вычисления r для простой линейной регрессии.
Вот мой контрольный код. Я вставил копию с ноутбука Jupyter (трудно не назвать его IPython Notebook ...), поэтому я прошу прощения, если что-то сломалось на пути. Команде% 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)