简单介绍使用 JAGS 进行蒙特卡罗模拟

Introduction to Monte Carlo simulation using JAGS

This practical gives an introduction to the Monte Carlo method using JAGS (Just Another Gibbs Sampler). First we have to install some R packages that act as frontenda for JAGS:

• R2jags allows fitting JAGS models from within R (depend on rjags);
• coda provides functions for output analysis and diagnostics for MCMC;
• lattice provides functions for visualisation of an MCMC object;
• MCHCvis in an R packages used to visualise, manipulate, and summarise MCMC output.

install.packages('R2jags', dependencies = TRUE) # install coda too
# install,packages('coda', dependencies = TRUE)
install.packages('lattice', dependencies = TRUE)
install.packages('MCMvis', dependencies = TRUE)

library(R2jags)
library(MCMCvis)
library(coda)
library(lattice)

Introduce

The main steps of fitting a Baysian model using R2jags are the following:

1. Define the JAGS model (as a function);
2. Prepare the data , which should be passed on as a list. Note that this list should include other known variables not in the data, like the length N of the data vector, when this N is used in the JAGS model definition (e.g. in a for loop);
3. Identify the parameters whose posterior distribution are of interest.
4. Define starting values for JAGS. We need starting values for all the stochastic nodes (unknoen parameters, missing data) in each chain. Initial values should be passed on as a list. Note that JAGS runs without user-defined starting values, in which cases it generates its own initial values. But this is not recommended, especially for complex models;
5. Fit the model in JAGS.
   The most important arguments of the #jags# function includes:
   data (data list), inits (list of initial values),
   parameters.to.save (parameters of interest),
   n.chains (number of MC chains - each gives a sequence of simulated values),
   n.iter (number of iterations in the sampling process),
   n.burnin (number of initial samples to discard),
   model.file (the JAGS model definition);
6. Convergence diagnostics. Update if convergence hasn’t reached;
7. Once we are satisfied that the chains have converged, we can start the statistical inference. E.g. we can plot the posterior density, get point estimates and interval estimates by extracting the posterior mean and credible intervals.

To demonstrate the structure of a Baysian model fit using R2jags we will consider some examples with increasing complexity.
For reproducibility we should always set the seed, e.g. set.seed(45621)

Coin example

Suppose we want to known the probability of getting 8 or more heads when we toss a fair coin 10 times.
When we toss a fair coin 10 times, the number of heads Y follows a binomial distribution with parameters n=10 and p=0.5. We are interested in $P(Y>=8)$.
First we represent this model using the BUGS language. Recall that ~ is used to represent stochastic dependent, and <- is used to represent logical dependence.

# model definition
jags.model.coin<- function(){
    Y ~ dbin(0.5,10) # our data model
    P8 <- ifelse(Y>7, 1, 0) # the probability of interest
}

For each iteration a sample from $Binom(10,0.5)$ is drawn, and it’s checked whether this sample is strictly greater than 7 or not (since the distribution is discrete, this is equivalent to $>=8$). For example, a run of length five could give us

$$Y=(5,5,8,4,9)$$
$$P8=(0,0,1,0,1)$$

In this simple example we have no data, and since this is a Monte Carlo sampling froma known distribution, there is no real need to provide initial va