다운로드
작성자: admin 작성일시: 2016-06-03 18:30:17 조회수: 4731 다운로드: 237
카테고리: 머신 러닝 태그목록:

확률론적 선형회귀모형

OLS(Ordinary Least Square) 방법을 사용하면 데이터에 대한 확률론적인 가정없이도 최적의 가중치를 계산할 수 있다. 그러나 이 경우에는 계산한 가중치가 어느 정도의 신뢰도 또는 안정성을 가지는지 확인할 수 있는 방법이 없다. 이를 확인하고자 하는 시도 중의 하나가 부트스트래핑(bootstrapping) 방법이다.

부트스트래핑

부트스트래핑(bootstrapping)은 회귀 분석에 사용한 데이터가 달라진다면 회귀 분석의 결과는 어느 정도 영향을 받는지를 알기 위한 방법이다.

데이터가 확률 변수로부터 생성된 표본이거나 혹은 더 큰 모집단 중에서 선택한 표본이라고 가정한다면 회귀 분석의 결과는 분석에 사용한 표본에 의존적임을 알 수 있다. 만약 추가적인 다른 표본을 얻어서 다시 회귀 분석에 사용한다면 회귀 분석 결과 즉, 가중치 벡터의 값은 달라질 것이다.

그러나 현실적으로는 데이터를 추가적으로 얻기가 힘들기 때문에 부트스트래핑 방법에서는 기존의 데이터를 재표본화(re-sampling)하는 방법을 선택한다. 재표본화는 기존의 $N$개의 데이터에서 다시 $N$개의 데이터를 선택하되 중복 선택도 가능하게 한다(unordered resampling with replacement). 이 경우 이론적으로는 ${2N-1}\choose{N}$ 개의 새로운 표본 집단을 얻을 수 있다.

직접 부트스트래핑을 실시해 보자. 우선 100개의 가상 데이터를 생성한다. 데이터의 모수는 각각 $w_0=0$, $w_1=42.3855$이다.

In:
from sklearn.datasets import make_regression

X0, y, coef = make_regression(n_samples=100, n_features=1, noise=20,
                              coef=True, random_state=0)
coef
Out:
array(42.38550485581797)

이 데이터를 기반으로 회귀 분석을 실시한다.

In:
dfX0 = pd.DataFrame(X0, columns=["X1"])
dfX = sm.add_constant(dfX0)
dfy = pd.DataFrame(y, columns=["y"])

model = sm.OLS(dfy, dfX)
result = model.fit()

추정된 계수 정보는 params 속성에 저장되어 있다. 실제 모수와 어느정도 오차를 가짐을 알 수 있다.

In:
result.params
Out:
const    -1.628364
X1       42.853356
dtype: float64

다음으로 이 데이터에서 중복을 허락하여 N개의 데이터를 선택한 후 다시 회귀 분석을 한다. 여기에서는 1,000번만 반복해 본다.

In:
%%time
N = 1000
params_c = np.zeros(N)
params_x1 = np.zeros(N)
for i in range(N):
    idx = np.random.choice(len(dfy), len(dfy), replace=True)
    dfX2 = dfX.iloc[idx, :]
    dfy2 = dfy.iloc[idx]
    r = sm.OLS(dfy2, dfX2).fit()
    params_c[i] = r.params.const
    params_x1[i] = r.params.X1
CPU times: user 1.84 s, sys: 0 ns, total: 1.84 s
Wall time: 1.85 s

전체 가중치 집합을 히스토그램으로 나타내면 다음과 같다.

In:
plt.figure(figsize=(8, 3))

ax1 = plt.subplot(121)
sns.distplot(params_c, ax=ax1)
plt.axvline(params_c.mean(), c='r')
plt.title("const parameter $\hat{w}_0$")

ax2 = plt.subplot(122)
sns.distplot(params_x1, ax=ax2)
plt.axvline(params_x1.mean(), c='r')
plt.title("x1 parameter $\hat{w}_1$")
plt.show()

평균과 분산은 다음과 같다.

In:
params_c.mean(), params_c.std()
Out:
(-1.5273502720200236, 2.1313819398826008)
In:
params_x1.mean(), params_x1.std()
Out:
(42.747134522988972, 1.9661344907376426)

가중치 중 상수항의 경우 평균은 -1.6이지만 표분 편차가 2보다 크기 때문에 값이 0일 가능성을 배제할 수 없다.

이 결과를 StatsModels 의 회귀 분석 보고서와 비교하자.

In:
print(result.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.803
Model:                            OLS   Adj. R-squared:                  0.801
Method:                 Least Squares   F-statistic:                     400.3
Date:                Sat, 03 Mar 2018   Prob (F-statistic):           2.21e-36
Time:                        09:29:06   Log-Likelihood:                -448.09
No. Observations:                 100   AIC:                             900.2
Df Residuals:                      98   BIC:                             905.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -1.6284      2.163     -0.753      0.453      -5.920       2.663
X1            42.8534      2.142     20.008      0.000      38.603      47.104
==============================================================================
Omnibus:                        3.523   Durbin-Watson:                   1.984
Prob(Omnibus):                  0.172   Jarque-Bera (JB):                2.059
Skew:                          -0.073   Prob(JB):                        0.357
Kurtosis:                       2.312   Cond. No.                         1.06
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

보고서의 std err 항목을 보면 표준 편차의 경우 2.163 이고 마지막의 신뢰 구간(confidence interval)이 -5.920 ~ 2.663 임을 보이고 있다. 부트스트래핑으로 얻은 결과와 유사하다. 이 결과는 다음에 설명할 확률론적 가정에 의해 계산된 값이다.

확률론적 선형회귀모형

확률론적 선형회귀모형을 사용하면 부트스트래핑 없이도 회귀분석 결과의 통계적 특성을 구할 수 있다. 확률론적 선형회귀모형에서는 데이터가 확률 변수로부터 생성된 표본이라고 가정한다. 구체적인 가정은 다음과 같다.

가정 1. 오차의 분포에 대한 가정

1-1. 선형 정규 분포 가정

선형 회귀 분석의 기본 가정은 종속 변수 $y$가 독립 변수 $x$의 선형 조합으로 결정되는 기댓값과 고정된 분산 $\sigma^2$을 가지는 정규 분포라는 것이다.

$$ y \sim \mathcal{N}(w^Tx, \sigma^2) $$

y의 확률 밀도 함수는 다음처럼 쓸 수 있다. 이 식에서 모수 벡터 $\theta=(w, \sigma^2)$이다.

$$ p(y \mid x, \theta) = \mathcal{N}(y \mid w^Tx, \sigma^2 ) $$

이 관계식을 오차(disturbance) $ \epsilon_i$ 개념으로 변환하면 더 간단하게 표현할 수 있다.

$$ \epsilon_i = y - w^Tx $$

오차의 확률 밀도 함수는 다음과 같다.

$$ p(\epsilon \mid \theta) = \mathcal{N}(0, \sigma^2 ) $$


여기에서 주의할 점은

x, y 중 그 어느 것도 그 자체로 정규 분포일 필요는 없다

는 것이다.

$y$도 $x$에 대해 조건부로 정규 분포를 이루는 것이지 $y$ 자체가 무조건부로 정규분포는 아니다.

1-2. 외생성(Exogeneity) 가정

오차 $\epsilon$의 기댓값은 독립 변수 $x$에 상관없이 항상 0이라고 가정한다.

$$ \text{E}[\epsilon \mid x] = 0$$


외생성 가정으로부터와 오차 $\epsilon$의 무조건부 기댓값이 0임을 증명할 수 있다.

$$ \text{E}[\epsilon] = 0$$

그리고 같은 가정으로부터 오차 $\epsilon$와 독립 변수 $x$가 무상관 관계라는 것도 증명할 수 있다.

$$ \text{E}[\epsilon x] = 0$$

1-3. 조건부 독립 가정

$i$번째 샘플의 오차 $\epsilon_i$와 $j$번째 샘플의 오차 $\epsilon_j$의 공분산 값이 $x$와 상관없이 항상 0이라고 가정한다.

$$ \text{Cov}[\epsilon_i, \epsilon_j \mid x] = 0 \;\; (i,j=1,2,\ldots,N) $$

이는 $i$번째 샘플의 오차 $\epsilon_i$와 $j$번째 샘플의 오차 $\epsilon_j$는 서로 독립이라는 가정과 같다.

이 가정과 위의 다른 가정들로부터 다음을 증명할 수 있다.

$$ \text{E}[\epsilon_i \epsilon_j] = 0 \;\; (i,j=1,2,\ldots,N) $$

이를 공분산 행렬의 형태로 표현하면 다음과 같다.

$$\text{Cov}[\epsilon] = \text{E}[\epsilon^{} \epsilon^T] = \sigma^2 I $$

가정 2. 독립 변수에 대한 가정

독립 변수의 공분산 행렬은 full rank 이어야 한다.

$$ \text{rank}(X^TX) = D \;\; (X \in \mathbf{R}^{N \times D})$$

즉, 독립 변수에서 서로 독립인 성분이 독립 변수의 갯수만큼 존재해야 한다. 만약 공분산 행렬이 full rank가 아니면 $X^TX$의 역행렬이 존재하지 않으므로 올바른 모수 추정이 어렵다. 또는 공분산 행렬의 조건이 나쁜(ill-conditioned) 경우에도 올바른 모수 추정이 되지 않는다.

MLE를 사용한 선형 회귀 분석

확률론적 선형 회귀모형의 가정과 MLE(Maximum Likelihood Estimation)을 사용하여 가중치 벡터 $w$의 값을 구해보자.

Likelihood는 다음과 같다.

$$ \begin{eqnarray} p(y_{1:N} \,\big|\, x_{1:N}, \theta) &=& \prod_{i=1}^N \mathcal{N}(y_i \,\big|\, w^T x_i , \sigma^2) \\ &=& \prod_{i=1}^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(y_i-w^T x_i)^2}{2\sigma^2} \right\} \\ \end{eqnarray} $$

계산을 쉽게하기 위해 Log를 적용하여 log-likelihood를 구한다.

$$ \begin{eqnarray} \text{LL} &=& \log p(y_{1:N} \,\big|\, x_{1:N}, \theta) \\ &=& \log \prod_{i=1}^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(y_i-w^T x_i)^2}{2\sigma^2} \right\} \\ &=& -\dfrac{1}{2\sigma^2} \sum_{i=1}^N (y_i-w^T x_i)^2 - \dfrac{N}{2} \log{2\pi}{\sigma^2} \\ \end{eqnarray} $$

이를 행렬로 표시하면 다음과 같다.

$$ \text{LL} = -C_1 (y - Xw)^T(y-Xw) - C_0 = -C_1(w^TX^TXw -2 y^TXw + y^Ty) - C_0 $$$$ C_1 = -\dfrac{1}{2\sigma^2} $$$$ C_0 = \dfrac{N}{2} \log{2\pi}{\sigma^2} $$

이를 최적화하면 OLS와 동일한 결과를 얻을 수 있다.

$$ \dfrac{\partial}{\partial w} \text{LL} = -C_1 \left( 2X^TX \hat{w} - 2X^Ty \right) = 0 $$$$ \hat{w} = (X^TX)^{-1}X^T y $$

여기에서도 OLS의 정규 방정식과 같은 형태의 정규 방정식을 얻을 수 있다.

$$ X^TX \hat{w} - X^Ty = 0 $$$$ X^T(X \hat{w} - y) = X^T(\hat{y} - y) = X^Te = 0 $$

잔차의 분포

위의 확률론적 선형 회귀모형에 따르면 회귀분석에서 생기는 잔차 $e = y - \hat{w}^Tx$ 도 정규 분포를 따른다. 이는 다음과 같이 증명할 수 있다.

확률론적 선형회귀모형의 오차 $\epsilon$와 잔차 $e$는 다음과 같은 관계를 가진다.

$$ \hat{y} = X\hat{w} = X (X^TX)^{-1}X^T y = Hy $$

이 행렬 $H$은 Hat 행렬 혹은 프로젝션(projection) 행렬 또는 영향도(influence) 행렬이라고 부르는 대칭 행렬이다.

Hat 행렬은 다음과 같은 성질을 가진다.

$$ HH = H $$

Hat 행렬을 이용하면 잔차는 다음처럼 표현된다.

$$ e = y - \hat{y}= y - Hy = (I - H) y$$

$M = I - H$이라고 정의하자. 이 행렬 $M$은 잔차(residual) 행렬 혹은 rejection 행렬이라고 부르고 $H$와 마찬가지로 대칭 행렬이다. 또한 $H$와 마찬가지로 다음 성질을 가진다.

$$ MM = M $$

이 식을 이용하면 다음과 같다.

$$ e = My = M (Xw + \epsilon) = MXw + M\epsilon $$

그런데

$$ X^TM = X^T(I - H) = X^T - X^TH = X^T - X^TX (X^TX)^{-1}X^T = 0 $$

$M$은 대칭이므로

$$ M^TX = MX = 0 $$

잔차식에 이 값을 다시 대입하면

$$ e = MXw + M\epsilon = M\epsilon $$$$ e = M\epsilon $$

즉, 잔차 $e$는 오차 $\epsilon$의 선형 변환(linear transform)이다. 정규 분포의 선형 변환은 마찬가지로 정규 분포이므로 잔차도 정규 분포를 다른다.

In:
sp.stats.probplot(result.resid, plot=plt)
plt.show()

잔차의 정규성은 다음과 같이 정규성 검정을 통해 살펴볼 수도 있다.

In:
test = sms.omni_normtest(result.resid)
for xi in zip(['Chi^2', 'P-value'], test):
    print("%-12s: %6.3f" % xi)
Chi^2       :  3.523
P-value     :  0.172
In:
test = sms.jarque_bera(result.resid)
for xi in zip(['Jarque-Bera', 'P-value', 'Skew', 'Kurtosis'], test):
    print("%-12s: %6.3f" % xi)
Jarque-Bera :  2.059
P-value     :  0.357
Skew        : -0.073
Kurtosis    :  2.312

또한 오차의 기댓값이 $x$와 상관없이 0이므로 잔차의 기댓값도 $x$와 상관없이 0이어야 한다. 다음은 $x$값이 달라짐에 따라 잔차의 분포가 어떻게 바뀌는 지를 살펴보기 위한 것이다. $x$값이 달라져도 분포의 형태가 크게 바뀌지 않음을 알 수 있다.

In:
plt.plot(X0, result.resid, 'o')
plt.axhline(y=0, c='k')
plt.xlabel("X1")
plt.ylabel("Residual")
plt.show()

다음은 $0 < x < 1$인 구간과 $1 < x < 2$인 구간으로 나누어 분포의 모양을 살펴본 것이다. 둘 다 기댓값이 0에 가깝고 분산의 크기가 비슷함을 알 수 있다.

In:
sns.kdeplot(result.resid[((1 < X0) & (X0 < 2)).flatten()], label="1 < x < 2")
sns.kdeplot(result.resid[((0 < X0) & (X0 < 1)).flatten()], label="0 < x < 1")
plt.legend()
plt.show()

회귀 계수의 표준 오차

가중치 $\hat{w}$ 도 정규 분포 확률 변수인 $y$의 선형 변환이므로 정규 분포를 따른다.

$$ \begin{eqnarray} \hat{w} &=& (X^TX)^{-1} X^T y \\ &=& (X^TX)^{-1} X^T (X w + \epsilon) \\ &=& w + (X^TX)^{-1} X^T \epsilon \\ \end{eqnarray} $$

따라서 가중치의 기댓값은 다음과 같다.

$$ \begin{eqnarray} \text{E}[\hat{w}] &=& \text{E}[ w + (X^TX)^{-1} X^T \epsilon ] \\ &=& w + (X^TX)^{-1} X^T \text{E}[ \epsilon ] \\ &=& w \end{eqnarray} $$

가중치의 (공)분산을 계산하면 다음과 같다.

$$ \begin{eqnarray} \text{Cov}[\hat{w}] &=& E[(\hat{w} - w)(\hat{w} - w)^T] \\ &=& E[((X^TX)^{-1} X^T \epsilon)((X^TX)^{-1} X^T \epsilon)^T] \\ &=& E[(X^TX)^{-1} X^T \epsilon \epsilon^T X(X^TX)^{−1} ] \\ &=& (X^TX)^{-1} X^T E[\epsilon \epsilon^T] X(X^TX)^{−1} \\ &=& (X^TX)^{-1} X^T (\sigma^2 I) X(X^TX)^{−1} \\ &=& \sigma^2 (X^TX)^{-1} \end{eqnarray} $$

$\sigma^2$의 값은 알지 못하므로 다음과 같이 추정한다.

$$ s^2 = \dfrac{e^Te}{N-K} = \dfrac{RSS}{N-K} $$

여기에서 $N$은 표본 데이터의 수, $K$는 모수의 갯수이다. 바이어스 즉, y절편을 고려하는 선형 모형이라면 모수의 갯수는 입력 데이터의 차원의 수 $D$에 1을 더한 값이 된다.

$$ K = D + 1 $$

따라서 $\hat{w}$의 (공)분산의 추정값은 다음과 같다.

$$ \text{Est.Cov}[\hat{w}] = s^2(X^TX)^{-1}$$

이 공분산 행렬에서 우리가 관심을 가지는 값은 $w_i$의 분산을 뜻하는 대각 성분이다.

$$\text{Est.Var}[\hat{w}_i] = \left( \text{Est.Cov}[\hat{w}] \right)_{ii} \;\; (i=0, \ldots, K-1)$$

이 값에서 구한 표준 편차를 회귀 계수의 표준 오차(Standard Error of Regression Coefficient)라고 한다.

$$ {se_i} = \sqrt{s^2 \big((X^TX)^{-1}\big)_{ii}} \;\; (i=0, \ldots, K-1)$$

실제 가중치 계수 $w_i$와 우리가 추정한 가중치 계수 $\hat{w}_i$의 차이, 즉 모수 오차는 기댓값이 0이고 표준 오차를 분산으로 가지는 스튜던트-t 분포를 따른다.

단일 계수 t-검정 (Single Coefficient t-test)

$\hat{w}$을 위에서 구한 표준 오차로 나눈 값은 자유도가 $N-K$인 표준 스튜던트 t 분포를 따른다.

$$ \dfrac{\hat{w}_i - w_i}{se_i} \sim t_{N-K} \;\; (i=0, \ldots, K-1)$$

이를 검정 통계량(test statistics)로 사용하면 특정 회귀 계수 $w_i$가 0 인지 아닌지에 대해 조사할 수 있다.

$$ H_0 : \;\; w_i = 0 \;\; (i=0, \ldots, K-1)$$

만약 유의확률이 0에 가깝게 나온다면 위의 귀무가설은 기각이므로 $w_i$ 값이 0일 가능성은 적다. 하지만 유의확률이 유의수준을 넘는 큰 값이 나온다면 반대로 $w_i$ 값이 0일 가능성이 크다. 즉, 해당 독립변수는 종속변수와 아무런 상관성이 없을 가능성이 있다는 뜻이다.

StatsModels summary 메서드가 출력하는 회귀분석 보고서에서 std err로 표시된 열이 모형계수의 표준오차, t로 표시된 열이 단일 계수 t-검정의 검정 통계량, 그리고 P>|t|로 표시된 열이 유의확률을 뜻한다.

아래의 보고서에서는 첫번째 모형계수는 실제 값이 0일 가능성 즉, 상수항을 가지지 않을 가능성이 높다고 할 수 있다.

In:
print(result.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.803
Model:                            OLS   Adj. R-squared:                  0.801
Method:                 Least Squares   F-statistic:                     400.3
Date:                Sat, 03 Mar 2018   Prob (F-statistic):           2.21e-36
Time:                        09:29:10   Log-Likelihood:                -448.09
No. Observations:                 100   AIC:                             900.2
Df Residuals:                      98   BIC:                             905.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -1.6284      2.163     -0.753      0.453      -5.920       2.663
X1            42.8534      2.142     20.008      0.000      38.603      47.104
==============================================================================
Omnibus:                        3.523   Durbin-Watson:                   1.984
Prob(Omnibus):                  0.172   Jarque-Bera (JB):                2.059
Skew:                          -0.073   Prob(JB):                        0.357
Kurtosis:                       2.312   Cond. No.                         1.06
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

RegressionResults 클래스 객체는 t test 를 위한 t_test 메서드를 제공한다. 이 메서드를 사용하면 계수 값이 0이 아닌 경우도 테스트할 수 있다. 예를 들어 위 예제에서 귀무 가설을 다음처럼 놓고

$$ H_0: w_1 = 40 $$

x1에 대한 계수값 $w_1$을 40으로 해도 되는지 테스트할 수도 있다.

In:
print(result.t_test("X1 = 40"))
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            42.8534      2.142      1.332      0.186      38.603      47.104
==============================================================================

이 검정 결과에 따르면 x1에 대한 계수를 40으로 한다고 해도 문제가 없음을 알 수 있다.

Loss-of-Fit 검정 (회귀 분석 F-검정)

개별 개수가 아닌 전체 회귀 계수가 모두 의미가 있는지 확인하는 경우에는 다음과 같은 귀무 가설을 생각할 수 있다.

$$ H_0 : w_0 = w_1 = \cdots = w_{K-1} = 0 $$

이는 전체 독립 변수 중 어느 것도 의미를 가진 것이 없다는 뜻이다. 대부분의 경우, 이 귀무가설은 기각된다. 다만 유의 확률이 얼마나 작은가에 따라서 기각되는 정도가 달라진다. 유의 확률이 작으면 작을수록 더 강력하게 기각된 것이므로 더 의미가 있는 모형이라고 할 수 있다. 따라서 여러 모형의 유의 확률을 비교하여 어느 모형이 더 성능이 좋은가를 비교할 때 이 유의 확률을 사용한다. 이러한 검정을 Loss-of-Fit 검정 또는 회귀 분석 F-검정(regression F-test)이라고 한다.

위 보고서에서 F-statistic라고 표시된 400.3이라는 값이 회귀분석 F-검정의 검정통계량이고 Prob (F-statistic)로 표시된 2.21e-36라는 값이 유의확률이다.

질문/덧글

확률론적 회귀모형 가정 중 y가 x에 조건부 정규분포? samk*** 2017년 11월 18일 2:57 오전

확률론적 회귀모형 가정1에서
강조해주신 y가 x에 조건부 정규분포이지 무조건적인 정규분포가 아니다 라는 내용이 어렵습니다.

바꾸어 말하면 어떤 x에서는 y가 정규분포가 아니게 되어 회귀모형이 성립안될 수도 있다는 가정인건가요?

확률론적 선형회귀모형 하위단원에 있는 가정(Assumptions)에 대한 문의 samk*** 2017년 11월 18일 3:11 오전

확률론적 선형회귀모형 하위단원에 있는 exogeneity, 독립 변수의 공분산 행렬은 full rank에 대한 가정(Assumptions)들이

위키에서는 OLS하위 단원에 나와있어서 질문드립니다.

https://en.wikipedia.org/wiki/Ordinary_least_squares

혹시 저 가정들은 모든 선형회귀모형의 가정인가요?(결정론적 선형회귀모형, 확률론적 선형회귀모형에 관계없이)

부트 스트랩 질문입니다. wjy5*** 2018년 3월 8일 1:21 오후

부트 스트랩 샘플 추출 시에 2^N이라고 하셨는데 N^N아닌가요?

답변: 부트 스트랩 질문입니다. 관리자 2018년 3월 9일 9:36 오전

정확한 수식은 ${2N-1}\choose{N}$ 입니다. [`https://www.probabilitycourse.com/chapter2/2_1_4_unordered_with_replacement.php`](https://www.probabilitycourse.com/chapter2/2_1_4_unordered_with_replacement.php)
를 참조하세요. 수정하겠습니다. 지적 감사합니다.