rjags를 이용한 베이시언Bayesian 입문 4

Introduction to Bayesian inference with JAGS   강의 정리 1

channel9.msdn.com/Events/useR-international-R-User-conferences/useR-International-R-User-2017-Conference/Introduction-to-Bayesian-inference-with-JAGS

Introduction to Bayesian inference with JAGS

 

Mean by Bayesian inference

평균이 3.7이고 표준편차가 1인 데이터(y)가 있다고 가정하고
이 데이터의 평균을 Bayesian inference를 이용하여 구한다.

 

가상의 데이터를 생성한다.

rm(list=ls())

## Data
set.seed(1355)
y = rnorm(50, mean = 3.7, sd = 1)

위 생성된 데이터의 평균을 Bayesian으로 추정

## JAGS model
modelstring = "
model{
# Likelihood
  for(i in 1:N) {
    y[i] ~ dnorm(theta, 1)
  }
  
# Prior
  theta ~ dnorm(0, 1)
}
"

## Estimate the mean of the data, y
library(rjags)
m = jags.model(textConnection(modelstring), data=list(y=y, N=length(y)))

# class(m)
# print(m)
# m$data() # identical with data, y

## extract samples
s = coda.samples(m, "theta", n.iter=1000)
s = unlist(s) 

# mean
mean(s)
round(mean(s), digits=2)

# standard error of the mean
sem = sd(s) / sqrt(length(s) - 1)

# posterior desnity plot
plot(density(s))
abline(v=mean(y), col="grey") # sample mean
abline(v = 3.7, col="black") # true mean

# check normality
qqnorm(s)
shapiro.test(s)

 

 

추가: 여러번 실행해도 같은 샘플이 나오도록 하기

inits1=list(".RNG.name"="base::Wichmann-Hill",
            ".RNG.seed"=7453210)
inits2=list(".RNG.name"="base::Wichmann-Hill",
            ".RNG.seed"=33543)

ma <- jags.model("location1.bug", data=list(y=y, N=length(y)), 
                 n.chains=2, inits=list(inits1, inits2))
mb <- jags.model("location1.bug", data=list(y=y, N=length(y)), 
                 n.chains=2, inits=list(inits1, inits2))

## extract samples
sa = coda.samples(ma, "theta", n.iter=1000)
sb = coda.samples(mb, "theta", n.iter=1000)
sa = unlist(sa)
sb = unlist(sb)

# mean
mean(sa) 
mean(sb)

평균 추정값은 아래와 같다.

> mean(sa) 
[1] 3.691397

> mean(sb)
[1] 3.691397

 

 

Shrinkage

prior를 Normal(0, 1)로 정해주었는데 이제 prior의 표준편차를 점점 줄여가면서 posterior 분포의 변화를 알아보자.

JAGS 언어에서 norm함수는 $norm( \mu, \tau) $와 같이 인자가 주어지는데 $\tau$는 $1 / \sigma^2$인 것을 주의해야 한다.

rm(list=ls())

## Data
set.seed(1355)
N=50
y = rnorm(N, mean = 3.7, sd = 1)

## JAGS model
modelstring = "
model{
# Likelihood
  for(i in 1:N) {
    y[i] ~ dnorm(theta, 1)
  }
  
# Prior
#   tau is supplied from jags.model() function
  theta ~ dnorm(0, tau)
}
"

## Estimate the mean of the data, y
library(rjags)
m = jags.model(textConnection(modelstring), data=list(y=y, N=N, tau=1))
s = coda.samples(m, "theta", n.iter=1e6)
s = unlist(s)

# posterior desnity plot
plot(density(s), col="blue", xlim = c(0, 4.5))

m2 = jags.model(textConnection(modelstring), data=list(y=y, N=N, tau=10))
s2 = coda.samples(m2, "theta", n.iter=1e6)
s2 = unlist(s2) 
lines(density(s2), col = "black")

m3 = jags.model(textConnection(modelstring), data=list(y=y, N=N, tau=20))
s3 = coda.samples(m3, "theta", n.iter=1e6)
s3 = unlist(s3) 
lines(density(s3), col = "red")

m4 = jags.model(textConnection(modelstring), data=list(y=y, N=N, tau=20000))
s4 = coda.samples(m4, "theta", n.iter=1e6)
s4 = unlist(s4) 
lines(density(s4), col = "orange")

prior의 표준편차를 점점 줄여감에따라 posterior의 평균이 prior의 평균에 가까워지고 true mean인 3.7에서 점점 벗어난다.