R의 rjags를 이용한 베이시언Bayesian 입문 1

패키지 설치

JAGS를 먼저 컴퓨터에 설치한다.
현재 버전은 4.3.0
http://mcmc-jags.sourceforge.net/

R에서 JAGS를 사용하기 위하여 rjags를 R에 설치한다.
R console에서 아래 스크립트를 실행

> install.packages("rjags")

데이터 내려받기

평균과 표준편차를 알지 못하는 가상의 데이터 example1.csv를 내려받아 R의 working directory에 저장한다.

example1.csv

JAGS 모델 작성

이제 평균과 표준편차를 모르는 표준정규분포를 따르는 example1.csv 데이터의 평균과 표준편차를 베이시언Bayesian 추론하여보자.

JAGS 모델은 다음과 같다.

1 2 3 4 5 6 7 8

model {     for (i in 1:N) {         x[i] ~ dnorm(mu, tau)     }     mu ~ dnorm(0, 1/100^2)     tau <- 1/simga^2     sigma ~ dunif(0, 100)  }

JAGS는 항상 모델이 무엇인지 말해줘야 하는데 2~4번째 줄에서 평균이 mu이고 표준편차가 1/tau인 정규분포라고 알려주었다. tau는 JAGS에서는 precision이라고 말하는데 단순이 1/(variance)이다.

5번째 줄부터 prior를 넣어주었다. 평균 mu는 평균이 0이고 표준편차가 100인 정규분포이고 분포가 너무 넓어서 non-informative prior이다. 6번째 줄은 tau가 단순히 분산의 역이라 알려주었는데 deterministic function이므로 "~" 대신 "<-"을 사용하여야 한다. 7번째 줄에서 sigma는 0부터 100까지 uniform prior라고 알려주었다.

위 모델을 알기 쉽게 prior를 먼저 쓰고  모델을 나중에 쓸수도 있다.

1 2 3 4 5 6 7 8

model {     mu ~ dnorm(0, 1/100^2)     sigma ~ dunif(0, 100)     tau <- 1/simga^2     for (i in 1:N) {         x[i] ~ dnorm(mu, tau)     }  }

텍스트 편집기를 이용하여 위 모델을 작성하고  example1.jags라고 R의 working directory 저장한다.

rjags을 이용한 R script 작성

R script editor를 사용하여 아래 스크립트를 만든다.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

setwd("d:/tmp/rcode/")  rm(list = ls())  x <- read.csv("example1.csv", header=FALSE)  x <- as.vector(x[,1])  N <- length(x)  library(rjags)  jags <- jags.model(    'example1.jags',    data = list('x' = x, 'N' = N),    n.chains = 4,    n.adapt = 100  )  update(jags, 1000)  jags.results <- jags.samples(   jags,   c('mu', 'tau'),   1000  )

1째줄: 작업디렉토리 바꾸기

2째줄: R의 모든 메모리 삭제

3째줄: 데이터 불러오기

4... : data.frame을 numeric vector로 바꾸기. 이걸 안하면 JAGS에서 인식을 못한다.

5... : 데이터 크기

6... : rjags 라이브러리 불러오기

7-12 : JAGS model object을 작성한다. example.jags는 조금전에 만들었던 JAGS모델 파일이름이고, 어느데이터를 이용할지 data로 넘겨주었다. 여기서는 꼭 list형으로 넘겨주어야만 한다. n.chains는 parallel chains의 갯수이고 n.adapt는 adaptive sampling을 할동안 얼마나 많은 샘플을 버리는지를 말해준다.

13... : 1000번 반복계산 (이와같은 간단한 문제에서는 꼭 필요하지는 않다)

14-18 : 1000개 샘플을 이용하여  mu와 tau를 구한다.

BayesTut1.R으로 위 스크립트를 R의 작업디렉토리에 저장한다.

R console에서 아래와 같이 입력한다.

> source('D:/tmp/rcode/Bayestut1.R')

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1000
   Unobserved stochastic nodes: 2
   Total graph size: 1011

Initializing model

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
  |**************************************************| 100%
  |**************************************************| 100%

위에서 d:/tmp/rcode/는 자신의 R working directory로 바꿔어야만 한다.

> names(jags.results)
[1] "mu"  "tau"

> jags.results$mu
mcarray:
[1] 0.09570689

Marginalizing over: iteration(1000),chain(4)

> jags.results$tau
mcarray:
[1] 0.04038518

Marginalizing over: iteration(1000),chain(4)

평균이 0.096이고 표준편차가 4.980의 데이터라고 추정하였다.

> 1/sqrt(jags.results$tau)
mcarray:
[1] 4.979855

사실 example1.csv는 R에서 rnorm(1000, 0, 5)를 이용하여 만든 데이터인데 1000개만 가지고 그 정확한 값을 추정하기란 무리가 있다.