MCMCWithinGibbs: Wrapper for the MCMC

Usage Arguments Value Author(s) Examples

View source: R/MCMCtools.R

Usage

1
MCMCWithinGibbs(theta0, data, hyper, nbatch, family, step.size = NULL, keepNM = FALSE)

Arguments

theta0

Starting value for the Markov chain. Must be of length three for the binomial model, and length four for the multbinom or doublebinom models. The elements must also be named correctly. See example usage below.

data

The data

hyper

A list containing the relevant hyperparameter information. hyper$a.m and hyper$b.m are the parameters in the Beta(a,b) prior on the mortality rate. hyper$alpha.lambda and hyper$beta.lambda are the parameters in a Gamma(alpha, beta) prior on lambda. hyper$a.p and hyper$b.p are the parameters in the Beta(a,b) prior on the probability p. hyper$mu.psi and hyper$sd.psi are the mean and standard deviation parameters for the N(mu, sd) prior on psi.

nbatch

Number of iterations in the MCMC

family

Either binomial, multbinom, or doublebinom depending on the model we wish to analyse.

step.size

A vector of length 2 giving the step size for the Metropolis-Hastings algorithm. Only needed for the multbinom and doublebinom models. Needs to be name correctly (see example usage below)

keepNM

logical; if TRUE, the Markov chains for N and M are returned.

Value

chain

Output of the Markov chain

accep
acc.vec
N.chain

The Markov chain output for the latent N values

M.chain

The Markov chain output for the latent M values

Author(s)

Richard Wilkinson

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#demo(MCMC, package="precision")

data(GlegneriSecondary)
###########################################################
#### Define hyperparameters/priors
####  - the four built in secondary datasets have priors for d and lambda built in, which load when you load the data
####  If you are using a different datasets you'll need to run the commands
# hyper<-list()
## prior for the mortality rate is Gamma(alpha, beta)
# hyper$a.m <- 6  # replace the value with something sensible for your data
# hyper$b.m <-  10

## prior for average clutch size, lambda, is Gamma(alpha, beta)
# hyper$alpha.lambda <- 4
# hyper$beta.lambda  <- 1

## prior for p is beta(a, b)
hyper$a.p <-  1 ## 20   ## Hyper parameters for p's prior distribution
hyper$b.p <-  1##100

# prior for psi - assumed to be Gaussian
hyper$mu.psi <- 0  
hyper$sd.psi <- 1
##########################################
### MCMC parameters

nbatch <- 1*10^3   ## number of MCMC samples to generate - usually we'll require at least 10^5 iterations 
thin = FALSE    
burnin <- 10^4 
thinby <- 3  

#######################################################################################################
# BINOMIAL MODEL
#######################################################################################################

## Define a start point for the MCMC chain - important to name the parameter and the column of the NM matrix
b.theta0 <-c(10, 0.1, 0.5) 
names(b.theta0) <-c("lambda", "p", "mort")

b.mcmc.out <- MCMCWithinGibbs( theta0=b.theta0,  data=GlegneriSecondary, hyper=hyper, nbatch=nbatch, family="binomial", keepNM=TRUE)


## Thin the chain or not
if(thin){
  b.mcmc.out.t  <- ThinChain(b.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
  b.mcmc.out.t <- b.mcmc.out
}
rm(b.mcmc.out)  ## delete as not needed and takes a lot of memory


#######################################################################################################
# MULTIPLICATIVE BINOMIAL MODEL
#######################################################################################################

## Define the Metropolis-Hastings random walk step size
m.step.size<-c(0.3,0.2)   ## 0.3 and 0.2 aree 
names(m.step.size) <- c("p.logit", "psi")


m.theta0 <-c(10, 0.1, 0, 0.1) 
names(m.theta0) <-c("lambda", "p", "psi", "mort")

m.mcmc.out <- MCMCWithinGibbs( theta0=m.theta0,  data=GlegneriSecondary, hyper=hyper, nbatch=nbatch,  family="multbinom", step.size=m.step.size, keepNM=TRUE)

if(thin){
  m.mcmc.out.t  <- ThinChain(m.mcmc.out, thinby=thinby, burnin=burnin)
}
if(!thin){
  m.mcmc.out.t <- m.mcmc.out
}
rm(m.mcmc.out)   ## to save space

rich-d-wilkinson/precision documentation built on May 27, 2019, 7:41 a.m.