posteriorMCMC {BMAmevt}

MCMC sampler for parametric spectral measures

Description

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

Usage

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 "".

Value

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.

See Also

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

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)