mcMAM: Optimization of the minimum (maximum) age model (using a Markov chain Monte Carlo method)

Description

Sampling from the joint-likelihood function of the minimum (maximum) age model using a Markov chain Monte Carlo (MCMC) method.

Usage

mcMAM(EDdata, ncomp = -1, addsigma = 0, iflog = TRUE, 
      nsim = 50000, inis = list(), control.args = list())

Value

Return an invisible list of S3 class object "mcAgeModels". See mcFMM for details.

Arguments

EDdata: matrix(required): a two-column matrix (i.e., equivalent dose values and
associated standard errors)
ncomp: integer(with default): number of components, ncomp=-1, ncomp=-2, ncomp=-3, or ncomp=-4 indicate fitting the "MAM3", "MAM4", "MXAM3", and "MXAM4", respectively
addsigma: numeric(with default): additional uncertainty, i.e., the sigmab value
iflog: logical(with default): transform equivalent dose values to log-scale or not
nsim: integer(with default): deseried number of iterations
inis: list(with default): initial state of parameters.
Example: inis=list(p=0.1,gamma=20,sigma=0.3) when ncomp=-1
control.args: list(with default): arguments used by the Slice Sampling algorithm, see function mcFMM for details

References

Galbraith RF, Roberts RG, Laslett GM, Yoshida H, Olley JM, 1999. Optical dating of single grains of quartz from Jinmium rock shelter, northern Australia. Part I: experimental design and statistical models. Archaeometry, 41(2): 339-364.

Neal RM, 2003. "Slice sampling" (with discussion). Annals of Statistics, 31(3): 705-767. Software is freely available at https://glizen.com/radfordneal/slice.software.html.

Peng J, Dong ZB, Han FQ, 2016. Application of slice sampling method for optimizing OSL age models used for equivalent dose determination. Progress in Geography, 35(1): 78-88. (In Chinese).

Examples

Run this code

  # Not run.
  # data(EDdata)
  # Construct a MCMC chain for MAM3.
  # obj<-mcMAM(EDdata$al3,ncomp=-1,addsigma=0.1,nsim=5000)
  # reportMC(obj,burn=1e3,thin=2)
  #
  # The convergence of the simulations may be diagnosed with 
  # the Gelman and Rubin's convergence diagnostic.
  # library(coda)   # Only if package "coda" has been installed.
  # args<-list(nstart=50)
  # inis1<-list(p=0.01,gamma=26,mu=104,sigma=0.01)
  # inis2<-list(p=0.99,gamma=100,mu=104,sigma=4.99)
  # obj1<-mcMAM(EDdata$al3,ncomp=-2,nsim=3000,inis=inis1,control.args=args)
  # obj2<-mcMAM(EDdata$al3,ncomp=-2,nsim=3000,inis=inis2,control.args=args)
  # chain1<-mcmc(obj1$chains)
  # chain2<-mcmc(obj2$chains)
  # chains<-mcmc.list(chain1,chain2)
  # gelman.plot(chains)