저는 Python과 Numpy를 사용하여 임의 차수의 최적 다항식을 계산하고 있습니다. x 값, y 값 및 내가 맞추려는 다항식의 정도 (선형, 2 차 등)의 목록을 전달합니다.
이 정도는 효과가 있지만 r (상관 계수)과 r 제곱 (결정 계수)도 계산하고 싶습니다. 내 결과를 Excel의 최적 추세선 기능 및 계산하는 r 제곱 값과 비교하고 있습니다. 이것을 사용하여 선형 최적에 대해 r- 제곱을 올바르게 계산하고 있음을 알고 있습니다 (차수가 1과 같음). 그러나 내 함수는 차수가 1보다 큰 다항식에서는 작동하지 않습니다.
Excel은이를 수행 할 수 있습니다. Numpy를 사용하여 고차 다항식에 대한 r 제곱을 어떻게 계산합니까?
내 기능은 다음과 같습니다.
import numpy
# Polynomial Regression
def polyfit(x, y, degree):
results = {}
coeffs = numpy.polyfit(x, y, degree)
# Polynomial Coefficients
results['polynomial'] = coeffs.tolist()
correlation = numpy.corrcoef(x, y)[0,1]
# r
results['correlation'] = correlation
# r-squared
results['determination'] = correlation**2
return results
답변
로부터 numpy.polyfit의 문서, 그것은 선형 회귀를 피팅 있습니다. 특히 차수가 ‘d’인 numpy.polyfit은 평균 함수가있는 선형 회귀에 적합합니다.
E (y | x) = p_d * x ** d + p_ {d-1} * x ** (d-1) + … + p_1 * x + p_0
따라서 해당 피팅에 대한 R- 제곱을 계산하면됩니다. 선형 회귀 에 대한 위키피디아 페이지 는 자세한 내용을 제공합니다. 몇 가지 방법으로 계산할 수있는 R ^ 2에 관심이 있습니다. 가장 쉬운 방법은
SST = Sum(i=1..n) (y_i - y_bar)^2
SSReg = Sum(i=1..n) (y_ihat - y_bar)^2
Rsquared = SSReg/SST
여기서는 ‘y_bar’를 y의 평균으로 사용하고 ‘y_ihat’을 각 점의 적합 값으로 사용합니다.
나는 numpy에별로 익숙하지 않습니다 (보통 R에서 일합니다), 아마도 R 제곱을 계산하는 더 깔끔한 방법이 있지만 다음은 정확해야합니다
import numpy
# Polynomial Regression
def polyfit(x, y, degree):
results = {}
coeffs = numpy.polyfit(x, y, degree)
# Polynomial Coefficients
results['polynomial'] = coeffs.tolist()
# r-squared
p = numpy.poly1d(coeffs)
# fit values, and mean
yhat = p(x) # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y) # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2) # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2) # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot
return results
답변
매우 늦게 답변했지만 누군가 이에 대한 준비 기능이 필요한 경우를 대비하여 :
즉
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
@Adam Marples의 답변에서와 같이.
답변
yanl (아직 다른 라이브러리)에서 sklearn.metrics
갖는 r2_score
기능;
from sklearn.metrics import r2_score
coefficient_of_dermination = r2_score(y, p(x))
답변
나는 이것을 성공적으로 사용하고 있는데, 여기서 x와 y는 배열과 같습니다.
def rsquared(x, y):
""" Return R^2 where x and y are array-like."""
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
return r_value**2
답변
나는 원래 추천 목적으로 아래 벤치 마크를 게시했으며 numpy.corrcoef
, 어리석게도 원래 질문이 이미 사용 corrcoef
하고 실제로 고차 다항식 적합에 대해 묻는다는 사실을 깨닫지 못했습니다 . 통계 모델을 사용하여 다항식 r- 제곱 질문에 실제 솔루션을 추가했으며 주제를 벗어 났지만 잠재적으로 누군가에게 유용 할 수있는 원래 벤치 마크를 남겼습니다.
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 테이블로 인쇄하거나 표시 할 수있는 적합 모델 요약도 살펴 봐야합니다. 결과 개체는뿐만 아니라 많은 유용한 통계 메트릭에 대한 액세스를 제공합니다 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)
답변
R- 제곱은 선형 회귀에만 적용되는 통계입니다.
기본적으로 선형 회귀로 설명 할 수있는 데이터의 변동 정도를 측정합니다.
따라서 평균에서 각 결과 변수의 총 제곱 편차 인 “총 제곱합”을 계산합니다. . .
\ sum_ {i} (y_ {i}-y_bar) ^ 2
여기서 y_bar는 y의 평균입니다.
그런 다음 FITTED 값이 평균과 얼마나 다른지 “회귀 제곱합”을 계산합니다.
\ sum_ {i} (yHat_ {i}-y_bar) ^ 2
그리고 그 둘의 비율을 찾으십시오.
이제 다항식 피팅을 위해해야 할 일은 해당 모델의 y_hat을 연결하는 것뿐입니다.하지만 r- 제곱이라고 부르는 것은 정확하지 않습니다.
여기 에 내가 찾은 링크가 있습니다.