[R] Nimble 패키지를 이용한 Bayesian 모수 추정

[1] Nimle 이란?

MCMC방식을 이용해 구하고자 하는 분포의 모수를 추정하는 방식인 베이지안 모수 추정 방식을 구해주는 R 패키지 중 하나임. (다른 것으론 STAN, JAGS도 있음)

 

[2] Nimble 사용방법

※ 본 내용에선 R(내지는 Rstuio)가 설치되어있음을 가정하고 진행합니다. 만일 R 자체가 설치되어있지 않다면 이 링크를 따라 진행해주세요!

 

전체 진행 순서

nimble 패키지 설치 > Rtools 설치 > 구현

 

1) R console창에 nimble 패키지 설치

install.packages("nimble")

 

2) R version확인

R.version

R의 버전 확인

Rtools를 설치하기전 현재 4.2.1의 R version을 확인해 이에 맞는 버전을 설치해야합니다.

 

3) Rtools 설치 방법

• Rtools설치 홈페이지(https://cran.r-project.org/bin/windows/Rtools/)에 접속합니다.

Rtools 4.2링크 클릭

• 현재 4.2.x 버전이므로 RTools 4.2를 클릭해 접속한 뒤, 아래 링크를 클릭해 설치파일 다운을 진행합니다.

• 설치파일 .exe 더블클릭해 실행을 눌러 설치를 진행합니다.

• Next 버튼을 쭉 누른뒤 Install을 누르면 완료됩니다.

 

 

실습 코드

 

library(nimble)
library(MCMCvis)


X = model.matrix(n~TypeCity+TypeCounty+TypeSchool+TypeTown+TypeVillage+C1+C2, data = data.train)


model <- nimbleCode({
  # likelihood
  for(j in 1:JJ){
    lambda[j] <-exp(inprod(X[j,1:KK], beta.hat[1:KK]))
    y[j] ~ dnegbin(prob = 2/(2+lambda[j]*R[Person[j]]), size = 2)
  }
  
  # prior of R
  for(n in 1:N.person){
    R[n] ~ dgamma(mean=1, sd=1/a.hat)
  }
  
  # prior of Betas (params 개수 만큼 추출)
  for(j in 1:KK){
    beta.hat[j] ~ dnorm(mean=0, sd=10)
  } 
  
  # prior of alpha
  a.hat ~ dunif(min=0, max=10)
  
})

my.data <- list(X = X, y = data.train$n)
my.constants <- list(JJ=dim(X)[1], KK=dim(X)[2], Person=person, N.person=max(person))
parameters.to.save <- c("beta.hat", "a.hat", "lambda")

pois_model <- nimbleModel(code = model, name = "pois_model", constants = my.constants, data = my.data)
Cpois_model  <- compileNimble(pois_model)
pois_modelConf <- configureMCMC(pois_model, print = TRUE)
pois_modelConf$addMonitors('beta.hat', "a.hat", "lambda")
MCMC <- buildMCMC(pois_modelConf)
CMCMC <- compileNimble(MCMC, project = pois_model)


n.iter <- 10000
n.burnin <- 4000
n.chains <- 2

set.seed(123)
samples <- runMCMC(CMCMC, niter = n.iter,
                   nburnin = n.burnin,
                   nchains = n.chains)

MCMCtrace(object = samples,
          pdf = FALSE, # no export to PDF
          ind = TRUE, # separate density lines per chain
          params = "beta.hat")

MCMCtrace(object = samples,
          pdf = FALSE, # no export to PDF
          ind = TRUE, # separate density lines per chain
          params = "a.hat")


str(samples[[1]])
result <- MCMCsummary(object = samples, round = 2)
dim(result)

# save estimated params
a.hat <- result[1,1] ;a.hat
beta.hat <-  result[2:9,1] ;beta.hat
lambda.hat <- result[10:1476,1] 
summary(lambda.hat)

unique(cbind(person, lambda.hat))