[python] Python과 Numpy를 사용하여 r-squared를 어떻게 계산합니까?

나는 원래 추천 목적으로 아래 벤치 마크를 게시했으며 numpy.corrcoef, 어리석게도 원래 질문이 이미 사용 corrcoef하고 실제로 고차 다항식 적합에 대해 묻는다는 사실을 깨닫지 못했습니다 . 통계 모델을 사용하여 다항식 r- 제곱 질문에 실제 솔루션을 추가했으며 주제를 벗어 났지만 잠재적으로 누군가에게 유용 할 수있는 원래 벤치 마크를 남겼습니다.

---

statsmodelsr^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

를 더 활용하려면 statsmodelsJupyter / 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.59ms
• 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를 계산하는 데 가장 좋은 도구 라고 생각합니다 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)