dsaddle {esaddle}

Evaluating the Extended Empirical Saddlepoint (EES) density

Description

Gives a pointwise evaluation of the EES density (and optionally of its gradient) at one or more locations.

Usage

dsaddle(
  y,
  X,
  decay,
  deriv = FALSE,
  log = FALSE,
  normalize = FALSE,
  control = list(),
  multicore = !is.null(cluster),
  ncores = detectCores() - 1,
  cluster = NULL
)

Arguments

y

points at which the EES is evaluated (d dimensional vector) or an n by d matrix, each row indicating a different position.

X

n by d matrix containing the data.

decay

rate at which the EES falls back on a normal density approximation, fitted to X. It must be a positive number, and it is inversely proportional to the complexity of the fit. Setting it to Inf leads to a Gaussian fit.

deriv

If TRUE also the gradient of the log-saddlepoint density is returned.

log

If TRUE the log of the saddlepoint density is returned.

normalize

If TRUE the normalizing constant of the EES density will be computed. FALSE by default.

control

A list of control parameters with entries:

• method the method used to calculate the normalizing constant. Either "LAP" (laplace approximation) or "IS" (importance sampling).

• nNorm if control$method == "IS", this is the number of importance samples used.

• tol the tolerance used to assess the convergence of the solution to the saddlepoint equation. The default is 1e-6.

• maxit maximal number of iterations used to solve the saddlepoint equation. The default is 100;

• ml Relevant only if control$method=="IS". n random variables are generated from a Gaussian importance density with covariance matrix ml*cov(X). By default the inflation factor is ml=2.

multicore

if TRUE the empirical saddlepoint density at each row of y will be evaluated in parallel.

ncores

number of cores to be used.

cluster

an object of class c("SOCKcluster", "cluster"). This allowes the user to pass her own cluster, which will be used if multicore == TRUE. The user has to remember to stop the cluster.

Value

A list with entries:

• llk the value of the EES log-density at each location y;

• mix for each location y, the fraction of saddlepoint used: 1 means that only ESS is used and 0 means that only a Gaussian fit is used;

• iter for each location y, the number of iteration needed to solve the saddlepoint equation;

• lambda an n by d matrix, where the i-th row is the solution of the saddlepoint equation corresponding to the i-th row of y;

• grad the gradient of the log-density at y (optional);

• logNorm the estimated log normalizing constant (optional);

Author(s)

Matteo Fasiolo <[e-mail]> and Simon N. Wood.

References

Fasiolo, M., Wood, S. N., Hartig, F. and Bravington, M. V. (2016). An Extended Empirical Saddlepoint Approximation for Intractable Likelihoods. ArXiv http://arxiv.org/abs/1601.01849.

Examples

library(esaddle)

### Simple univariate example
set.seed(4141)
x <- rgamma(1000, 2, 1)

# Evaluating EES at several point
xSeq <- seq(-2, 8, length.out = 200)
tmp <- dsaddle(y = xSeq, X = x, decay = 0.05, log = TRUE)  # Un-normalized EES
tmp2 <- dsaddle(y = xSeq, X = x, decay = 0.05,             # EES normalized by importance sampling
                normalize = TRUE, control = list("method" = "IS", nNorm = 500), log = TRUE)

# Plotting true density, EES and normal approximation
plot(xSeq, exp(tmp$llk), type = 'l', ylab = "Density", xlab = "x")
lines(xSeq, dgamma(xSeq, 2, 1), col = 3)
lines(xSeq, dnorm(xSeq, mean(x), sd(x)), col = 2)
lines(xSeq, exp(tmp2$llk), col = 4)
suppressWarnings( rug(x) )
legend("topright", c("EES un-norm", "EES normalized", "Truth", "Gaussian"), 
        col = c(1, 4, 3, 2), lty = 1)