#===============================================================================================================
# Load packages
#===============================================================================================================
sapply(c("ggplot2"
, "rjags"
, "R2jags"
, "hdrcde"
, "runjags"
, "mcmcplots"
, "KernSmooth"), library, character.only = TRUE)
#===============================================================================================================
# Model file
#===============================================================================================================
cat("
model {
# Likelihood
for (i in 1:N) {
x[i] ~ dbin(theta, n)
}
# Prior
n ~ dpois(mu)
lambda ~ dgamma(0.005, 0.005)
# lambda ~ dunif(0, 1000)
mu <- lambda/theta
theta ~ dunif(0, 1)
}
", file="jags_model_binomial.txt")
#===============================================================================================================
# Data
#===============================================================================================================
data.list <- list(x = c(53, 57, 66, 67, 72, NA), N = 6) # Waterbuck example from Raftery (1988)
#===============================================================================================================
# Inits
#===============================================================================================================
jags.inits <- function() {
list(
n = sample(max(data.list$x, na.rm = TRUE):1000, size = 1)
, theta = runif(1, 0, 1)
, lambda = runif(1, 1, 10)
# , cauchy = runif(1, 1, 1000)
# , mu = runif(1, 0, 5)
)
}
#===============================================================================================================
# Run the chains
#===============================================================================================================
# Parameters to store
params <- c("n"
, "theta"
, "lambda"
, "mu"
, paste("x[", which(is.na(data.list[["x"]])), "]", sep = "")
)
# MCMC settings
niter <- 500000 # number of iterations
nburn <- 20000 # number of iterations to discard (the burn-in-period)
nchains <- 5 # number of chains
# Run JAGS
out <- jags(
data = data.list
, parameters.to.save = params
, model.file = "jags_model_binomial.txt"
, n.chains = nchains
, n.iter = niter
, n.burnin = nburn
, n.thin = 50
, inits = jags.inits
, progress.bar = "text")
# Load packages
#===============================================================================================================
sapply(c("ggplot2"
, "rjags"
, "R2jags"
, "hdrcde"
, "runjags"
, "mcmcplots"
, "KernSmooth"), library, character.only = TRUE)
#===============================================================================================================
# Model file
#===============================================================================================================
cat("
model {
# Likelihood
for (i in 1:N) {
x[i] ~ dbin(theta, n)
}
# Prior
n ~ dpois(mu)
lambda ~ dgamma(0.005, 0.005)
# lambda ~ dunif(0, 1000)
mu <- lambda/theta
theta ~ dunif(0, 1)
}
", file="jags_model_binomial.txt")
#===============================================================================================================
# Data
#===============================================================================================================
data.list <- list(x = c(53, 57, 66, 67, 72, NA), N = 6) # Waterbuck example from Raftery (1988)
#===============================================================================================================
# Inits
#===============================================================================================================
jags.inits <- function() {
list(
n = sample(max(data.list$x, na.rm = TRUE):1000, size = 1)
, theta = runif(1, 0, 1)
, lambda = runif(1, 1, 10)
# , cauchy = runif(1, 1, 1000)
# , mu = runif(1, 0, 5)
)
}
#===============================================================================================================
# Run the chains
#===============================================================================================================
# Parameters to store
params <- c("n"
, "theta"
, "lambda"
, "mu"
, paste("x[", which(is.na(data.list[["x"]])), "]", sep = "")
)
# MCMC settings
niter <- 500000 # number of iterations
nburn <- 20000 # number of iterations to discard (the burn-in-period)
nchains <- 5 # number of chains
# Run JAGS
out <- jags(
data = data.list
, parameters.to.save = params
, model.file = "jags_model_binomial.txt"
, n.chains = nchains
, n.iter = niter
, n.burnin = nburn
, n.thin = 50
, inits = jags.inits
, progress.bar = "text")
#===============================================================================================================
# Inspect results
#===============================================================================================================
print(out
, digits = 2
, intervals = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975))
# Inspect results
#===============================================================================================================
print(out
, digits = 2
, intervals = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975))
Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 50
n.sims = 48000 iterations saved
mu.vect sd.vect 2.5% 10% 25% 50% 75% 90% 97.5% Rhat n.eff
lambda 62.90 5.18 53.09 56.47 59.45 62.74 66.19 69.49 73.49 1 48000
mu 521.28 968.41 92.31 113.02 148.00 232.87 467.10 1058.17 3014.82 1 1600
n 521.73 968.54 95.00 114.00 148.00 233.00 467.00 1060.10 3028.00 1 1600
theta 0.29 0.18 0.02 0.06 0.13 0.27 0.42 0.55 0.66 1 1600
x[6] 63.03 7.33 49.00 54.00 58.00 63.00 68.00 72.00 78.00 1 36000
deviance 34.88 1.53 33.63 33.70 33.85 34.34 35.34 36.81 39.07 1 48000
Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 50
n.sims = 48000 iterations saved
mu.vect sd.vect 2.5% 10% 25% 50% 75% 90% 97.5% Rhat n.eff
lambda 62.90 5.18 53.09 56.47 59.45 62.74 66.19 69.49 73.49 1 48000
mu 521.28 968.41 92.31 113.02 148.00 232.87 467.10 1058.17 3014.82 1 1600
n 521.73 968.54 95.00 114.00 148.00 233.00 467.00 1060.10 3028.00 1 1600
theta 0.29 0.18 0.02 0.06 0.13 0.27 0.42 0.55 0.66 1 1600
x[6] 63.03 7.33 49.00 54.00 58.00 63.00 68.00 72.00 78.00 1 36000
deviance 34.88 1.53 33.63 33.70 33.85 34.34 35.34 36.81 39.07 1 48000
The posterior mean of $N$$N$ is $522$$522$ and the posterior median is $233$$233$. I calculated the posterior mode of $N$$N$ on the log-scale and back-transformed the estimate:
jagsfit.mcmc <- as.mcmc(out)
jagsfit.mcmc <- combine.mcmc(jagsfit.mcmc)
hpd.80 <- hdr.den(log(as.vector(jagsfit.mcmc[, "n"])), prob = c(80), den = bkde(log(as.vector(jagsfit.mcmc[, "n"])), gridsize = 10000))
exp(hpd.80$mode)
[1] 149.8161
jagsfit.mcmc <- as.mcmc(out)
jagsfit.mcmc <- combine.mcmc(jagsfit.mcmc)
hpd.80 <- hdr.den(log(as.vector(jagsfit.mcmc[, "n"])), prob = c(80), den = bkde(log(as.vector(jagsfit.mcmc[, "n"])), gridsize = 10000))
exp(hpd.80$mode)
[1] 149.8161
And the 80%-HPD of $N$$N$ is:
(hpd.ints <- HPDinterval(jagsfit.mcmc, prob = c(0.8)))
lower upper
deviance 33.61011007 35.677810
lambda 56.08842502 69.089507
mu 72.42307587 580.027182
n 78.00000000 578.000000
theta 0.01026193 0.465714
x[6] 53.00000000 71.000000
(hpd.ints <- HPDinterval(jagsfit.mcmc, prob = c(0.8)))
lower upper
deviance 33.61011007 35.677810
lambda 56.08842502 69.089507
mu 72.42307587 580.027182
n 78.00000000 578.000000
theta 0.01026193 0.465714
x[6] 53.00000000 71.000000
