laplace method

[jgabry]

Submission Checklist

• Run unit tests
• Declare copyright holder and agree to license (see below)

Summary

Closes #760
Builds off of PR #799

Implements new method CmdStanModel$laplace() with new fitted model class CmdStanLaplace.
Follows the design from @WardBrian in #760 (comment).

Copyright and Licensing

Please list the copyright holder for the work you are submitting
(this will be you or your assignee, such as a university or company):
Columbia University

By submitting this pull request, the copyright holder is agreeing to
license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[jgabry]

Example usage

file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
mod <- cmdstan_model(file)
mod$print()

stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))

# pass CmdStanMLE from optimize to 'mode' argument of laplace
fit_mode <- mod$optimize(data = stan_data, jacobian = TRUE)
fit_laplace <- mod$laplace(data = stan_data, mode = fit_mode)
fit_laplace$summary()
fit_laplace$draws("theta")

# can also pass CSV file to 'mode' argument
fit_laplace <- mod$laplace(data = stan_data, mode = fit_mode$output_files())
fit_laplace$summary()

# if mode isn't specified optimize is run internally first
# can pass arguments to optimize via opt_args
fit_laplace <- mod$laplace(data = stan_data, opt_args = list(iter = 200))
fit_laplace$summary()

[codecov-commenter]

Codecov Report

Merging #800 (8a9ae96) into master (17678d5) will increase coverage by 0.30%.
The diff coverage is 98.54%.

❗ Current head 8a9ae96 differs from pull request most recent head 5c2dbbd. Consider uploading reports for the commit 5c2dbbd to get more accurate results

@@            Coverage Diff             @@
##           master     #800      +/-   ##
==========================================
+ Coverage   88.19%   88.49%   +0.30%     
==========================================
  Files          12       12              
  Lines        4218     4347     +129     
==========================================
+ Hits         3720     3847     +127     
- Misses        498      500       +2     

Files Changed	Coverage Δ
R/example.R	97.56% <ø> (ø)
R/fit.R	96.06% <88.88%> (-0.16%)	⬇️
R/model.R	92.31% <98.41%> (+0.49%)	⬆️
R/args.R	97.87% <100.00%> (+0.11%)	⬆️
R/csv.R	96.62% <100.00%> (+0.19%)	⬆️
R/run.R	93.90% <100.00%> (+0.01%)	⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[jgabry]

I think this is ready for review. I've implemented everything that I can think of. Anything I need to do when adding a new method that I'm forgetting?

[jgabry]

Any ideas why everything is passing except with the WSL backend?

[avehtari]

I'm still in progress of testing, but one observation now. The laplace method has argument draws but the variational method has output_samples and furthermore the posterior package resample_draws has ndraws. Could we have some common argument name?

[jgabry]

I agree. Those are the names CmdStan uses, but I actually added the draws argument to variational in a recent commit to match Laplace. output_samples will still work but we can deprecate it eventually.

[avehtari]

I guess this should be done in CmdStan, but mentioning now here. It would be nice to be able to be able to parallelize the computation for the draws from the approximation similarly way as $sampling() can do. Often the optimization is fast, but also often obtaining 4000 draws and computing lp__ and lp_approx__ takes time that is similar to getting 1000 MCMC draws.

[avehtari]

Here's an example to do resample importance sampling

file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
mod <- cmdstan_model(file)
mod$print()

stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,0))
fit_mode <- mod$optimize(data = stan_data, jacobian = TRUE)
# 1000 draws from the approximation
fit_laplace <- mod$laplace(data = stan_data, mode = fit_mode, refresh=0)
draws_laplace <- fit_laplace$draws(variables="theta")
# Compute Pareto smoothed importance sampling weights
(psis_laplace <- psis(log_ratios=fit_laplace$draws("lp__")-fit_laplace$draws("lp_approx__"), r_eff=1))
# Resample 1000 draws using simple resampling without resampling
draws_psis <- resample_draws(draws_laplace,
                             weights = exp(psis_laplace$log_weights),
                             ndraws = 1000,
                             method = 'simple')

summarize_draws(draws_laplace, default_summary_measures())
summarize_draws(draws_psis, default_summary_measures())

library(patchwork)
p1 <- mcmc_areas(draws_laplace) + xlim(c(0,0.7))
p2 <- mcmc_areas(draws_psis) + xlim(c(0,0.7))
p1 / p2

[avehtari]

Is it possible to get draws without evaluating lp__ and lp_approx__? In the cases when I would trust the normal approximation, there is unnecessary but computationally potentially very costly evaluation of lp__.

[jgabry]

I guess this should be done in CmdStan, but mentioning now here. It would be nice to be able to be able to parallelize the computation for the draws from the approximation similarly way as $sampling() can do. Often the optimization is fast, but also often obtaining 4000 draws and computing lp__ and lp_approx__ takes time that is similar to getting 1000 MCMC draws.

I agree. But yeah I think this is more of a CmdStan issue than a CmdStanR issue.

[jgabry]

Is it possible to get draws without evaluating lp__ and lp_approx__? In the cases when I would trust the normal approximation, there is unnecessary but computationally potentially very costly evaluation of lp__.

Not by changing anything in CmdStanR unfortunately. All we're doing is providing access to what CmdStan has already computed and written to CSV, so this would need to happen in CmdStan.

[jgabry]

Here's an example to do resample importance sampling

file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
mod <- cmdstan_model(file)
mod$print()

stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,0))
fit_mode <- mod$optimize(data = stan_data, jacobian = TRUE)
# 1000 draws from the approximation
fit_laplace <- mod$laplace(data = stan_data, mode = fit_mode, refresh=0)
draws_laplace <- fit_laplace$draws(variables="theta")
# Compute Pareto smoothed importance sampling weights
(psis_laplace <- psis(log_ratios=fit_laplace$draws("lp__")-fit_laplace$draws("lp_approx__"), r_eff=1))
# Resample 1000 draws using simple resampling without resampling
draws_psis <- resample_draws(draws_laplace,
                             weights = exp(psis_laplace$log_weights),
                             ndraws = 1000,
                             method = 'simple')

summarize_draws(draws_laplace, default_summary_measures())
summarize_draws(draws_psis, default_summary_measures())

library(patchwork)
p1 <- mcmc_areas(draws_laplace) + xlim(c(0,0.7))
p2 <- mcmc_areas(draws_psis) + xlim(c(0,0.7))
p1 / p2

Cool, thanks! Maybe we should put this in the documentation somewhere.

[avehtari]

Looks good. Added one suggestion

[jgabry]

@andrjohns When you have time would you be able to check if you can figure out why this is failing on WSL only? Thank you!

[andrjohns]

Ack sorry for forgetting about this! I'll have time to dig in on Friday

[jgabry]

No worries, and thanks!

[andrjohns]

Just a heads up that I'm still working on this, I'm on an M1 Mac as my daily driver now so I'm just ironing out some kinks on getting WSL running in my windows VM

[jgabry]

Ok thanks for working on this @andrjohns

[andrjohns]

@jgabry ended up being a super minor issue of path handling, all good to go now!

[jgabry]

Awesome, thanks @andrjohns! Glad it turned out to be a simple fix. Everything is passing so I'll go ahead and merge.