• Home

베이지언 회귀모델 시뮬레이션

1885년 Galton이 조사한 부모의 키와 자식의 키의 데이터를 기반으로 몇가지 회귀모델을 작성해 보았으며, 개인적으론 전통적인 회귀모델과 베이지언 관점의 모델의 비교를 해보고자 했다.

작업 방식은 아래와 같다.

  1. 전통적인 least squares기반의 회귀모형 생성 그리고 플로팅
  2. MCMCpack을 이용한 Bayesian Linear Regression
  3. rjags를 이용한 모델링 및 시뮬레이션
    • Likelihood를 Normal 분포로
    • Likelihood를 T분포로..

그리고 MCMC결과물로 prediction을 하고 prediction 결과분포를 도식화 해보는 실습을 해보았다.

library(UsingR)
library(ggplot2)
library(MCMCpack)
library(rjags)
library(hdrcde)
library(ggmcmc)
data(galton)

#EDA 
ggplot(galton, aes(child, parent)) + geom_point(position=position_jitter()) + stat_smooth(method='lm')

plot of chunk init

#Linear Regression 
summ <- lm(parent ~ child, data=galton)
summary(summ)
## 
## Call:
## lm(formula = parent ~ child, data = galton)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.670 -1.170 -0.147  1.132  4.272 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  46.1353     1.4122    32.7   <2e-16 ***
## child         0.3256     0.0207    15.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.59 on 926 degrees of freedom
## Multiple R-squared:  0.21,   Adjusted R-squared:  0.21 
## F-statistic:  247 on 1 and 926 DF,  p-value: <2e-16

전통적인 방식의 Linear Regression 결과를 도출하고…

#using MCMCpack 
summary(MCMCregress(parent ~ child, data=galton))
## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean     SD Naive SE Time-series SE
## (Intercept) 46.137 1.4117 0.014117       0.014117
## child        0.326 0.0207 0.000207       0.000207
## sigma2       2.531 0.1181 0.001181       0.001181
## 
## 2. Quantiles for each variable:
## 
##               2.5%    25%    50%   75%  97.5%
## (Intercept) 43.365 45.184 46.147 47.09 48.908
## child        0.285  0.312  0.325  0.34  0.367
## sigma2       2.309  2.450  2.527  2.61  2.769

MCMCregress기반 시뮬레이션 요약을 확인한 다음 ….

#Normal Likelihood 기반 회귀모형 
model <- "
model {
    for(i in 1:length(x) ) {
        y[i] ~ dnorm( mu[i] , tau )
        mu[i] <- beta0 + beta1 * x[i]
    }
    beta0 ~ dnorm( 0 , 1.0E-12 )
    beta1 ~ dnorm( 0 , 1.0E-12 )
    tau ~ dgamma( 0.001 , 0.001 )
}
"

#초기 모델 생성 체인은 3개 
jagsModel = jags.model(textConnection(model) , data=list(x = galton$child ,y = galton$parent),
                        n.chains=3 , n.adapt=10000, quiet=TRUE)
#burnin 과정 
update(jagsModel, n.iter=10000 ,progress.bar = 'none')

#각 체인당 10만개의 샘플을 가져온다. 
codaSamples <- coda.samples(jagsModel, variable.names=c("beta0" , "beta1" , "tau")   , 
                            n.iter=500000 , thin=5, progress.bar = 'none')



ggs_density(ggs(codaSamples))

plot of chunk jags1

#각 체인에서 200쌍의  파라멘터 기반으로 도식화, 체인의 색은 다르게....
xseq <- seq(range(galton$child)[1],range(galton$child)[2], length=100)

sampseq <- seq(1,nrow(codaSamples[[1]]), by=2000)

basep <- ggplot() + geom_point(data=galton, aes(child,parent ),position=position_jitter())

colours <- c("hotpink", "skyblue", "lightyellow")
for(j in 1:3){
  for(i in sampseq){
    predy <- codaSamples[[j]][i, 'beta0'] + codaSamples[[j]][i, 'beta1'] *  xseq
    basep <- basep %+% geom_line(data=data.frame(x=xseq, y=predy), aes(x,y), colour=colours[j], alpha=0.5)
  }
}

#각 파라메터의 분포가 대칭이나 평균값을 대표값으로 사용 대표 모델을 도식화 함 
predy <- mean(codaSamples[[1]][, 'beta0']) + mean(codaSamples[[1]][, 'beta1']) * xseq
basep <- basep %+% geom_line(data=data.frame(x=xseq, y=predy), aes(x,y), colour="black", size=1)

#Beta0
print(mean(codaSamples[[1]][, 'beta0']))
## [1] 46.08
#Beta1
print(mean(codaSamples[[1]][, 'beta1']))
## [1] 0.3265
#10개의 선택된 x값의 y에 대한 95% HDR 분포 도식  
xseqpred <- seq(range(galton$child)[1],range(galton$child)[2], length=10)

dfres <- data.frame(x=c(), ymin=c(), ymax=c())

for(i in xseqpred){
  predy <- codaSamples[[1]][, 'beta0'] + codaSamples[[1]][, 'beta1'] * i
  hdres <- hdr(predy)
  dfres <- rbind(dfres, data.frame(x=i,ymin=hdres$hdr[2,][1], ymax=hdres$hdr[2,][2]))
}

basep %+% geom_errorbar(data=dfres,aes(x=x, ymin=ymin, ymax=ymax),width=0.2, linetype=5, colour='red') + ggtitle("Simple Linear Regression with MCMC")

plot of chunk jags1

아래는 Likelihood를 T분포를 사용한 로버스트 회귀 모형이다.

model2 <- "
model {
    for( i in 1 : length(x) ) {
        y[i] ~ dt( mu[i] , tau , tdf)
        mu[i] <- beta0 + beta1 * x[i]
    }
    beta0 ~ dnorm( 0 , 1.0E-12 )
    beta1 ~ dnorm( 0 , 1.0E-12 )
    tau ~ dgamma( 0.001 , 0.001 )
    udf ~ dunif(0, 1)
    tdf <- 1 - tdfGain * log(1-udf)
}
"


jagsModel2 = jags.model(textConnection(model2) , 
                        data= list(x = galton$child, y = galton$parent, tdfGain = 100) ,
                        n.chains=3 , n.adapt=10000, quiet=TRUE)


update( jagsModel2 , n.iter=10000, progress.bar='none')

codaSamples2 = coda.samples( jagsModel2 , variable.names= c("beta0" , "beta1" , "tau", "tdf") , 
                            n.iter=50000 , thin=4, progress.bar = 'none')

ggs_density(ggs(codaSamples2))

plot of chunk jags2

xseq <- seq(range(galton$child)[1],range(galton$child)[2], length=100)

sampseq <- seq(1,nrow(codaSamples2[[1]]), by=50)


basep <- ggplot() + geom_point(data=galton, aes(child, parent),position=position_jitter())

colours <- c("hotpink", "skyblue", "lightyellow")
for(j in 1:3){
  for(i in sampseq){
    predy <- codaSamples2[[j]][i, 'beta0'] + codaSamples2[[j]][i, 'beta1'] *  xseq
    basep <- basep %+% geom_line(data=data.frame(x=xseq, y=predy), aes(x,y), colour=colours[j], alpha=0.5)
  }
}

predy <- mean(codaSamples2[[1]][, 'beta0']) + mean(codaSamples2[[1]][, 'beta1']) * xseq
basep <- basep %+% geom_line(data=data.frame(x=xseq, y=predy), aes(x,y), colour="black", size=1)

print(mean(codaSamples2[[1]][, 'beta0']))
## [1] 46.57
print(mean(codaSamples2[[1]][, 'beta1']))
## [1] 0.3193
print(basep)

plot of chunk jags2

xseqpred <- seq(range(galton$child)[1],range(galton$child)[2], length=10)

dfres <- data.frame(x=c(), ymin=c(), ymax=c())

for(i in xseqpred){
  predy <- codaSamples2[[1]][, 'beta0'] + codaSamples2[[1]][, 'beta1'] * i
  hdres <- hdr(predy)
  dfres <- rbind(dfres, data.frame(x=i,ymin=hdres$hdr[2,][1], ymax=hdres$hdr[2,][2]))
}

basep %+% geom_errorbar(data=dfres,aes(x=x, ymin=ymin, ymax=ymax),width=0.2, linetype=5, colour='red') + ggtitle("Robust Regression with MCMC")

plot of chunk jags2

여러 방식으로 단순 회귀 모형을 만들어본 결과 베이지언 방식의 특징인 사후 분포(y)를 직접 확인 가능하다는 특징을 잘 보여주고 있음을 알 수 있다. 엄밀히 따지면 베이지언 회귀모형은 베이지언 계층 모델의 한 예일 뿐이다. 모델의 모든 부분이 prior, posterior 그리고 hyperprior로 해석 될 수 있으니 말이다.

사전 분포로 정해준 값들에서 모두 의미 있게 떨어져 있는 값이기 때문에 분포 자체가 의미가 있으며 전통적인 회귀분석 방법론의 추정 변수의 유의미성을 검정하는 F검정과 같은 방식이 단일값 best fit값의 유의미성 검정이라는 것을 볼때 우리가 사후 분포로 추정한 모델의 분포를 알게 됨으로써 이 모델을 동적으로 업데이틀 할 수 있는 여지가 생기게 된다. 따라서 베이지언으로 어떤 모델을 만들던 새로운 데이터를 기반으로 기존 모델을 업데이트 할 수 있다는 것인데, 이는 모델 자동화의 큰 포인트 중에 하나이며, 시계열 적인 특징을 가장 자연스럽게 적용할 수 있는 통로가 된다. 심지어 위에서 만든 베이지언 회귀모형도 새로운 데이터를 기반으로 과거 모수를 업데이트 할 수 있다.

전통적인 회귀분석 역시 신뢰구간이라는 것을 가질 수 있으나 모수가 해당 영역에 존재할 확률이 95%라는 것과 모수의 분포를 가질 수 있는것과는 큰 의미적 차이와 관점의 차이가 있다. 이런 이유때문에 베이지언은 알고리즘이 아니라 페러다임이라는 이야기가 맞는거 같다.

개인적인 생각으로는 빈도주의자들의 신뢰구간은 모수가 특정 (단일)값이라는 가정이 기본적으로 깔려 있고, 베이지언의 경우 모수는 특정 분산을 가진 분포라는 가정이 가장 기본적인 차이가 아닐까하는 생각도 해본다. 아무리 데이터를 넣어도 단일 값으로 떨어질 수 없는 것이 베이지언 모수라는 이야기다.

베이지언 책에 대한 이야기

뭔가 새로운걸 배울땐 여러 책을 구입해 겹쳐읽기를 하는 패턴으로 학습을 하곤한다. 이 방법의 장점은 같은 주제의 설명을 다른 저자로부터 들을 수 있다는 것이고 이 덕분에 빨리 제대로된 이해를 할 수 있게 된다. 

이번 전반기에 구입한 베이지언 책들은 크게 국내서와 해외서로 나눌 수 있는데 개인적으로 가장 만족감을 느꼈던 국내서 한권과 해외서 한권을 소개한다. 물론 개인의 여건에 따라 다른 책들이 더 의미있게 다가올 수 있을거라 생각하며 전적으로 개인적인 평임을 강조하는 바이다. 

국내서 한권은 오만숙 교수님의 베이지안 통계추론이다. 다른 베이지언 책과 유사하게 수식이 상당히 많으나 모두 필요한 설명이며, 적절한 코드 예시는 실제 구현을 하는데 큰 도움이 되었다.  필자는 이 책의 연습문제까지 풀어볼 정도로 한때 이 책을 끼고 살았으며 상당히 만족하면서 책을 탐독했다.  뒷 부분으로 갈수록 앞에서 느꼈던 간략하고 명확한 설명이 부족한 측면이 있었으나 그런 부분은 개인적으로 더 공부해야될 부분들이라 생각에 크게 신경쓰지 않았다.  베이지언 계층(Hierarchical)모델 이전까지의 기초를 확실하게 잡는데 가장 탁월한 책이라 생각한다. 기회가 된다면 교수님 강의도 들어보고 싶다는 강한 욕구도 느낄 수 있었다. 

외국의 경우도 베이지언 책에 대한 추천 논의가 있는데, 이 리스트에서 가장 많은 표를 얻은 Doing Bayesian Data Analysis를 본인 역시 추천한다. 

이 책의 경우 통계학 개론 정도만 뗀 사람이라도 차근차근 탐독한다면 충분히 섭렵할 만한 내용과 난이도를 가지고 있다. 그렇다고 책이 커버하는 영역이 한정적인 것은 절대 아니다. 대부분 베이지언에서 언급되는 내용은 모두 포함하고 있고, 특히나 jags나 Bugs 코드가 제공되고 있어 실제 코드레벨에서의 구현을 바로 해볼 수 있는 장점이 있다. 개인적으로는 수식이 주는 간결함과 명확함이 좋기도 하지만 실제 코드를 실행하고 데이터를 넣어보면서 체내화 되는 그 알고리즘의 느낌 역시 중요하다고 생각한다.  많은 내용을 쉽게 설명해 주려는 저자의 노력이 많은 부분에 숨겨져 있어 보는 내내 감동 비슷한 느낌을 받기도 한다. 

외서의 경우 많은 사람들이 추천하는 Gelman의 책은 정말 이쪽 영역의 바이블같은 느낌을 주기에 충분하나 연구자가 아닌 실무영역에서 베이지언을 사용하는데 너무 많은 내용을 포함하고 있는게 사실이다. 이 책이 주는 그 가치(특히 예제)가 정말 어마하다는 것을 알고 있으나 개념에 대한 직관적인 체내화를 하는데 너무 긴 시간이 걸리는 단점이 있다. 따라서 이 책은 베이지언에 대해서 특정 내용을 찾아볼 수 있는 레퍼런스 북 같은 역할을 하고 있다. 거의 베이지언에 대한 모든 것을 다룬 책이라 베이지언을 하다보면 반드시 옆에 두고 봐야될 책이라 생각하나 처음부터 이 책을 기반으로 공부하는건 아마도 큰 각오를 해야될 것이다. 

책이라 하는건 시간이 지나고 본인의 경험이 쌓여 가면서 선호도가 갈리는 특징을 가지고 있다. 따라서 지금 현재 자신에게 가장 맞는 책을 만나는게 뭔가를 처음 배울때 가장 중요한 부분이라 생각한다. 따라서 지금 올리는 글 조차도 내일이 되면 다른 느낌으로 다가올지 모르겠지만 뭔가를 배우는 사람 입장에서는 오히려 반겨야 될 상황이 아닐까 한다. 

몇가지 베이지언 계산 방법 정리

지금까지 알고 있는 몇가지 방식의 베이지언 계산법 정리를 해볼 필요가 있어서 같은 문제를 여러 방법으로 살펴봤다.

상세한 모델 설계를 할 수 있는 jags 또는 좀더 빠르다고 하는 stan을 좀더 익숙하게 쓸 수 있도록 좀 살펴볼 필요가 있을거 같다. 그리고 여러 실 활용 예들을 찾아서 책을 좀 뒤적이도록 하자!

Grid Approximation

library(rjags)
library(MCMCpack)
library(extrafont)
library(ggthemes)
library(ggmcmc)
library(knitr)
#Grid Approximation 

m <- 0.5
sd <- 0.1
samps <- rnorm(10000, m, sd)

tbl <- table(round(samps, 3))

tbl <- tbl[as.numeric(names(tbl)) >  0 & as.numeric(names(tbl)) <   1]

theta <- as.numeric(names(tbl))

#theta의 사전 분포 
pTheta <- as.numeric(tbl/sum(tbl))

plot(theta, pTheta, main='prior', pch=20)

plot of chunk unnamed-chunk-3

#14번 던져 번 11앞면 
N <- 14
h_t <- c(1,1,1,1,1,1,1,1,1,1,1,0,0,0)
z <- length(which(h_t == 1))
#likelihood
pDataGivenTheta <- theta^z * (1-theta)^(N-z)

plot(theta, pDataGivenTheta, pch=20, main='likelihood')

plot of chunk unnamed-chunk-3

posterior_grid <- (pTheta * pDataGivenTheta)/sum(pTheta * pDataGivenTheta)
par(family='Bitstream Vera Sans Mono')
plot(theta, posterior_grid, main="Posterior", pch=20)

plot of chunk unnamed-chunk-3

Beta Prior and Posterior 기반

#using beta conjugate prior 
prior_beta_alpha <- m*(m*(1-m)/(sd^2) - 1)
prior_beta_beta <- (1-m)* (m*(1-m)/(sd^2) - 1)

posterior_conjugate <- dbeta(theta, z + prior_beta_alpha,N - z + prior_beta_beta)

plot(theta, posterior_conjugate, main="using beta prior", pch=20)

plot of chunk unnamed-chunk-4

MCMCpack 이용 시뮬레이션

posterior_mc <- MCbinomialbeta(z, N,alpha =prior_beta_alpha, beta =prior_beta_beta,   mc = 100000)

plot(density(posterior_mc), main="using MCbinomialbeta function")

plot of chunk unnamed-chunk-5

jags 이용 시뮬레이션

#jags 방식 
model_str <- "
model{
  for(i in 1:N){
    y[i] ~ dbern(theta)
  }
  theta ~ dnorm(m,1/pow(sd,2))T(0.0001, 0.99)
}
"

mdl <- jags.model(file=textConnection(model_str), 
                  data=list(N=N, y=h_t, m=m, sd=sd),quiet = T)

# Burn-in:
update(mdl, n.iter=10000,progress.bar = 'none')

sampled <- coda.samples(mdl, thin = 10,n.iter = 300000,variable.names = 'theta',progress.bar = 'none')


ggs_density(ggs(sampled)) %+% ggtitle("jags 기반 시뮬레이션") %+% theme_solarized_2(base_family = "Bitstream Vera Sans Mono")

plot of chunk unnamed-chunk-6