Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

laplace method #800

Merged
merged 28 commits into from
Sep 25, 2023
Merged

laplace method #800

merged 28 commits into from
Sep 25, 2023

Conversation

jgabry
Copy link
Member

@jgabry jgabry commented Jul 30, 2023

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:

@jgabry
Copy link
Member Author

jgabry commented Jul 30, 2023

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
Copy link

codecov-commenter commented Jul 30, 2023

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 jgabry marked this pull request as ready for review July 31, 2023 17:00
@jgabry
Copy link
Member Author

jgabry commented Jul 31, 2023

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
Copy link
Member Author

jgabry commented Jul 31, 2023

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

@avehtari
Copy link
Contributor

avehtari commented Aug 3, 2023

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
Copy link
Member Author

jgabry commented Aug 3, 2023

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
Copy link
Contributor

avehtari commented Aug 4, 2023

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
Copy link
Contributor

avehtari commented Aug 4, 2023

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
Copy link
Contributor

avehtari commented Aug 4, 2023

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
Copy link
Member Author

jgabry commented Aug 4, 2023

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
Copy link
Member Author

jgabry commented Aug 4, 2023

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
Copy link
Member Author

jgabry commented Aug 4, 2023

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.

Copy link
Contributor

@avehtari avehtari left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good. Added one suggestion

R/fit.R Outdated Show resolved Hide resolved
@jgabry jgabry changed the title Add laplace method laplace method Sep 14, 2023
@jgabry
Copy link
Member Author

jgabry commented Sep 19, 2023

@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
Copy link
Collaborator

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

@jgabry
Copy link
Member Author

jgabry commented Sep 20, 2023

No worries, and thanks!

@andrjohns
Copy link
Collaborator

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
Copy link
Member Author

jgabry commented Sep 25, 2023

Ok thanks for working on this @andrjohns

@andrjohns
Copy link
Collaborator

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

@jgabry
Copy link
Member Author

jgabry commented Sep 25, 2023

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

@jgabry jgabry merged commit a63e418 into master Sep 25, 2023
12 checks passed
@jgabry jgabry deleted the laplace-sample branch September 25, 2023 18:39
@jgabry jgabry mentioned this pull request Sep 25, 2023
2 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add method laplace_sample
4 participants