"Quantile Regression - Theory and Applications"(Cristina Davino 외 2명) 2장을 공부했다.
저번 포스팅에서는 Quantile Regression(QR)의 기본적인 내용을 소개하고 파이썬 예제를 살펴보았다면 이번 포스팅에서는 먼저 QR 파라미터를 어떻게 추정하는지 알아보고 이를 파이썬을 이용하여 구현해보고자 한다(밑바닥부터 구현하는 건 아니다 ㅎㅎ). 그리고 여러 가지 실험을 통하여 QR을 왜 사용하는지에 대해서 고찰해보고자 한다.
1. QR 파라미터 추정법
1. 추정 방법
먼저 데이터 $(y_i, \tilde{x}_i^t), i=1, \ldots, n$ 가 있다고 하자. 여기서 $\tilde{x}_i^t = (1, x_{i1}, \ldots, x_{ip})$이다. QR 파라미터는 아래의 식을 최소화하는 값으로 추정하게 된다.
$$\DeclareMathOperator*{\argminB}{argmin} \argminB_{\beta}\sum_{i=1}^n \rho_\theta(y_i-\tilde{x}_i^t\beta)\tag{1}$$
여기서 $\rho_\theta= \theta [y]_+ + (1-\theta) [-y]_+$이고 $[y]_+=\max \{0, y\}$ 이다. 식 (1)을 최소화하는 $\beta$는 Linear Programming을 이용하여 구할 수 있다. 식 (1)을 최소화하는 문제를 Linear Programming으로 바꿔보자.
먼저 $X$를 $i$ 번째 행이 $\tilde{x}_i^t$인 행렬, $y = (y_1, \ldots, y_n)$이라 할 때 $s_1 = [y-X\beta]_+, s_2 = [X\beta-y]_+$라 하자. 여기서 $[\cdot]_+$는 스칼라 연산이 아니고 벡터 연산이다. 즉, 원소별로 연산한 것이다. 그렇다면 다음을 알 수 있다.
$$y - X\beta = s_1-s_2$$
그리고 $\beta = [\beta]_+ - [-\beta]_+$라 하자. 그렇다면 식 (1)을 최소화하는 문제는 다음과 같이 Linear Programming 문제가 된다.
$$\DeclareMathOperator*{\argminB}{argmin} \argminB_{\beta} \{ \theta 1^ts_1+(1-\theta)1^ts_2 | y=X\beta + s_1 - s_2, s_j \in \mathbf{R}_+^{n} , j=1, 2\}\tag{2}$$
이를 좀 더 바꿔보자. $I$는 $n\times n$ 단위행렬이라 할 때
$B = [X -X I -I]$, $\psi = ([\beta]_+^t, [-\beta]_+^t, s_1^t, s_2^t)^t$, $d = (0_p^t, 0_p^t, \theta 1_n^t, (1-\theta)1_n^t)^t$라 한다면 식 (2)는 다음과 같이 표현할 수 있다.
$$\begin{align} & \text{minimize}_\beta \;\; d^t\psi \\ & \text{subject to } B\psi = y, \psi \geq 0 \end{align}$$
2. 파이썬 구현
이제 실제 데이터를 이용하여 Linear Programming으로 파라미터를 추정해보자. 데이터를 다운받아준다.
필요한 모듈을 임포트하고 데이터를 불러온다.
import pandas as pd
import numpy as np
import statsmodels.api as sm
from scipy.optimize import linprog
car_df = pd.read_csv('../dataset/Cars93.csv')
Price와 Origin의 관계를 보기 위해 필요한 칼럼을 추출하고 Origing을 더미화시켜준다.
car_df = car_df[['Price','Origin']]
car_df = pd.get_dummies(car_df,columns=['Origin'],drop_first=True)
Linear Programming은 scipy 모듈에서 linprog를 이용하여 풀 수 있다. 아래는 이를 이용하여 5개의 분위수에 대한 파라미터를 추정하는 코드이다. 특히 12~14번째 줄을 주목하자. 이는 Linear Programming을 정의하고 linprog의 인자로 넣어주는 아주 중요한 부분이다. linprog의 대한 설명은 여기를 참고하기 바란다.
y = car_df['Price'].values
X = car_df['Origin_non-USA']
X = sm.add_constant(X)
X = X.values
n = X.shape[0] ## the number of data
p = X.shape[1] ## the number of parameters
qs = [0.1, 0.25, 0.5, 0.75, 0.9]
for q in qs:
A_eq = np.hstack([X, -X, np.identity(n), -np.identity(n)])
b_eq = y
c = np.array([0]*(2*p) + [q]*n + [1-q]*n)
res = linprog(c, A_eq=A_eq, b_eq=b_eq,method='highs')
x_plus = res.x[:p]
x_minus = res.x[p:2*p]
x = x_plus-x_minus
print(f'q : {q}, intercept : {x[0]}, slope : {round(x[1],3)}')
위 코드를 실행하면 다음과 같은 결과가 나오며 이는 저번 포스팅에서 statsmodel을 이용해서 구했던 파라미터 값과 동일하다.
2. QR을 어떻게 그리고 왜 사용하는가?
여기서는 실험을 통하여 QR을 사용해야 하는 이유에 대해서 알아보고자 한다.
먼저 다음과 같은 10개의 모델을 고려하며 각 모델마다 10,000개의 샘플을 취한다.
$$\begin{align} & y_j = 1+2x + e_j, j=1,2,\ldots, 10, x \sim N(10,1) \\& e_1 \sim N(0, 1) \\ & e_2 \sim LN(0, 0.25) \\ & e_3 \sim LN(0, 0.5) \\ & e_4 \sim LN(0.1.25) \\ & e_5 = -e_4 \\& e_6 = (1+x)e_1 \\ & e_{AR(1)[\rho]} \rightarrow e_i = \rho e_{i-1} + a_i, e_0 = 0, a_i \sim N(0,1) \\& e_7 = e_{AR(1)[-0.2]} \\& e_8 = e_{AR(1)[0.2]} \\& e_9 = e_{AR(1)[-0.5]} \\& e_{10} = e_{AR(1)[0.5]} \end{align}$$
아래 코드는 위 모델을 구현한 것이다.
np.random.seed(10)
n_sample = 10000 ## 샘플 개수
x = np.random.normal(10,1,n_sample)
error1 = np.random.normal(0,1,n_sample)
error2 = np.random.lognormal(0,0.25,n_sample)
error3 = np.random.lognormal(0,0.5,n_sample)
error4 = np.random.lognormal(0,1.25,n_sample)
error5 = -error4
error6 = (1+x)*error1
a_list = np.random.normal(0,1,n_sample)
rhos = [-0.2,0.2,-0.5,0.5]
errors = []
for rho in rhos:
error = []
e_cur = 0
for a in a_list:
e_next = rho*e_cur + a
error.append(e_next)
e_cur = e_next
errors.append(error)
y1 = 1+2*x+error1
y2 = 1+2*x+error2
y3 = 1+2*x+error3
y4 = 1+2*x+error4
y5 = 1+2*x+error5
y6 = 1+2*x+error6
y7 = 1+2*x+errors[0]
y8 = 1+2*x+errors[1]
y9 = 1+2*x+errors[2]
y10 = 1+2*x+errors[3]1. QR for Homogeneous and Heterogeneous Models
먼저 모델 1은 정규분포를 따르는 등분산(가장 이상적인) 데이터이고 모델 6은 이분산 데이터이다. 여기서는 등분산 데이터와 이분산 데이터에 대하여 QR이 어떻게 추정되는지 살펴보자.
아래 코드는 각 데이터의 6개 분위수에 대한 QR 직선을 추정하고 시각화한 것이다. 최소 제곱 회귀 직선도 포함시켰다. 여기서는 위에서 구현한 방법이 아닌 statsmodels에서 제공하는 QuantReg를 사용했다.
X = np.expand_dims(x,axis=1)
X = sm.add_constant(X)
ys = [y1, y6]
qs = [0.05,0.1,0.25,0.75,0.9,0.95]
fig = plt.figure(figsize=(6,12))
fig.set_facecolor('white')
for j in range(len(ys)):
ax_idx = 210+j+1
ax = fig.add_subplot(ax_idx)
y = ys[j]
ax.scatter(x,y,color='y')
for q in qs:
model = sm.QuantReg(y,X).fit(q=q)
params = model.params
ax.plot(x,params[0]+params[1]*x)
ols = sm.OLS(y,X).fit()
ax.plot(x,ols.params[0]+ols.params[1]*x,color='k',linestyle='--')
plt.show()
위 2개의 그림에서 상단은 등분산 데이터에 대한 QR 직선을 그린 것이고 하단은 이분산 데이터에 대한 QR 직선을 그린 것이다. 먼저 등분산 데이터에 대해서 QR 직선 기울기 추정량은 대략 1.98~2.0 사이로 거의 동일했으며 값이 조금씩 차이가 나는 이유는 Sampling Variance 때문이다. 실제 True 모델에서 오차항이 표준 정규분포를 따르므로 QR 직선은 기존 최소 제곱 추정법으로 얻은 직선을 평행 이동하여 얻을 수 있다.
아래 그림을 보면 QR 직선의 기울기가 각 분위수 마다 다른 것을 알 수 있으며 데이터의 분포가 커지는 방향으로 직선의 방향이 넓어지는 것을 알 수 있다. 기존 최소 제곱 직선은 오로지 반응 변수의 조건부 기대값만을 추정하므로 다시 말하면 하나의 직선만 있으므로 이분산 데이터를 제대로 설명할 수 없다. 하지만 QR은 이분산 데이터에서 적절한 분위수를 가지고 이분산 데이터의 특징을 잘 나타낼 수 있다.
이번엔 모델 1과 모델 2, 3을 비교해보고자 한다. 이는 오차항이 등분산 정규 분포를 따를 때, 그리고 등분산 이지만 정규 분포가 아닌 데이터에 대해서 QR을 비교해보고자함이다. 아래 코드를 실행해보자.
from scipy.stats import norm
import seaborn as sns
X = np.expand_dims(x,axis=1)
X = sm.add_constant(X)
ys = [y1, y2, y3]
qs = [0.25,0.5,0.9]
colors = sns.color_palette('hls',len(qs))
fig = plt.figure(figsize=(6,12))
fig.set_facecolor('white')
for j in range(len(ys)):
ax_idx = 310+j+1
ax = fig.add_subplot(ax_idx)
y = ys[j]
ax.scatter(x,y,color='y')
ols = sm.OLS(y,X).fit()
fitted = ols.fittedvalues
ax.plot(x,fitted,color='k')
for i, q in enumerate(qs):
model = sm.QuantReg(y,X).fit(q=q)
params = model.params
ax.plot(x,params[0]+params[1]*x, color = colors[i], label=str(q)+'-quantile')
if q != 0.5:
ax.plot(x,fitted+norm.ppf(q),color='k', linestyle=(0,(2,4)))
ax.legend()
plt.show()
먼저 첫 번째 그림은 오차항이 등분산 정규 분포를 따르는 데이터이다. 이때 0.25 분위수에 대한 QR 직선은Ordinary Least Square(OLS) 직선을 -0.67(표준 정규 분포에 4분위수에 해당하는 값) 만큼 평행이동하면 얻을 수 있다. 마찬가지로 +1.28(표준 정규 분포에 0.9 분위수에 해당하는 값) 만큼 평행이동하면 0.9 분위수에 대응하는 QR 직선을 얻을 수 있다(물론 Sampling Variance에 따른 오차가 약간 있다).
그러나 실제 오차항의 분포가 정규분포가 아닌 경우 더이상 OLS를 평행이동하여 QR 직선을 얻을 수 없다. 그에 반해 QR 직선은 오차항의 분포의 상관없이 해당 분위수에 대응하는 QR 직선을 잘 추정하는 것을 알 수 있다. 즉, QR은 오차에 대한 사전 정보 없이 조건부 반응 변수의 분포를 OLS 보다 잘 추정할 수 있다는 것이다.
지금까지 비교 실험을 통해 하나 더 알 수 있는 것은 QR은 오차항의 이분산성과 정규성에 Robust하다는 것이다. 하지만 QR 추정량은 설명 변수의 이상치에 민감하다고 알려져 있다.
2. QR prediction intervals
QR 추정량을 이용하여 예측 구간을 계산할 수 있다. 이 방법으로 얻은 예측 구간은 정규 분포 가정에 Robust하다는 것이 알려져 있다. $\hat{q}(\theta,x)$를 $x$가 주어졌을 때 조건부 $\theta$ 분위수 추정값이라고 하자. 이때 구간 $[\hat{q}(\theta_1,x), \hat{q}(\theta_2,x)]$는 $x$가 주어졌을 때 반응 변수의 $(\theta_2-\theta_1)$% 예측 구간이 된다. 즉, QR 직선만을 이용하여 간단하게 예측 구간을 구할 수 있다. 간단하게 구한 예측 구간이지만 굉장히 안정적(Stable)이라고 한다.
실제로 그런지 비교실험을 진행해보고자 한다. $\alpha$를 Nominal Coverage라 할 때 실제 예측값을 포함하는 예측 구간의 비율이 Nominal Coverage와 같은지 살펴볼 것이다.
먼저 각 모델에서 오차항의 실제 분포를 알고 있으니(실제 사례에서는 그럴 일이 없겠지만 ㅎㅎ) 이론적인 예측구간을 구할 수 있다. 이때 모델 1~6은 실제 분포를 알고 있으니 이를 통하여 0.1 분위수와 0.9 분위수를 구할 수 있지만 모델 7~10은 직접 계산하기 어려우므로 반복 샘플링을 이용하여 0.1 분위수와 0.9 분위수를 계산하였다. 코드와 그 결과는 다음과 같다.
np.random.seed(100)
new_x = [8,9,10,11,12]
alpha = 0.2
models = 10
z_coef = norm.ppf(1-0.5*alpha)
replication = 100000
header = [f'x = {j}' for j in new_x]
header = ' '.join(header)
res = dict()
print(header)
for m in range(models):
text = 'model'+str(m+1)
res_val = []
for nx in new_x:
pred_val = 1+2*nx
if m == 0:
lower = pred_val - z_coef
upper = pred_val + z_coef
elif m == 1:
lower = lognorm.ppf(q=0.1,s=0.25,loc=pred_val)
upper = lognorm.ppf(q=0.9,s=0.25,loc=pred_val)
elif m == 2:
lower = lognorm.ppf(q=0.1,s=0.5,loc=pred_val)
upper = lognorm.ppf(q=0.9,s=0.5,loc=pred_val)
elif m == 3:
lower = lognorm.ppf(q=0.1,s=1.25,loc=pred_val)
upper = lognorm.ppf(q=0.9,s=1.25,loc=pred_val)
elif m == 4:
lower = -lognorm.ppf(q=0.9,s=1.25,loc=-pred_val)
upper = -lognorm.ppf(q=0.1,s=1.25,loc=-pred_val)
elif m == 5:
lower = pred_val + norm.ppf(0.1,scale=1+nx)
upper = pred_val + norm.ppf(0.9,scale=1+nx)
else:
temp_lower_list = []
temp_upper_list = []
for j in tqdm(range(replication),total=replication):
temp_data = 1+2*nx+generate_error(m,nx,1000)
temp_lower = np.quantile(temp_data,0.1)
temp_upper = np.quantile(temp_data,0.9)
temp_lower_list.append(temp_lower)
temp_upper_list.append(temp_upper)
lower = np.mean(temp_lower_list)
upper = np.mean(temp_upper_list)
text += f' [{round(lower,2)} {round(upper,2)}] '
res_val.append([round(lower,2), round(upper,2)])
res['model'+str(m+1)] = res_val
print(text)
예를 들어 모델1에서$x=8$인 경우 예측 구간은 (15,72, 18.28)이 된다. 이번엔 예측 구간을 추정해보자. 아래 코드는 예측 구간을 추정한 코드이다. 이때 OLS의 예측 구간도 포함시켰다.
np.random.seed(100)
n_sample = 1000
x = np.random.normal(10,1,n_sample)
replication = 1000
models = 10
theta1 = 0.1
theta2 = 0.9
alpha = 0.2
z_coef = norm.ppf(1-0.5*alpha)
new_x = [8,9,10,11,12]
header = [f'x = {j}' for j in new_x]
header = ' '.join(header)
print(header)
for m in range(models):
qr_text = 'model'+str(m+1)+' QR '
ls_text = 'model'+str(m+1)+' LS '
error = generate_error(m,x,n_sample)
y = 1+2*x+error
X = np.expand_dims(x,axis=1)
X = sm.add_constant(X)
ols = sm.OLS(y,X).fit()
model1 = sm.QuantReg(y,X).fit(q=theta1)
params1 = model1.params
model2 = sm.QuantReg(y,X).fit(q=theta2)
params2 = model2.params
for i, nx in enumerate(new_x):
lower = round(params1[0]+params1[1]*nx,2)
upper = round(params2[0]+params2[1]*nx,2)
qr_text += str([lower,upper]) + ' '
std = np.sqrt(ols.mse_resid*(1+1/n_sample+np.square(nx-np.mean(x))/np.sum(np.square(x-np.mean(x)))))
pred_val = ols.params[0] + ols.params[1]*nx
lp = pred_val-z_coef*std
up = pred_val+z_coef*std
ls_text += str([round(lp,2),round(up,2)])+' '
print(ls_text)
print(qr_text)
OLS와 QR의 예측 구간을 비교해보면 대체로 구간 길이는 비슷하다. 이 책에 저자는 이러한 예측 구간 추정을 반복하여 실제 이론적인 예측 구간을 포함하는 비율을 계산하여 그 비율이 0.8에 가까운지 실험했다고 한다. 나도 똑같이 구현하려고 했으나 결과가 책과 많이 달랐고(아마 내가 잘못한 것 같다) 반복 횟수가 몇 번인지 그리고 정확한 실험 내용이 소개되어 있지 않아 이 내용은 책에 있는 것을 캡처했다.
먼저 OLS 예측 구간의 Coverage rate를 보면(LS 행) 모델 1에서는 좋고 나머지 모델에서는 별로 좋지 못하다. 모델 4~5는 OLS 예측 구간의 성능이 매우 저조함을 알 수 있다. 반면 QR(QR 행) 예측 구간은 어떤 모델이든 간에 꾸준히 0.8에 가까운 Coverage rate를 보인다는 것을 알 수 있다. 즉, QR의 예측 구간은 오차항의 분포에 상관없이 매우 안정적으로 주어진 유의 수준만큼의 Coverage rate를 달성하는 것을 알 수 있다.
아래 표는 OLS와 QR을 이용하여 예측 구간 평균 길이를 계산한 것이다.
위 테이블에서 알 수 있듯이 일반적으로 QR의 예측 구간 평균 길이가 더 짧다는 것을 알 수 있다. 이는 정보력 측면에서 QR 예측 구간이 더 좋다는 뜻이 된다. 즉, QR 예측 구간은 OLS에 비해 더 안정적이고 정보력 측면에서 유리하기 때문에 QR 예측 구간을 사용하자!
3. 추정해야 할 Quantile 개수와 Quantile 간 거리
지금까지 QR 추정량에서 사용한 Quantile은 기껏해야 6개(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)였다. 그렇다면 궁금증이 생긴다. 추정해야할 Quantile는 몇 개로 해야 할까? 개수가 정해졌다고 하면 각 Quantile 간에 거리를 동등하게 해야할까? 다르게 해야할까?
이 책에 저자는 데이터 개수가 많아질수록 변수의 개수가 많아질수록 추정할 수 있는 유니크한 Quantile 개수는 많아지는 것을 실험적으로 보였다. 즉, 데이터 개수가 많고 변수의 개수가 많을수록 Quantile의 개수도 늘려서 추정해야 한다는 것이다. 또한 오차항의 분포 특성에 따라서도 추정된 Quantile의 개수가 달라진다고 한다. 하지만 구체적으로 몇 개의 Quantile을 추정해야 하는지에 대해서는 언급이 없었다.
그리고 인접한 Quantile 간 거리에 대해서는 데이터가 많아질수록 무시해도 된다고 한다. 즉, Quantile 간 거리를 동등하게 할지 아닐지는 신경 안 써도 된다. 물론 여기서도 단순히 데이터가 많으면 신경 안써도 된다고 했지 정확히 데이터가 몇개 이상일때 신경 안써도 되는지에 대해서는 나와있지 않다.
4. 요약
1. QR은 오차에 대한 사전 정보 없이, 다시 말해 이분산이든 등분산이든, 정규분포를 따르던 말던 조건부 반응 변수의 분포를 (OLS 보다) 잘 추정할 수 있다는 것이다.
2. QR 예측 구간은 OLS에 비해 더 안정적이고 정보력 측면에서 유리하다.
(즉, 1~2번의 이유로 QR을 써야 한다는 것이다.)
3. 하지만 QR 추정량은 설명 변수의 이상치에 민감하다고 알려져 있다.
'통계 > 머신러닝' 카테고리의 다른 글
| 9. 의사결정나무(Decision Tree) 에 대해서 알아보자 with Python (3) | 2021.06.10 |
|---|---|
| 8. 연관 규칙 분석(Association Rule Analysis) with Python (0) | 2021.05.23 |
| [Quantile Regression] 1. A visual introduction to quantile regression (0) | 2021.04.24 |
| 7. 이상치 탐지(Outlier Detection) - 통계적 검정과 여러가지 판별법 with Python (2) | 2021.02.17 |
| 6. Least Absolute Deviation Regression에 대해서 알아보자 with Python (0) | 2021.01.26 |
댓글