[data analysis] 최대우도 추정: AIC와 BIC

최대우도 추정: AIC와 BIC

AIC, BIC 등 확률적 모델 선택을 위한 평가 기준은 최대우도추정(maximum kikelihood estimation)을 근거로 작동합니다. 최대우도 추정은 관찰치를 근거로 반응변수(라벨)를 추정하는데 우도(likelihood, 가능도)를 최대로 하기 위한 모델의 매개변수(paramenter)를 발견하는 것입니다. 예를 들어 사건의 발생(1) 확률이 μ, 발생하지 않을(0) 확률이 1-μ라 하고 샘플 (0, 1, 0, 0, 1, 0)이 이항분포를 따른다고 할 경우 최대우도 추정(L(μ))는 식 1과 같이 계산됩니다.

\begin{align}P(x=1)&=\mu, \quad P(x=0)=1-\mu\\ \tag{식 1}L(\mu)&=P(x=0)\cdot P(x=1) \cdot P(x=0)\cdot P(x=0)\cdot P(x=1)\cdot P(x=0)\\ &=(1-\mu)\cdot \mu \cdot (1-\mu)\cdot(1-\mu)\cdot\mu\cdot(1-\mu)\cdot\\ &=(1-\mu)^4\mu^2 \end{align}

수학적 편의를 위해 식 1은 2와 같이 양변에 로그화를 진행하고 최대값을 산출하기 위해 μ에 대한 미분의 극값을 계산합니다. 즉, $\frac{\partial \log(\mu)}{\partial \mu} = 0$

\begin{align}&\begin{aligned}\log(L(\mu)) &= \log((1-\mu)^4\mu^2)\\ &=4\log(1-\mu)+2\log(\mu)\end{aligned}\\ &\tag{식 2} \frac{\partial \log(\mu)}{\partial \mu} = 0\\ & \Rightarrow 4\frac{1}{1-\mu}(-1)+2\frac{1}{\mu}=0\\ & \Rightarrow -4\mu+2-2\mu=0\\ & \Rightarrow \mu=\frac{1}{3} \end{align}

위에서 계산한 극점이 극소값인지 극대값인지를 결정해야 합니다. 이것은 식 3과 같이 2에서 계산한 미분을 다시 미분한 값으로 알 수 있습니다. 즉, 2차 미분한 값이 음의 값일 경우 극대값이 됩니다.(1차와 2차 미분계수의 극값참조)

\begin{align}\frac{\partial^2 \log(L(\mu))}{\partial^2 \mu}&=-4\frac{1}{(1-\mu)^2}-\frac{2}{\mu^2}\\\tag{식 3} &=-4\frac{1}{\left(1-\frac{1}{3}\right)^2}-\frac{2}{\left(\frac{1}{3}\right)^2}\\&=-27 \end{align}

위 과정으로 $\mu =\frac{1}{3}$에서 L(μ)이 극대값이 됩니다. 그러므로 극대값은 다음과 같습니다.(식 4)

\begin{align}\tag{식 4}L(\mu)&=(1-\mu)^4\mu^2\\&=\left(1-\frac{1}{3}\right)^4\left(\frac{1}{3}\right)^2\\&=0.021948\end{align}

위 최대우도추정을 로그우도 함수로 대체한다면 식 5와 같습니다.

\begin{align}\tag{식 5}\ln(L(\mu))&=\ln(0.021948)\\&=-3.1819 \end{align}

AIC

AIC는 정보이론과 빈도를 근거한 추론을 토대로 식 AIC_1과 같이 계산되며 모델 선택 도구로 사용됩니다.

$$\tag{식 6} \text{AIC} =-2 \cdot \ln(L)+2k $$

위 식에서 L은 우도(likelihood), k는 변수의 수를 나타냅니다. 이식에서 우도가 증가할수록 AIC가 감소됩니다. 또한 동일한 우도에서 변수의 수가 작은 경우에 비해 큰 경우 AIC는 증가하게 됩니다. 즉, 모델의 복잡도의 증가는 AIC에 의한 모델선택에 불리하게 됩니다. 결과적으로 이 통계량은 모델 성능과 모델 복잡도를 모두 고려합니다.

연속형 변수의 경우 우도는 mse로 판단할 수 있습니다. 즉, 우도가 증가할수록 mse는 감소할 것입니다. 식 6는 7과 같이 수정된 형태로 계산될 수 있습니다.

$$\tag{식 7} \text{AIC} =n\cdot \ln(MSE)+2k $$

식 7에서 n은 샘플의 크기입니다.

BIC

AIC와 유사하지만 모델 복잡도에 샘플의 수를 고려하는 차이를 가집니다. 즉, 모델 선택에 있어 모델의 복잡도의 영향을 증가시킨것으로 변수 선택에 따른 모델의 영향이 클 경우 BIC 통계량이 우선적으로 사용합니다.

$$\tag{식 8} \text{BIC} =-2\cdot \ln(L)+\ln(n)\cdot k $$

AIC와 같이 연속형 자료일 경우 MSE를 식 BIC_2와 같이 수정하여 적용합니다.

$$\tag{식 9} \text{BIC} =n\cdot \ln(MSE)+\ln(n) \cdot k$$

예)

다음 데이터에서 변수의 수에 따른 AIC와 BIC를 계산합니다. 데이터는 일일 kospi 지수 자료로서 "Close", "Open", "High", "Low", "Volume", "Change" 변수로 구성되어 있습니다. 각 변수에 대해 5일, 20일, 75일 이동평균을 계산해 변수로 첨가하였으며 "Close"를 반응변수로 하는 회귀분석을 실시하였습니다. 반응변수는 독립변수에 비해 1일 후 데이터 입니다.

	Open	High	Low	Close	Volume	5	20	75
Date
2021-04-21	3214.239990	3214.239990	3165.489990	3171.659912	0.363293	3196.830029	3127.433008	3086.599066
2021-04-22	3174.520020	3196.969971	3174.520020	3177.520020	-0.008149	3193.468018	3135.892505	3089.706667
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

import numpy as np
import pandas as pd
import yfinance as yf
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

st=pd.Timestamp(2021, 1,2)
et=pd.Timestamp(2024, 9, 30)
da=yf.download('^KS11', st, et)
da=da.drop(columns="Adj Close")
da["Volume"]=da["Volume"].pct_change()
for i in [5, 20, 75]:
    da[str(i)]=da["Close"].rolling(i).mean()
da1=da.dropna()

위 자료로 부터 독립변수와 반응변수를 구분하고 변수간 값의 scale과 편차의 정도가 크기 때문에 안정적인 모델의 구축을 위해서는 표준화를 진행합니다. MinMaxScler()를 적용합니다. 또한 훈련과 검증 세트를 구분합니다.

ind=da1.drop(columns="Close")[:-1]
de=da1["Close"][1:]
final=da1.drop(columns="Close").tail(1)
scal=MinMaxScaler().fit(ind)
X=scal.transform(ind)
finalN=scal.transform(final)
Xtr, Xte, ytr, yte=train_test_split(X,de, test_size=0.3)
finalN

array([[0.44870017, 0.4427595 , 0.44371789, 0.40882925, 0.4017102 ,
        0.37083346, 0.42339721]])

모든 변수들을 포함한 회귀모델을 생성하고 그 모델에 의한 MSE를 계산합니다.

mod=LinearRegression().fit(Xtr, ytr)
mseTr=mean_squared_error(ytr, mod.predict(Xtr))
mseTr.round(3)

885.434

위 결과 MSE를 기반으로 AIC, BIC를 계산하기 위한 함수(식 7와 9)를 작성하였습니다.

def AicBic(mse, featureN, sampleN):
    aic=sampleN*np.log(mse)+2*featureN
    bic=sampleN*np.log(mse)+np.log(sampleN)*featureN
    return(aic, bic)

aicT, bicT=AicBic(mseTr, Xtr.shape[0], Xtr.shape[1])
print("aic: %.3f, bic: %.3f" %(aicT, bicT))

aic: 1227.503, bic: 1195.590

위 데이터의 설명변수의 변화에 따른 AIC와 BIC의 변화는 다음과 같습니다.

aic=[]
bic=[]
for i in range(2, 7):
    mod1=LinearRegression().fit(Xtr[:, :i], ytr)
    mse1=mean_squared_error(ytr, mod1.predict(Xtr[:,:i]))
    r, c=Xtr[:,:i].shape
    a, b=AicBic(mse1, r, c)
    aic.append(a)
    bic.append(b)

fig, ax=plt.subplots(figsize=(4, 3))
ax.plot([2,3,4,5,6], aic, color="b", label="AIC")
ax2=ax.twinx()
ax2.plot([2,3,4,5,6], bic, color="r", label="BIC")
ax.set_xlabel("변수의 수")
ax.set_ylabel("AIC")
ax2.set_ylabel("BIC")
ax.legend(loc=(0.1, 0.8), frameon=False)
ax2.legend(loc=(0.1, 0.65), frameon=False)
plt.show()

위 결과는 변수의 수에 따른 영향은 BIC에서 더욱 민감하게 나타납니다.