First JAGS model in R

[GaryBAYLOR]

The JAGS can perform MCMC and R package rjags can be used to access JAGS in R. The process is pretty easy. First we write the model file in R (the file can be in .R format), save the file. Then invoke the model using function in the package 'rjags'. The following is an easy example. Step 1. Save the following code with name "myFirstJags.R".

model {
  for(i in 1:J) {
    y[i] ~ dbin(theta[i], 100)
    theta[i] ~ dbeta(alpha, beta)
  }
  alpha ~ dgamma(1, 1)
  beta ~ dgamma(1, 1)
}

Step 2. Run the following code.

set.seed(1001)
theta <- rbeta(10, 1, 5)
set.seed(1002)
y <- rbinom(10, 100, theta)

jagresult <- jags.model("myFirstJags.R", data = list("y" = y, "J" = 10), n.adapt = 1000)
thesamps <- coda.samples(jagresult, c("alpha", "beta"), n.iter = 10000, thin = 1, progress.bar = "gui")

We can use plot function to get the graph of alpha and beta.

plot(thesamps)

We will get the following plot.

But sadly, as I know, JAGS cannot use improper prior!

[GaryBAYLOR]

Now we can compare the results by running JAGS and that by running Metropolis Hasting algorithm from scratch. The following is the result by running MH algorithm. We can see that the two methods give us very close results. The code for MH algorithm for this example is as follows.

set.seed(1001)
theta<-rbeta(10, 1, 5)
theta2<-theta
set.seed(1002)
y <- rbinom(10, 100, theta)
n <- 100
J <- 10

##############################################################

log.prior <- function(alpha,beta) {
  {-2.5}*log(alpha + beta)
}

draw.thetas <- function(alpha,beta) {
  return(rbeta(J,alpha+y,beta+n-y))
}

draw.alpha <- function(alpha,beta,theta,prop.sd) {
  alpha.star <- rnorm(1,alpha,prop.sd)
  num <- J*(lgamma(alpha.star+beta) - lgamma(alpha.star)) +
    alpha.star*sum(log(theta)) + log(dgamma(alpha.star, 1, 1))
  den <- J*(lgamma(alpha+beta)      - lgamma(alpha)) +
    alpha     *sum(log(theta)) + log(dgamma(alpha, 1, 1))
# print(c(alpha,alpha.star,num,den))
  acc <- ifelse((log(runif(1))<=num - den)&&(alpha.star>0),1,0)
  alpha.acc <<- alpha.acc + acc
  return(ifelse(acc,alpha.star,alpha))
}

draw.beta <- function(alpha,beta,theta,prop.sd) {
  beta.star <- rnorm(1,beta,prop.sd)
  num <- J*(lgamma(alpha+beta.star) - lgamma(beta.star)) +
    beta.star*sum(log(1-theta)) + log(dgamma(beta.star,1, 1))
  den <- J*(lgamma(alpha+beta)      - lgamma(beta)) +
    beta     *sum(log(1-theta)) + log(dgamma(beta, 1, 1))
# print(c(beta,beta.star,num,den))
  acc <- ifelse((log(runif(1))<=num - den)&&(beta.star>0),1,0)
  beta.acc <<- beta.acc + acc

  return(ifelse(acc,beta.star,beta))
}

################################################################  

B <- 100
M <- 10000

MM <- B + M

alpha <- numeric(MM)
beta <- alpha
theta <- matrix(NA,nrow=MM,ncol=J)

# Metropolis tuning parameters
alpha.prop.sd <-  0.25
beta.prop.sd <-   3

# Initial values for the chain
alpha[1] <- 1
beta[1] <- 1
theta[1,] <- draw.thetas(alpha[1],beta[1]) # or theta[1,] <- (y+.5/(n+.5)

# Monitor acceptance frequency
alpha.acc <- 0
beta.acc <- 0

# MCMC simulation
for (m in 2:MM) {
  alpha[m] <- draw.alpha(alpha[m-1],beta[m-1],theta[m-1,],alpha.prop.sd)
  beta[m] <- draw.beta(alpha[m],beta[m-1],theta[m-1,],beta.prop.sd)
  theta[m,] <- draw.thetas(alpha[m],beta[m])
}

good <- (B+1):MM

alpha.mcmc <- alpha[good]
beta.mcmc <- beta[good]
theta.mcmc <- theta[good,]

par(mfrow=c(2,2))
plot(alpha.mcmc,type="l", main = "Trace of alpha")
plot(density(alpha.mcmc), type = "l", main = "density of alpha")
plot(beta.mcmc,type="l", main = "Trace of beta")
plot(density(beta.mcmc), type = "l", main = "density of beta")