(파이썬 코드는 맨 밑에 있습니다. 참고하세요)
이전 글에서 AIC와 BIC에 대해 다루었습니다. 이들은 일반적으로 널리 사용되는 최대가능도법 maximum likelihood 으로 모수치를 추정했을 때 사용할 수 있는 방법으로, 최대화된 로그가능도에 -2를 곱한 값에 일종의 페널티를 붙인 뒤, 그 값이 가장 작은 모델을 선택하는 방식으로 작동한다고 했습니다. 이 원칙 자체는 베이지안 추론에서도 똑같이 적용됩니다. 하지만 한 가지 문제가 있습니다. 베이지안 추론에서는 일반적으로 최대가능도법을 사용하는 대신, 사후분포 전체를 추정해서 여기서 평균을 계산한다든지 합니다 (물론 maximum a posterior 같은 방법을 사용하지 않는 건 아닙니다만). 그러면 AIC나 BIC를 계산하는 데 필요한 최대가능도 추정치 (MLE) 가 없습니다. 이럴 땐 어떻게 해야 할까요?
두 가지 해법이 있습니다. 첫 번째 해법은 MLE 대신 모수치의 사후평균 posterior mean 을 이용하는 것이고, 두 번째 해법은 이와 달리 사후분포로부터 얻어진 모수치들에 대해 각각 로그가능도를 구한 뒤, 그 평균을 내는 것입니다. 그러니까 두 방법의 차이는 먼저 평균을 내고 이것을 로그가능도에 대입할 것이냐, 아니면 로그가능도에 일단 모수치를 대입한 뒤 그 평균을 낼 것이냐 하는 차이입니다. 첫 번째 해법은 deviance information criterion (DIC) 이라 불리고, 두 번째 해법은 Watanabe-Akaike information criterion (WAIC)라 불립니다. [1] 전자가 먼제 제안되었고 더 널리 쓰이긴 하지만, 후자가 이론적으로 더 나은 성질들을 갖고 있다고 알려져 있습니다. 그리고 WAIC는 DIC와 달리 fully Bayesian한 방식이라고 소개되는데, 그 이유는 하나의 값, 즉 모수치의 점추정치에 의존하는 값이 아니라 사후분포 전체에 대해 평균이 취해진 값이기 때문입니다.
DIC의 공식은 다음과 같습니다 (사실 보시면 AIC와 매우 유사하다는 것을 아실 수 있을겁니다):
여기서 \(\hat{\theta}\) 는 사후분포로부터 얻은 모수치의 사후평균 posterior mean 을 의미하며, \(p_{DIC}\)는 모수치의 ‘갯수’를 의미합니다. 하지만 여기서 ‘갯수’는 AIC에서와 같은 단순한 갯수를 의미하는 것은 아닙니다. 대신 ‘유효한 모수치의 갯수’ effective number of parameters 인데, 이는 모수치들의 추정치 사이에 상관이 있음을 감안하기 위한 것입니다. 이를테면 어떤 모델의 모수치들의 추정치 사이에 0.99 가량의 상관이 존재한다면, 이들을 독립적인 모수치로 간주하는 것은 다소 말이 안 되겠죠. 사실 어쩌면 이들은 한 몸처럼 함께 움직이는 모수치들일 수도 있고, 따라서 이들을 하나로 간주하는 게 더 옳을지도 모릅니다. (특히 이런 상황은 모수치들 간 의존도가 심한 hierarchical model에서 흔히 발생합니다.) \(p_{DIC}\)는 이런 사실을 반영하여, 모델에 ‘실질적인’ 파라미터가 몇 개가 있는지 계산한 수치입니다. 자세한 계산법은 참고문헌 또는 아래 코드를 보세요.
WAIC도 비슷한 아이디어에 의존하지만, 공식이 복잡하므로 여기서 소개하지는 않겠습니다. (굳이 계산과정이 궁금하신 분들은 참고문헌과 맨 밑에 있는 별도의 코드를 보세요.) 다만 \(\hat{\theta}\)를 계산하여 대입하는 대신 사후분포를 이용하여 먼저 density를 계산한 다음 평균을 낸 것을 가지고 로그 가능도를 계산한다는 것과, effective number of parameters의 계산법이 조금 다르다는 것만 지적하겠습니다.
여기서는 WAIC를 실제로 계산해 보겠습니다. 이젠 익숙하시겠지만, 또 iris 자료의 3번째 열 (Petal.Length) 을 첫/두번째 열 (Sepal.Length, Sepal.Width) 로 예측하는 선형회귀모형을 만들어 보겠습니다. 핏팅에는 rstan 패키지를 사용하며, 코드는 아래와 같습니다. 이번에는 맨 마지막에 generated quantities 블록이 새로 생긴 걸 볼 수 있는데, 이렇게 log-likelihood를 데이터 포인트마다 저장해주어야 WAIC의 계산이 가능합니다.
library(rstan)
data <- list(N = nrow(iris),
X = cbind(1,iris[,1:2]),
y = iris[,3])
model_code <- '
data {
int N;
matrix[N,3] X;
vector[N] y;
}
parameters {
real<lower=0> sigma2;
vector[3] beta;
}
transformed parameters {
real<lower=0> sigma;
sigma = sqrt(sigma2);
}
model {
sigma2 ~ inv_gamma(.001, .001);
beta ~ normal(0, 1000);
y ~ normal(X*beta, sigma);
}
generated quantities {
vector[N] log_lik;
for(i in 1:N) log_lik[i] = normal_lpdf(y[i] | X[i,]*beta, sigma);
}
'
fit <- stan(model_code=model_code, data=data)
별도의 패키지인 loo는 rstan의 fit 오브젝트로부터 waic를 계산할 수 있는 waic() 함수를 제공합니다. 이것으로 waic를 계산하면 다음과 같습니다:
library(loo)
log_lik <- extract_log_lik(fit)
waic <- waic(log_lik)
waic
Computed from 4000 by 150 log-likelihood matrix
Estimate SE
elpd_waic -149.8 8.2
p_waic 3.6 0.5
waic 299.6 16.4
\(p_{WAIC}\)는 WAIC가 사용하는 effective number of parameters인데, 보면 그 추정치가 3.6입니다. 이것이 실제 모수치 갯수인 4개 (회귀계수 세 개, 분산 한 개) 보다 작음에 주목하시기 바랍니다. 이것은 모수치들 간에 상관이 존재하며, 따라서 완전히 별도의 모수치로 볼 수 없다는 것을 반영합니다. 계산된 WAIC는 299.6입니다.
모형비교의 예를 들기 위해, 위 모델에서 2번째 열을 predictor에서 제외하고 다시 모형을 만들어 WAIC를 계산하면 다음과 같습니다:
data <- list(N = nrow(iris),
X = cbind(1,iris[,1]),
y = iris[,3])
model_code2 <- '
data {
int N;
matrix[N,2] X;
vector[N] y;
}
parameters {
real<lower=0> sigma2;
vector[2] beta;
}
transformed parameters {
real<lower=0> sigma;
sigma = sqrt(sigma2);
}
model {
sigma2 ~ inv_gamma(.001, .001);
beta ~ normal(0, 1000);
y ~ normal(X*beta, sigma);
}
generated quantities {
vector[N] log_lik;
for(i in 1:N) log_lik[i] = normal_lpdf(y[i] | X[i,]*beta, sigma);
}
'
fit2 <- stan(model_code=model_code2, data=data)
log_lik2 <- extract_log_lik(fit2)
waic2 <- waic(log_lik2)
waic2
Computed from 4000 by 150 log-likelihood matrix
Estimate SE
elpd_waic -193.3 8.3
p_waic 2.4 0.4
waic 386.6 16.6
이번에도 \(p_{WAIC}\)는 2.4로, 실제 모수치 갯수인 3개보다 적습니다. 앞 모형의 WAIC를 비교해 보면 앞 모형의 WAIC는 299.6, 이번에 만든 모형의 WAIC는 386.6입니다. WAIC가 더 작은 쪽, 즉 첫 번째 모형을 선택하면 되겠습니다.
지금까지 DIC와 WAIC에 대해 살펴보았습니다. 생각보다 간단하게(?) WAIC를 계산할 수 있었죠? 사실 이 외에도 stan과 loo 패키지는 leave-one-out cross-validation (LOO-CV) 과 같은 다른 모형비교 기준들도 제공하기 때문에, 편리하게 사용하실 수 있습니다. 계산량이 작을 때는 LOO-CV, 클 때는 그 근삿값으로 WAIC를 사용할 것이 권장됩니다.
[부록] 여기서부터는 DIC와 WAIC의 구체적인 계산 과정이 궁금하신 분들을 위한 부분입니다. 위 참고문헌 [1]에 나온 공식을 참조하여 DIC와 WAIC를 계산하여 원래의 값과 비교해 보았습니다. 관심 있으신 분들은 참조하시기 바랍니다.
# R
library(loo)
library(rstan)
# 데이터 준비
data <- list(N = nrow(iris),
X = cbind(1,iris[,1:2]),
y = iris[,3])
# 모델 핏팅
model_code <- '
data {
int N;
matrix[N,3] X;
vector[N] y;
}
parameters {
real<lower=0> sigma2;
vector[3] beta;
}
transformed parameters {
real<lower=0> sigma;
sigma = sqrt(sigma2);
}
model {
sigma2 ~ inv_gamma(.001, .001);
beta ~ normal(0, 1000);
y ~ normal(X*beta, sigma);
}
generated quantities {
vector[N] log_lik;
for(i in 1:N) log_lik[i] = normal_lpdf(y[i] | X[i,]*beta, sigma);
}
'
fit <- stan(model_code=model_code, data=data)
# likelihood와 log-likelihood 추출
post_parms <- extract(fit)
post_loglik <- extract(fit)$log_lik
post_lik <- exp(post_loglik)
# 모수치들의 사후평균 계산 (DIC 계산용)
postmean_beta <- colMeans(post_parms$beta)
postmean_sigma <- mean(post_parms$sigma)
loglik_postmean <- sum(dnorm(data$y, as.matrix(data$X) %*% postmean_beta, postmean_sigma, log=T))
# p_DIC1이나 p_DIC2 중 하나를 사용할 것
# 논문 섹션 3.3 참고
#p_DIC1 <- 2*(loglik_postmean - mean(apply(post_loglik, 1, sum)))
p_DIC2 <- 2*var(apply(post_loglik, 1, sum))
elpd_DIC <- loglik_postmean - p_DIC2
DIC <- -2*elpd_DIC
# p_WAIC1이나 p_WAIC2 중 하나를 사용할 것
# 논문 섹션 3.4 참고
#p_WAIC1 <- 2*(sum(log(colMeans(post_lik)) - colMeans(post_loglik)))
p_WAIC2 <- sum(apply(post_loglik, 2, var))
lppd <- sum(log(apply(post_lik, 2, mean)))
elppd_WAIC <- lppd - p_WAIC2
WAIC <- -2*elppd_WAIC
# loo 패키지의 계산 결과와 비교
log_lik <- extract_log_lik(fit)
waic <- waic(log_lik)
waic
DIC
[1] 299.8967
WAIC
[1] 299.4765
waic
Computed from 4000 by 150 log-likelihood matrix
Estimate SE
elpd_waic -149.7 8.2
p_waic 3.6 0.5
waic 299.5 16.4
# Python
import numpy as np
from scipy.stats import norm
import pystan
from sklearn.datasets import load_iris
# 데이터 불러오기
iris = load_iris()
X = np.hstack((np.ones(150).reshape(150,1),iris['data'][:,[0,1]]))
y = iris['data'][:,2]
data = {'N' : len(X), 'X': X, 'y': y}
# Stan 코드
model_code = """
data {
int N;
matrix[N,3] X;
vector[N] y;
}
parameters {
real<lower=0> sigma2;
vector[3] beta;
}
transformed parameters {
real<lower=0> sigma;
sigma = sqrt(sigma2);
}
model {
sigma2 ~ inv_gamma(.001, .001);
beta ~ normal(0, 1000);
y ~ normal(X*beta, sigma);
}
generated quantities {
vector[N] log_lik;
for(i in 1:N) log_lik[i] = normal_lpdf(y[i] | X[i,]*beta, sigma);
}
"""
# 모델 핏팅하기
model = pystan.StanModel(model_code=model_code)
fit = model.sampling(data=data)
# 사후표본 추출하기
post = fit.extract()
# likelihood와 log-likelihood 추출
post_loglik = post['log_lik']
post_lik = np.exp(post_loglik)
# 모수치들의 사후평균 계산 (DIC 계산용)
postmean_beta = np.mean(post['beta'], 0)
postmean_sigma = np.mean(post['sigma'])
loglik_postmean = np.sum(np.log(norm.pdf(data['y'], np.matmul(data['X'], postmean_beta), postmean_sigma)))
# p_DIC1이나 p_DIC2 중 하나를 사용할 것
# 참고문헌 [1] 섹션 3.3 참고
p_DIC1 = 2*(loglik_postmean - np.mean(np.sum(post_loglik, 1)))
p_DIC2 = 2*np.var(np.sum(post_loglik, 1))
elpd_DIC = loglik_postmean - p_DIC2
DIC = -2*elpd_DIC
print(DIC)
# p_WAIC1이나 p_WAIC2 중 하나를 사용할 것
# 위 논문 섹션 3.4 참고
p_WAIC1 = 2*(np.sum(np.log(np.mean(post_lik,0)) - np.mean(post_loglik,0)))
p_WAIC2 = np.sum(np.apply_along_axis(np.var, 0, post_loglik))
lppd = np.sum(np.log(np.mean(post_lik, 0)))
elppd_WAIC = lppd - p_WAIC2
WAIC = -2*elppd_WAIC
print(WAIC)
299.8930147646932
299.6463113027089