Usage Arguments Value Author(s) Examples
1 | MCMCWithinGibbs(theta0, data, hyper, nbatch, family, step.size = NULL, keepNM = FALSE)
|
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. |
nbatch |
Number of iterations in the MCMC |
family |
Either |
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. |
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 |
Richard Wilkinson
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
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.