[data analysis] 회귀계수의 추정: 최소제곱법(Least Square method)

최소자승법(Least Square method)

식 1의 형태인 회귀선(회귀방정식)은 설명변수에 대응하는 반응변수의 관계를 설명하기 위해 통계적으로 추정된 방정식입니다. 이 모형으로부터의 예측값 역시 통계적으로 추정된 값으로 실제로 관측된 값과는 차이가 존재합니다. 그러므로 추정값과 관측값과의 차이를 평가하여 모형의 적합성 여부를 결정할 필요가 있습니다.

y = β0 + β1x + ε	(식 1)
x: 설명변수
y: 반응변수
β0 : 편차
β1 : 회귀계수(가중치)
ε : 오차

모형의 구성요소인 편차 β₀와 회귀계수 β₁은 모집단의 회귀모형에 대한 것으로 미지의 값(unknown value)인 모수이므로 표본의 통계량으로부터 추정되어야 합니다. 이를 구분하기 위해 표본집단의 편차와 회귀계수를 각각 b₀와 b₁로 나타내며 모수를 추정하기 위한 불편추정치(unbiased estimator)로 사용합니다. 이 추정치들 중 관측치와 실측치의 차이인 오차(error, e) 또는 잔차(residual)는 식 2와 같이 계산됩니다.

\begin{align} e & = y − (b_0+ b_1x)\\ &=y − \hat{y}\\& e: \text{오차 또는 잔차}\\ &\hat{y}: \text{추정치}\end{align}

(식 2)

자료에서 발생하는 각 샘플의 오차는 음수와 양수 모두 가능하므로 그들의 합은 0에 근접하기 때문에 회귀모형의 적합성을 위한 판단근거로 사용할 수 없습니다. 대신에 각 오차의 절대값이나 제곱값들의 합을 사용합니다. 절대값을 사용하는 경우는 1차 식이되며 제곱을 적용한 경우는 2차 식이 됩니다.

회귀모형을 나타내는 회귀선(기울기와 편차)은 다양하게 나타낼 수 있습니다(회귀분석의 정의와 가정의 그림 2 참조). 오차는 회귀식과 설명변수에 반응하는 결과로서 궁극적으로 최소의 오차를 생성하는 회귀식이 최적의 모형이 됩니다. b₀와 b₁을 미지수로 하는 오차함수(식 2)의 극소점은 그 함수의 미분으로 결정할 수 있으며 미분을 적용하기 위한 함수는 1차 보다는 2차식이 유리합니다. 그러므로 오차 함수로 절대값 보다는 식 3과 같이 제곱합(SS)을 사용하며 이를 오차제곱합(Sum of Square Error, SSE) 또는 잔차제곱합(Residual of Square Sum, RSS)라고 합니다.

\begin{align}\text{SSE}& = \sum^n_{i=1}e^2\\& = \sum^n_{i=1}(y-\hat{y})^2 \\&= \sum^n_{i=1}\left(y-(b_0+b_1x)^2\right)^2\\n, \hat{y}:\;&\text{표본의 크기, 추정치}\end{align}

(식 3)

식 3으로부터 SSE를 최소로 하기 위한 b₀와 b₁를 계산할 수 있습니다(식 4). 이 최소점은 그림 1에서 나타낸 것과 같이 회귀계수와 편차항의 계수에 대한 SSE의 미분이 0이 되는 극점에서 발견될 것입니다.

그림 1. SSE 함수의 미분.

x=np.linspace(-2, 2, 100)
y=x**2
x1=np.linspace(-2, 0, 50)
y1=-2*x1-1
x2=np.linspace(0, 2, 50)
y2=2*x2-1
plt.figure(figsize=(4,3))
plt.plot(x, y, color="g", label=r"SSE=$(y_i-bx_i)^2$")
plt.plot(x1, y1, color="brown", ls="--", label=r"$\frac{d\;SSE}{dx}>0$" )
plt.plot(x2, y2, color="b",  ls="--", label=r"$\frac{d\;SSE}{dx}<0$" )
plt.hlines(0, -2, 2, ls="--", color="k", label=r"$\frac{d\;SSE}{dx}=0$")
plt.xlabel("x")
plt.ylabel("SSE")
plt.legend(loc="best",prop={'size':10}, labelcolor='linecolor', frameon=False)
plt.show()

그림 1의 모델은 편차항을 고려하지 않은 것으로 편차항을 고려한다면 미분은 각 계수에 따라 이루어지므로 식 4와 같이 편미분으로 각 계수를 결정할 수 있습니다.

\begin{align}\frac{\partial{\text{sse}}}{\partial b_0} & =0, \quad \frac{\partial{\text{sse}}}{\partial b_1} =0\\ b_1& = \frac{\sum^n_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^n_{i=1}(x_i-\bar{x})^2} \\&=\frac{S_{xy}}{S_{xx}}\\ b_0& = \frac{y-b_1x}{n} \\ & n: \text{자료의 크기}\\& x, y: \text{설명변수, 반응변수}\end{align}

(식 4)

식 4와 같이 오차가 최소가 되도록 회귀계수를 추정하는 것을 최소자승법(Least Square method) 또는 Ordinary Least Square(OLS)이라 합니다.

예 1)

kospi 지수의 일일 주가 자료중 시가(Open)을 설명변수로 하여 종가(Close)를 추정하는 회귀모델을 작성합니다.

	Open	Close
0	2874.50	2944.45
1	2943.67	2990.57
2	2993.34	2968.21
⋮	⋮	⋮

다음 코드는 분석을 위한 자료를 호출하기 위한 것입니다.

st=pd.Timestamp(2021,1, 1)
et=pd.Timestamp(2024, 5, 10)
kos=fdr.DataReader('KS11',st, et)[["Open","Close"]]
kos.index=range(len(kos))
kos.head(3).round(2) 

	Open	Close
0	2201.21	2175.17
1	2192.58	2176.46
2	2154.97	2155.07

위 자료를 표준화합니다.

X=kos.values[:,0].reshape(-1,1)
y=kos.values[:,1].reshape(-1,1)
from sklearn.preprocessing import StandardScaler 
#독립변수 정규화(표준화)
xScaler=StandardScaler().fit(X)
X_n=xScaler.transform(X)
#반응변수 정규화(표준화)
yScaler=StandardScaler().fit(y)
y_n=yScaler.transform(y)

b1=(np.sum((X_n-X_n.mean())*(y_n-y_n.mean())))/np.sum((X_n-X_n.mean())**2)
b0=np.sum(y_n-b1*X_n)/len(y_n)
print("편차: %.3f, 회귀계수: %.3f" %(b0, b1))

편차: -0.000, 회귀계수: 0.997

위 결과는 OLS()와 LinearRegression() 클래스를 사용한 결과와 같습니다(단순회귀 분석의 예 1 참조).

1개의 설명변수와 1개의 반응변수를 가지는 단순회귀모형일 경우 위와 같이 계산이 쉽게 이루어지지만 설명변수가 2개 이상일 경우 변수를 배열구조 특히 행렬로 전환하는 것이 계산이나 모형을 이해하는데 도움이 됩니다. 간단한 예로 단순회귀모형에서 설명변수와 종속변수를 행렬로 전환하면 식 5와 같이 나타낼 수 있습니다.

\begin{align}\begin{bmatrix}y_1\\y_2\\ \vdots \\y_n \end{bmatrix}& =\begin{bmatrix}1& x_1\\1& x_2\\ \vdots& \vdots\\1 & x_n \end{bmatrix}\begin{bmatrix}b_0\\b_1 \end{bmatrix}\\ \Leftrightarrow &\; y=Xb \end{align}

(식 5)

식 6은 연립방정식을 행렬 시스템으로 나타낸 것으로 행렬 X가 역행렬이 존재하는 가역행렬일 경우 미지의 벡터 a를 결정할 수 있습니다.

\begin{align}\begin{aligned}x_{11}a_1+x_{12}a_2=y_1\\x_{21}a_1+x_{22}a_2=y_2 \end{aligned}&\Rightarrow \begin{bmatrix}x_{11}&x_{12}\\x_{21} &x_{22}\end{bmatrix} \begin{bmatrix}a_1\\a_2 \end{bmatrix} = \begin{bmatrix}y_1\\y_2 \end{bmatrix}\\ &\Rightarrow X\cdot a=y \\ &\Rightarrow a = X^{-1}y \end{align}

(식 6)

행과 열의 수가 같은 정방행렬 A의 역행렬은 A^-1과 같이 위첨자 -1을 사용하여 나타내며 AA^-1 = I가 성립합니다. I는 항등행렬입니다.

식 5는 기지(known)의 설명변수와 반응변수 그리고 미지(unknown)의 b₀, b₁으로 구성되는 선형시스템(연립방정식)으로 식 6과 같은 형태로 나타낼 수 있습니다. 그러므로 설명변수의 역행렬을 적용하여 회귀계수 b₀, b₁를 결정할 수 있지만 설명변수가 정방행렬이 아니므로 역행렬을 결정할 수 없습니다.

m×n 형태의 행렬로부터 정방행렬은 식 7과 같이 그 행렬의 전치행렬과의 행렬곱 연산에 의해 생성됩니다.

a × aT: (m×n) × (n×m) = (m×m)	(식 7)
또는
aT × a : (n×m) × (m×n) = (m×m)

식 7을 적용하여 식 8과 같이 벡터 b를 결정합니다.

\begin{align}&y =Xb\\ &\qquad \Downarrow\\& X^Ty=X^TXb\\ &\qquad \Downarrow\\&\begin{aligned}(X^TX)^{-1}X^Ty&=(X^TX)^{-1}X^TXb\\&=b\end{aligned} \\& \therefore\; b=(X^TX)^{-1}X^Ty \end{align}

(식 8)

예 1에 대한 회귀계수를 식 8의 행렬 연산을 적용하여 결정하면 다음과 같습니다.

xtx=np.dot(X_n0.T, X_n0)
xty=np.dot(X_n0.T, y_n)
b=np.linalg.solve(xtx, xty)
print(np.around(b,3))

[[-0.   ]
 [ 0.997]]