iterLap: Iterated Laplace Approximation

In iterLap: Approximate Probability Densities by Iterated Laplace Approximations

Iterated Laplace Approximation

Description

Iterated Laplace Approximation

Usage

iterLap(post, ..., GRobj = NULL, vectorized = FALSE, startVals = NULL,
        method = c("nlminb", "nlm", "Nelder-Mead", "BFGS"), control = NULL,
        nlcontrol = list())

Arguments

post	log-posterior density
...	additional arguments to log-posterior density
GRobj	object of class mixDist, for example resulting from a call to GRApprox
vectorized	Logical determining, whether post is vectorized
startVals	Starting values for GRApprox, when GRobj is not specified. Vector of starting values if dimension=1 otherwise matrix of starting values with the starting values in the rows
method	Type of optimizer to be used.
control	List with entries: gridSize Determines the size of the grid for each component delta Stopping criterion based on the maximum error on the grid maxDim Maximum number of components allowed (default 20) eps Stopping criterion based on normalization constant of approximation info How much information should be displayed during iterations: 0 - none, 1 - minimum information, 2 - maximum information
nlcontrol	Control list for the used optimizer.

Value

Produces an object of class mixDist: A list with entries
weights Vector of weights for individual components
means Matrix of component medians of components
sigmas List containing scaling matrices
eigenHess List containing eigen decompositions of scaling matrices
dets Vector of determinants of scaling matrix
sigmainv List containing inverse scaling matrices

Author(s)

Bjoern Bornkamp

References

Bornkamp, B. (2011). Approximating Probability Densities by Iterated Laplace Approximations, Journal of Computational and Graphical Statistics, 20(3), 656–669.

Examples

  #### banana example
  banana <- function(pars, b, sigma12){
    dim <- 10
    y <- c(pars[1], pars[2]+b*(pars[1]^2-sigma12), pars[3:dim])
    cc <- c(1/sqrt(sigma12), rep(1, dim-1))
    return(-0.5*sum((y*cc)^2))
  }

  ###############################################################
  ## first perform multi mode Laplace approximation
  start <- rbind(rep(0,10),rep(-1.5,10),rep(1.5,10))
  grObj <- GRApprox(banana, start, b = 0.03, sigma12 = 100)
  ## print mixDist object
  grObj
  ## summary method
  summary(grObj)
  ## importance sampling using the obtained mixDist object 
  ## using a mixture of t distributions with 10 degrees of freedom
  isObj <- IS(grObj, nSim=1000, df = 10, post=banana, b = 0.03,
              sigma12 = 100)
  ## effective sample size
  isObj$ESS
  ## independence Metropolis Hastings algorithm
  imObj <- IMH(grObj, nSim=1000, df = 10, post=banana, b = 0.03,
               sigma12 = 100)
  ## acceptance rate
  imObj$accept

  ###############################################################
  ## now use iterated Laplace approximation
  ## and use Laplace approximation above as starting point
  iL <- iterLap(banana, GRobj = grObj, b = 0.03, sigma12 = 100)
  isObj2 <- IS(iL, nSim=10000, df = 100, post=banana, b = 0.03,
               sigma12 = 100)
  ## residual resampling to obtain unweighted sample
  samples <- resample(1000, isObj2)
  ## plot samples in the first two dimensions
  plot(samples[,1], samples[,2], xlim=c(-40,40), ylim = c(-40,20))
  ## independence Metropolis algorithm
  imObj2 <- IMH(iL, nSim=1000, df = 10, post=banana, b = 0.03,
                sigma12 = 100)
  imObj2$accept
  plot(imObj2$samp[,1], imObj2$samp[,2], xlim=c(-40,40), ylim = c(-40,20))

  ## IMH and IS can exploit multiple cores, example for two cores
  ## Not run: 
  isObj3 <- IS(iL, nSim=10000, df = 100, post=banana, b = 0.03,
               sigma12 = 100, cores = 2)
  
## End(Not run)

iterLap documentation built on Oct. 1, 2023, 1:06 a.m.