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

상세한 모델 설계를 할 수 있는 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)

#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')

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

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)

MCMCpack 이용 시뮬레이션

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

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

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")

몇가지 베이지언 계산 방법 정리 by from __future__ import dream is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.