Laplace sampling

The laplace method produces a sample from a normal approximation centered at the mode of a distribution in the unconstrained space. If the mode is a maximum a posteriori (MAP) estimate, the samples provide an estimate of the mean and standard deviation of the posterior distribution. If the mode is a maximum likelihood estimate (MLE), the sample provides an estimate of the standard error of the likelihood. In general, the posterior mode in the unconstrained space doesn’t correspond to the mean (nor mode) in the constrained space, and thus the sample is needed to infer the mean as well as the standard deviation. (See this case study for a visual illustration.)

This is computationally inexpensive compared to exact Bayesian inference with MCMC. The goodness of this estimate depends on both the estimate of the mode and how much the true posterior in the unconstrained space resembles a Gaussian.

Configuration

This method takes several arguments:

  • mode - Input file of parameters values on the constrained scale. When Stan’s optimize method is used to estimate the modal values, the value of boolean argument jacobian should be false if optimize was run with default settings, i.e., the input is the MLE estimate; if optimize was run with argument jacobian=true, then the laplace method default setting, jacobian=true, should be used.

  • jacobian - Whether or not the Jacobian adjustment should be included in the gradient. The default value is true (include adjustment). (Note: in optimization, the default value is false, for historical reasons.)

  • draws - How many total draws to return. The default is \(1000\).

  • calculate_lp - Whether to calculate the log probability of the model at each draw. If this is false, the log_p__ column of the output will be entirely nan. The default value is true.

CSV output

The output file consists of the following pieces of information:

  • The full set of configuration options available for the laplace method is reported at the beginning of the output file as CSV comments.

  • Output columns log_p__ and log_q__, the unnormalized log density and the unnormalized density of the Laplace approximation, respectively. These can be used for diagnostics and importance sampling.

  • Output columns for all model parameters on the constrained scale.

Diagnostic file outputs

If requested with output diagnostic_file=, a JSON file will be created which contains the log density, the gradient, and the Hessian of the log density evaluated at the mode.

Example

To get an approximate estimate of the mode and standard deviation of the example Bernoulli model given the example dataset:

  • find the MAP estimate by running optimization with argument jacobian=true

  • run the Laplace estimator using the MAP estimate as the mode argument.

Because the default output file name from all methods is output.csv, a more informative name is used for the output of optimization. We run the commands from the CmdStan home directory. This results in a sample with mean 2.7 and standard deviation 0.12. In comparison, running the NUTS-HMC sampler results in mean 2.6 and standard deviation 0.12.

./examples/bernoulli/bernoulli optimize jacobian=1 \
  data file=examples/bernoulli/bernoulli.data.json \
  output file=bernoulli_optimize_lbfgs.csv random seed=1234


./examples/bernoulli/bernoulli laplace mode=bernoulli_optimize_lbfgs.csv \
 data file=examples/bernoulli/bernoulli.data.json random seed=1234

The header and first few data rows of the output sample are shown below.

# stan_version_major = 2
# stan_version_minor = 35
# stan_version_patch = 0
# model = bernoulli_model
# start_datetime = 2022-12-20 01:01:14 UTC
# method = laplace
#   laplace
#     mode = bernoulli_lbfgs.csv
#     jacobian = true (Default)
#     draws = 1000 (Default)
#     calculate_lp = true (default)
# id = 1 (Default)
# data
#   file = examples/bernoulli/bernoulli.data.json
# init = 2 (Default)
# random
#   seed = 875960551 (Default)
# output
#   file = output.csv (Default)
#   diagnostic_file =  (Default)
#   refresh = 100 (Default)
#   sig_figs = -1 (Default)
#   profile_file = profile.csv (Default)
# num_threads = 1 (Default)
# stanc_version = stanc3 v2.31.0-7-g20444266
# stancflags =
log_p__,log_q__,theta
-9.4562,-2.33997,0.0498545
-6.9144,-0.0117349,0.182898
-7.18171,-0.746034,0.376428
...
Back to top