posteriorMCMC function

MCMC sampler for parametric spectral measures

Generates a posterior parameters sample, and computes the posterior mean and component-wise variance on-line.

posteriorMCMC(
  prior = function(type = c("r", "d"), n, par, Hpar, log, dimData) {
     NULL
 },
  proposal = function(type = c("r", "d"), cur.par, prop.par, MCpar, log) {
     NULL
 },
  likelihood = function(x, par, log, vectorial) {
     NULL
 },
  Nsim,
  dat,
  Hpar,
  MCpar,
  Nbin = 0,
  par.start = NULL,
  show.progress = floor(seq(1, Nsim, length.out = 20)),
  seed = NULL,
  kind = "Mersenne-Twister",
  save = FALSE,
  class = NULL,
  name.save = NULL,
  save.directory = "~",
  name.dat = "",
  name.model = ""
)

Arguments

• prior: The prior distribution: of type

  function(type=c("r","d"), n ,par, Hpar, log, dimData ), where dimData is the dimension of the sample space (e.g., for the two-dimensional simplex (triangle), dimData=3. Should return either a matrix with n rows containing a random parameter sample generated under the prior (if type == "d"), or the density of the parameter par (the logarithm of the density if log==TRUE. See prior.pb and prior.nl for templates.

• proposal: The proposal function: of type

  function(type = c("r","d"), cur.par, prop.par, MCpar, log). Should return the (logarithm of) the proposal density for the move cur.par --> prop.par if type == "d". If type =="r", proposal must return a candidate parameter, depending on cur.par, as a vector. See proposal.pb or proposal.nl

  for templates.

• likelihood: The likelihood function. Should be of type

  function(x, par, log, vectorial), where log and vectorial are logical flags indicating respectively if the result is to be returned on the log-scale, and if the value is a vector of length nrow(x) or a single number (the likelihood, or the log-likelihood, for the data set x). See dpairbeta or dnestlog

  for templates.

• Nsim: Total number of iterations to perform.

• dat: An angular data set, e.g., constructed by cons.angular.dat: A matrix which rows are the Cartesian coordinates of points on the unit simplex (summing to one).

• Hpar: A list containing Hyper-parameters to be passed to prior.

• MCpar: A list containing MCMC tuning parameters to be passed to proposal.

• Nbin: Length of the burn-in period.

• par.start: Starting point for the MCMC sampler.

• show.progress: An vector of integers containing the times (iteration numbers) at which a message showing progression will be printed on the standard output.

• seed: The seed to be set via

  set.seed.

• kind: The kind of random numbers generator. Default to "Mersenne-Twister". See set.seed for details.

• save: Logical. Should the result be saved ?

• class: Optional character string: additional class attribute to be assigned to the result. A predefined class "PBNLpostsample" exists, for which a method performing convergence diagnostics is defined (see diagnose )

• name.save: A character string giving the name under which the result is to be saved. If NULL (default), nothing is saved. Otherwise, the result is saved in file

  paste(save.directory,"/", name.save,".rda",sep=""). A "log" list is also saved, named

  paste(name.save, ".log", sep=""), in file

  paste(save.directory,"/", name.log,".rda",sep="").

• save.directory: A character string giving the directory where the result is to be saved (without trailing slash).

• name.dat: A character string naming the data set used for inference. Default to "".

• name.model: A character string naming the model. Default to "".

Returns

A list made of

• stored.vals: A (Nsim-Nbin)*d matrix, where d

  is the dimension of the parameter space.

• llh A vector of size (Nsim-Nbin) containing the loglikelihoods evaluated at each parameter of the posterior sample.

• lprior A vector of size (Nsim-Nbin) containing the logarithm of the prior densities evaluated at each parameter of the posterior sample.

• elapsed: The time elapsed, as given by proc.time between the start and the end of the run.

• Nsim: The same as the passed argument

• Nbin: idem.

• n.accept: The total number of accepted proposals.

• n.accept.kept: The number of accepted proposals after the burn-in period.

• emp.mean The estimated posterior parameters mean

• emp.sd The empirical posterior sample standard deviation.

Examples

data(Leeds)
data(pb.Hpar)
data(pb.MCpar)
postsample1 <- posteriorMCMC(Nsim=1e+3,Nbin=500,
         dat= Leeds,
         prior = prior.pb,
         proposal = proposal.pb,
         likelihood = dpairbeta,
         Hpar=pb.Hpar,
         MCpar=pb.MCpar)

dim(postsample1[[1]])
postsample1[-1]

## Not run:

## a more realistic one:

postsample2 <- posteriorMCMC(Nsim=50e+3,Nbin=15e+3,
         dat= Leeds,
         prior = prior.pb,
         proposal = proposal.pb,
         likelihood = dpairbeta,
         Hpar=pb.Hpar,
         MCpar=pb.MCpar)
dim(postsample2[[1]])
postsample2[-1]
## End(Not run)

See Also

posteriorMCMC.pb, posteriorMCMC.pb for specific uses in the PB and the NL models.