최소자승법
Contents
3.4. 최소자승법#
이 절에서는 선형회귀분석을 실제로 수행하여 가중치와 예측치를 계산하는 방법의 하나인 최소자승법(Least Squares)에 대해 설명한다. 최소자승법은 최적화(optimization) 방법을 사용하여 잔차제곱합 RSS(residual sum of squares)를 최소화하는 가중치 벡터 \(\hat{w}\)를 찾아내는 방법이다. 최소자승법에도 여러가지 종류가 있지만 가장 기본적인 방법은 OLS(Ordinary Least Squares)라고 한다. OLS 이외에도 GLS(Generalized Least Squares), WLS(Weighted Least Squares) 등의 방법이 있으면 GLS, WLS에 대해서는 추후 설명하도록 하겠다.
OLS 방법에 의해 찾아낸 최적의 가중치 \(\hat{w}\)는 다음 수식으로 계산한다.
이 식이 성립하려면 행렬 \(X^TX\)의 역행렬 \((X^TX)^{-1}\)이 존재해야 한다. 이 역행렬은 항상 존재하는 것은 아니고 경우에 따라 존재하지 않을 수도 있다. 이 역행렬이 존재하기 위한 조건에 대해서는 추후 자세히 설명하기로 하고 일단은 역행렬이 존재하여 최적 가중치 \(\hat{w}\)를 계산할 수 있다고 가정한다.
3.4.1. 최소자승법의 수학적 유도#
최소자승법을 수학적으로 유도하려면 수학적 최적화(optimzation) 이론과 벡터 및 행렬의 미분을 알고 있어야 한다. 이에 대한 자세한 내용은 저자의 “데이터사이언스스쿨: 수학편”을 참조한다.
예측치 벡터 \(\hat{y}\)는 다음 수식과 같이 독립변수행렬 \(X\)와 가중치 벡터 \(w\)의 곱으로 계산된다.
잔차 벡터 \(e\)는 다음과 같이 정의된다.
이 식을 사용하여 잔차제곱합(RSS:residual sum of squares)을 유도하면 다음과 같다.
최적화의 기본 원리에서 잔차제곱합을 가장 작게하는 가중치는 잔차제곱합 계산식의 기울기 즉, 가중치에 대한 잔차제곱합의 그레디언트(gradient) 벡터를 0으로 하는 가중치를 구하는 것이다.
행렬의 미분을 사용하면 가중치에 대한 잔차제곱합의 그레디언트 벡터는 다음과 같이 계산된다.
따라서 최적가중치 \(\hat{w}\) 에서 다음과 같은 방정식이 성립한다.
또는
이 방정식은 선형회귀분석의 정규방정식(normal equation)이라고 부르는 중요한 방정식이다. 선형회귀분석의 다양한 성질이 이 정규방정식으로부터 유도된다.
만약 \(𝑋^𝑇𝑋\) 행렬의 역행렬이 존재한다면 정규방정식에서 다음처럼 최적가중치 벡터 \(\hat{w}\)를 구할 수 있다.
3.4.2. numpy 패키지를 사용한 최소자승법 가중치 계산과 예측#
앞 절에서 설명한 바와 같이 \(X\) 행렬과 \(y\) 벡터가 준비되어 있다면 numpy 패키지를 사용하여 다음과 같이 최소자승법 가중치를 계산할 수 있다.
여기에서는 팁 데이터에 대해 앞 절에서 준비한 데이터를 사용한다.
import seaborn as sns
from patsy import dmatrix
import numpy as np
tips = sns.load_dataset("tips")
X = np.asarray(dmatrix("scale(total_bill) + scale(size) + C(sex)", tips))
y = tips["tip"].values.reshape(-1, 1)
이 때 독립변수행렬 X를 얻기 위해 DesignMatrix 인스턴스에 numpy asarray 함수를 적용하였고 y의 경우 1차원 벡터를 2차원 행렬 형태로 변형하기 위해 reshape 메서드를 사용하였다.
이제 최적가중치 계산식
을 다음과 같이 적용한다. 이 코드에서 numpy 패캐지의 linalg 서브패키지가 제공하는 역행렬 계산함수 inv를 사용하였다.
from numpy.linalg import inv
hat_w = inv(X.T @ X) @ X.T @ y
hat_w
array([[2.98885891],
[0.02641868],
[0.82551828],
[0.18279436]])
이 값을 사용하여 예측치 \(\hat{y}\)를 다음과 같이 계산한다.
hat_y = X @ hat_w
이 값을 원래의 데이터와 비교해하여 잔차를 구해보자
import pandas as pd
e = y - hat_y
df_result = pd.DataFrame({"y": y.flatten(), "hat_y": hat_y.flatten(), "e": e.flatten()})
df_result
| y | hat_y | e | |
|---|---|---|---|
| 0 | 1.01 | 2.645766 | -1.635766 |
| 1 | 1.66 | 2.194015 | -0.534015 |
| 2 | 3.50 | 3.185475 | 0.314525 |
| 3 | 3.31 | 3.240984 | 0.069016 |
| 4 | 3.61 | 3.737136 | -0.127136 |
| ... | ... | ... | ... |
| 239 | 5.92 | 3.930696 | 1.989304 |
| 240 | 2.00 | 3.592624 | -1.592624 |
| 241 | 2.00 | 3.147135 | -1.147135 |
| 242 | 1.75 | 2.696471 | -0.946471 |
| 243 | 3.00 | 2.812093 | 0.187907 |
244 rows × 3 columns
잔차에 대한 서술적 통계를 계산하면 평균오차는 0에 가깝지만 최대오차는 약 4달러인 것을 알 수 있다.
df_result.e.describe()
count 2.440000e+02
mean -1.401429e-16
std 1.009249e+00
min -2.921205e+00
25% -5.603338e-01
50% -8.777778e-02
75% 5.061815e-01
max 4.045499e+00
Name: e, dtype: float64
3.4.3. statsmodels 패키지를 사용한 최소자승법 가중치 계산과 예측#
statsmodels 패키지는 최소자승법을 포함한 각종 통계계산을 위한 전문 패키지다. statsmodels 패키지를 사용하면 위에서 수행했던 최소자승법 계산이나 예측을 편리하게 할 수 있는 것은 물론이고 이후에 설명할 가중치 검정을 포함한 각종 추가적인 분석도 자동으로 수행해 준다.
statsmodels 패키지를 기반으로 최소자승법을 수행하는 방법은 다음과 같다.
우선 X, y 데이터를 사용하여 OLS 클래스 객체를 만든다. 이 때 데이터의 순서에 주의하라. OLS 클래스 객체를 모형(model) 객체라고도 한다.
import statsmodels.api as sm
model = sm.OLS(y, X)
model
<statsmodels.regression.linear_model.OLS at 0x153832cfe50>
모형 객체는 선형회귀 분석에 필요한 다양한 정보를 속성으로 가지고 있다. 예를 들어 종속변수와 독립변수는 각각 endog_names, exog_names 속성에 저장된다.
model.endog_names
'y'
model.exog_names
['const', 'x1', 'x2', 'x3']
다음으로 fit 메서드를 실행하여 분석결과 객체를 만든다. 이 객체는 OLSResults 클래스 객체다.
result = model.fit()
result
<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x153833549a0>
분석결과 객체의 summary 메서드를 사용하면 가중치를 포함한 분석결과요약을 출력할 수 있다.
print(result.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.468
Model: OLS Adj. R-squared: 0.461
Method: Least Squares F-statistic: 70.36
Date: Fri, 29 Jul 2022 Prob (F-statistic): 1.11e-32
Time: 08:35:14 Log-Likelihood: -347.97
No. Observations: 244 AIC: 703.9
Df Residuals: 240 BIC: 717.9
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 2.9889 0.081 36.737 0.000 2.829 3.149
x1 0.0264 0.137 0.193 0.847 -0.244 0.297
x2 0.8255 0.082 10.104 0.000 0.665 0.986
x3 0.1828 0.081 2.253 0.025 0.023 0.343
==============================================================================
Omnibus: 25.084 Durbin-Watson: 2.098
Prob(Omnibus): 0.000 Jarque-Bera (JB): 47.089
Skew: 0.550 Prob(JB): 5.95e-11
Kurtosis: 4.850 Cond. No. 2.88
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
위 결과요약은 다양한 정보를 표시한다. 이 정보를 모두 설명하기 위해서는 몇가지 이론을 학습해야 하므로 우선 간단히 알 수 있는 부분만 설명하고 나머지는 해당 이론을 설명하면서 같이 설명하도록 하겠다.
우선 가장 상단은 분석제목으로 최소자승법을 사용하는 경우 OLS Regression Results 문자열이 나온다.
다음 부분은 회귀분석에 사용한 데이터와 성능을 표시한다. 상단 왼쪽 열의 정보에서 각각의 제목이 뜻하는 바는 다음과 같다.
Dep. Variable: 종속변수. 이 경우에는y변수다.Model: 사용한 분석모형. 이 경우에는 OLS(Orderary Lease Square)을 사용하였다.Method: 사용한 분석모형. 이 경우에는 최소자승법(Least Squares)을 사용하였다.Date: 분석을 실시한 날짜.Time: 분석을 실시한 시간.No. Observations: 데이터의 개수. 이 경우에는 244개다.Df Residuals: “데이터의 개수 - 독립변수의 개수 - 1”로 계산된 수.Df Model: 사용한 독립변수의 개수. 이 경우에는 total_bill, size, sex 3개다.Covariance Type: 공분산 유형.nonrobust문자열은 이분산성(heteroskedasticity) 모형을 사용하지 않는다는 뜻이로 이분산성에 대해서는 추후 설명한다.
상단의 오른쪽 열은 선형회귀분석의 성능을 표시하는 부분이다. 이 부분에 대해서는 추후 설명한다.
그 다음으로 나오는 것이 가장 중요한 부분으로 추정된 가중치 \(\hat{w}\)를 표시한다.
=======================
coef
-----------------------
const 2.9889
x1 0.0264
x2 0.8255
x3 0.1828
=======================
이 부분은 각 행이 \(\hat{w}_0, \hat{w}_1, \hat{w}_2, \hat{w}_3\)를 표시한다. 즉 다음과 같은 의미다.
result 객체는 가중치를 포함한 다양한 정보를 내포하고 있다. 예를 들어 가중치 값은 params 속성에 들어있다.
result.params
array([2.98885891, 0.02641868, 0.82551828, 0.18279436])
또한 분석에 사용한 모형 객체도 model 속성에 포함하고 있다.
result.model
<statsmodels.regression.linear_model.OLS at 0x153832cfe50>
3.4.4. from_formula 클래스 메서드를 이용한 최소자승법#
위에서 설명한 방법을 쓰려면 먼저 앞절에서 설명한 patsy 패키지와 전처리 방법을 이용하여 X 행렬과 y 벡터를 계산해야 한다.
하지만 statsmodels 패키지는 patsy 패키지를 내장하고 있으며 이를 사용하기 위한 OLS 클래스의 from_formula 클래스 메서드를 제공한다.
OLS 클래스의 from_formula 클래스 메서드를 이용하면 다음과 데이터의 전처리 과정과 최소자승법 계산을 한꺼번에 할 수 있다.
from_formula 클래스 메서드는 모형 문자열과 데이터프레임을 입력으로 가진다.
이 때 모형 문자열은 앞절에서 설명한 것과 비슷하지만 왼쪽에 종속변수 이름과 ~ 기호를 추가해주어야 한다.
model = sm.OLS.from_formula("tip ~ scale(total_bill) + scale(size) + C(sex)", tips)
result = model.fit()
print(result.summary())
OLS Regression Results
==============================================================================
Dep. Variable: tip R-squared: 0.468
Model: OLS Adj. R-squared: 0.461
Method: Least Squares F-statistic: 70.36
Date: Fri, 29 Jul 2022 Prob (F-statistic): 1.11e-32
Time: 08:35:14 Log-Likelihood: -347.97
No. Observations: 244 AIC: 703.9
Df Residuals: 240 BIC: 717.9
Df Model: 3
Covariance Type: nonrobust
=====================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------
Intercept 2.9889 0.081 36.737 0.000 2.829 3.149
C(sex)[T.Female] 0.0264 0.137 0.193 0.847 -0.244 0.297
scale(total_bill) 0.8255 0.082 10.104 0.000 0.665 0.986
scale(size) 0.1828 0.081 2.253 0.025 0.023 0.343
==============================================================================
Omnibus: 25.084 Durbin-Watson: 2.098
Prob(Omnibus): 0.000 Jarque-Bera (JB): 47.089
Skew: 0.550 Prob(JB): 5.95e-11
Kurtosis: 4.850 Cond. No. 2.88
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.