Chapter 1. 단순회귀분석 (python)

이 글은 Stat 512: Applied Regression Analysis (Siddhartha Nandy) (purdue.edu)  사이트의 Lecture 3 의 SAS 코드를 참고 하여 Python으로 작성된 글입니다.

---

회귀모형을 벡터로 표현하면 다음과 같다.

 

$$ Y = X \beta + \varepsilon, \varepsilon~N_n(0,\sigma^2I_N) $$

 

이때 \(Y = (Y_1,...,Y_N)^T, x_ik = (X)_{ik\bullet}\)

 

Chapter 1. 단순회귀분석 에서 언급한 최소제곱법으로 \(\beta\)를 추정할 때 잔차제곱합(residual sum of squares)을 최소화하는 \(\beta\)를 찾는다. 이 내용을 벡터로 표현하면 다음과 같다.

 

$$ RSS(\beta) = \sum_{i=1}^n \varepsilon_i^2 = (y-X\beta)^T(y-X\beta) $$

 

위 식을 \(\beta\)에 대해 미분하고 0을 만족하는 식을 만들면 정규방정식(normal equation)이 된다.

 

$$ X^TX\beta = X^TY $$

 

X가 full column-rank를 가질때, \(\hat{\beta}\)은 다음과 같다.

 

$$ \hat{\beta} = (X^TX)^{-1}X^Ty $$

 

위에서 X full column-rank 라는것은 X에서 모든 열들이 서로 선형독립임을 의미한다.

 

이 내용을 이용해서 Python 라이브러리로 \(\beta\)를 추정한 값과 직접 계산한 값을 비교해 보도록하자.

1. 패키지 불러오기

import numpy as np
import numpy.linalg as lin
import matplotlib.pyplot as plt
import pandas as pd  
import statsmodels.api as sm  

2. 데이터 입력하기

 

이 데이터는 diamond의 무게와 가격을 이용하여 회귀선을 적합하고 다이아몬드 캐럿의 무게가 0.43g 일때 의 가격을 예측하는 단순회귀분석 문제이다.

X = np.array([.17,.16,.17,.18,.25,.16,.15,.19,
              .21,.15,.18,.28,.16,.20,.23,.29,
              .12,.26,.25,.27,.18,.16,.17,.16,
              .17,.18,.17,.18,.17,.15,.17,.32,
              .32,.15,.16,.16,.23,.23,.17,.33,
              .25,.35,.18,.25,.25,.15,.26,.15])
 
Y = np.array([355,328,350,325,642,342,322,485,
              483,323,462,823,336,498,595,860,
              223,663,750,720,468,345,352,332,
              353,438,318,419,346,315,350,918,
              919,298,339,338,595,553,345,945,
              655,1086,443,678,675,287,693,316])
 
X = X.reshape(48,1)
Y = Y.reshape(48,1)

3. 단순회귀적합

X1 = sm.add_constant(X, has_constant = "add") 
 
model = sm.OLS(Y, X1)
 
fitted_model = model.fit()
 
fitted_model.summary()

결과 값을보면

 

$$ \hat{y} = -259.6259 + 3721.0249 x_1 +\varepsilon  $$

 

으로 적합되는 것을 볼 수 있다.

 

산점도 그리기

pred = fitted_model.predict(X1)
 
plt.scatter(X,Y,label="Ture")                                      
plt.plot(X,pred,label = "predict")         
                                                    
plt.legend()                                        
plt.show()

 

적합된 회귀식과 직접 계산한 \(\hat{\beta} = (X^TX)^-1X^Ty\) 의 값을 비교한결과를 보면 같은것을 확인 할 수 있다..

 

 

 

코드 한눈에 보기

import numpy as np
import numpy.linalg as lin
import matplotlib.pyplot as plt
import pandas as pd  
import statsmodels.api as sm  
 
X = np.array([.17,.16,.17,.18,.25,.16,.15,.19,
              .21,.15,.18,.28,.16,.20,.23,.29,
              .12,.26,.25,.27,.18,.16,.17,.16,
              .17,.18,.17,.18,.17,.15,.17,.32,
              .32,.15,.16,.16,.23,.23,.17,.33,
              .25,.35,.18,.25,.25,.15,.26,.15])
 
Y = np.array([355,328,350,325,642,342,322,485,
              483,323,462,823,336,498,595,860,
              223,663,750,720,468,345,352,332,
              353,438,318,419,346,315,350,918,
              919,298,339,338,595,553,345,945,
              655,1086,443,678,675,287,693,316])
 
X = X.reshape(48,1)
Y = Y.reshape(48,1)
 
X1 = sm.add_constant(X, has_constant = "add") 
 
model = sm.OLS(Y, X1)
 
fitted_model = model.fit()
 
fitted_model.summary()
 
fitted_model.params  
 
pred = fitted_model.predict(X1)
 
plt.scatter(X,Y,label="Ture")                                      
plt.plot(X,pred,label = "predict")         
                                                    
plt.legend()                                        
plt.show()
 
np.dot(np.dot(lin.inv(np.dot(X1.T,X1)),X1.T),Y)