Dh$?v@eG4@1AjFOC=>Gw_V%!U>d*BayZJC7ctoN+x08;!ikkCYpf3F@^
zLgY%|5q^|Pc}2+e_kc2>yHM;`(Hb`612X3iGxW`9+H@yd@!3|f81s6 This vignette demonstrates how to use the OpenCL capabilities of
-CmdStan with CmdStanR. The functionality described in this vignette
-requires CmdStan 2.26.1 or newer. As of version 2.26.1, users can expect speedups with OpenCL when
-using vectorized probability distribution functions (functions with the
- The actual speedup for a model will depend on the particular
- In order to establish the bottlenecks in your model we recommend
-using profiling,
-which was introduced in Stan version 2.26.0. OpenCL is supported on most modern CPUs and GPUs. In order to use
-OpenCL in CmdStanR, an OpenCL runtime for the target device must be
-installed. A guide for the most common devices is available in the
-CmdStan manual’s chapter
-on parallelization. In case of using Windows, CmdStan requires the
- By default, models in CmdStanR are compiled without OpenCL
-support. Once OpenCL support is enabled, a CmdStan model will make use
-of OpenCL if the functions in the model support it. Technically no
-changes to a model are required to support OpenCL since the choice of
-using OpenCL is handled by the compiler, but it can still be useful to
-rewrite a model to be more OpenCL-friendly by using vectorization as
-much as possible when using probability distributions. Consider a simple logistic regression with parameters
- Some fake data will be useful to run this model: In this model, most of the computation will be handled by the
- To build the model with OpenCL support, add
- Running models with OpenCL requires specifying the OpenCL platform
-and device on which to run the model (there can be multiple). If the
-system has one GPU and no OpenCL CPU runtime, the platform and device
-IDs of the GPU are typically both On an Ubuntu system with both CPU and GPU OpenCL support,
- On this system the GPU is platform ID 0 and device ID 0, while the
-CPU is platform ID 1, device ID 0. These can be specified with the
- We’ll also run a version without OpenCL and compare the run
-times. The speedup of the OpenCL model is: This speedup will be determined by the particular GPU/CPU used, the
-input problem sizes (data as well as parameters) and if the model uses
-functions that can be run on the GPU or other OpenCL devices. This vignette is intended to be read after the Getting
-started with CmdStanR vignette. Please read that first for
-important background. In this document we provide additional details
-about compiling models, passing in data, and how CmdStan output is saved
-and read back into R. We will only use the The Subsequently, if you create a Internally, It is also possible to delay compilation when creating the
- If you are using CmdStan version 2.24 or later and CmdStanR version
-0.2.1 or later, you can run a pedantic check for your model. CmdStanR
-will always check that your Stan program does not contain any invalid
-syntax but with pedantic mode enabled the check will also warn you about
-other potential issues in your model, for example: For the latest information on the checks performed in pedantic mode
-see the Pedantic
-mode chapter in the Stan Reference Manual. Pedantic mode is available when compiling the model or when using the
-separate To turn on pedantic mode at compile time you can set
- To turn on pedantic mode separately from compilation use the
- Using If using CmdStan 2.27 or newer, you can obtain the names, types and
-dimensions of the data, parameters, transformed parameters and generated
-quantities variables of a Stan model using the The Each variable is represented as a list containing the type
-information (currently limited to By default, the executable is created in the same directory as the
-file containing the Stan program. You can also specify a different
-location with the There are three data formats that CmdStanR allows when fitting a
-model: Like the RStan interface, CmdStanR accepts a named list of R objects
-where the names correspond to variables declared in the data block of
-the Stan program. In the Bernoulli model the data is Because CmdStan doesn’t accept lists of R objects, CmdStanR will
-first write the data to a temporary JSON file using
- If you already have your data in a JSON file you can just pass that
-file directly to CmdStanR instead of using a list of R objects. For
-example, we could pass in the JSON file we created above using
- Finally, it is also possible to use the R dump file format. This is
-not recommended because CmdStan can process JSON faster than R
-dump, but CmdStanR allows it because CmdStan will accept files created
-by When fitting a model, the default behavior is to write the output
-from CmdStan to CSV files in a temporary directory. These files will be lost if you end your R session or if you remove
-the To save these files to a non-temporary location there are two
-options. You can either specify the With the exception of some diagnostic information, the CSV files are
-not read into R until their contents are requested by calling a method
-that requires them (e.g., After we call a method that requires the draws then if we reexamine
-the structure of the object we will see that the For models with many parameters, transformed parameters, or generated
-quantities, if only some are requested (e.g., by specifying the
- Internally, the If you need to manually create fitted model objects from CmdStan CSV
-files use This is pointless in our case since we have the original
- If CmdStanR does not yet provide a special method for processing these
-files but they can be read into R using R’s standard CSV reading
-functions. The column Our model has a single parameter CmdStanR can of course be used for developing other packages that
-require compiling and running Stan models as well as using new or custom
-Stan features available through CmdStan. You may compile a Stan model at runtime (e.g. just before sampling),
-or you may compile all the models inside the package file system in
-advance at installation time. The latter avoids compilations at runtime,
-which matters in centrally managed R installations where users should
-not compile their own software. To pre-compile all the models in a package, you may create top-level
-scripts When developing or testing new features it might be useful to have
-more information on how CmdStan is called internally and to see more
-information printed when compiling or running models. This can be
-enabled for an entire R session by setting the option
- CmdStanR (Command Stan R) is a lightweight interface to Stan for R users that provides an
-alternative to the traditional RStan interface. See the Comparison with RStan section
-later in this vignette for more details on how the two interfaces
-differ. Using CmdStanR requires installing the cmdstanr R
-package and also CmdStan, the command line interface to Stan. First we
-install the R package by running the following command in R. We can now load the package like any other R package. We’ll also load
-the bayesplot and posterior packages
-to use later in examples. CmdStanR requires a working installation of CmdStan,
-the shell interface to Stan. If you don’t have CmdStan installed then
-CmdStanR can install it for you, assuming you have a suitable C++
-toolchain. The requirements are described in the CmdStan Guide: To double check that your toolchain is set up properly you can call
-the If your toolchain is configured correctly then CmdStan can be
-installed by calling the Before CmdStanR can be used it needs to know where the CmdStan
-installation is located. When the package is loaded it tries to help
-automate this to avoid having to manually set the path every
-session: If the environment variable If no environment variable is found when loaded but any directory
-in the form If neither of these applies (or you want to subsequently change the
-path) you can use the To check the path to the CmdStan installation and the CmdStan version
-number you can use The The object The Stan program can be printed using the The path to the compiled executable is returned by the
- The There are many more arguments that can be passed to the
- The The The To convert an existing draws object to a different format use the
- In general, converting to a different draws format in this way will
-be slower than just setting the appropriate format initially in the call
-to the The vignette Working with
-Posteriors has more details on posterior draws, including how to
-reproduce the structured output RStan users are accustomed to getting
-from The The We see the number of divergences for each of the four chains, the
-number of times the maximum treedepth was hit for each chain, and the
-E-BFMI for each chain. In this case there were no warnings, so in order to demonstrate the
-warning messages we’ll use one of the CmdStanR example models that
-suffers from divergences. After fitting there is a warning about divergences. We can also
-regenerate this warning message later using
- CmdStanR also supports running Stan’s optimization algorithms and its
-algorithms for variational approximation of full Bayesian inference.
-These are run via the We can find the (penalized) maximum likelihood estimate (MLE) using
- Here’s a plot comparing the penalized MLE to the posterior
-distribution of For optimization, by default the mode is calculated without the
-Jacobian adjustment for constrained variables, which shifts the mode due
-to the change of variables. To include the Jacobian adjustment and
-obtain a maximum a posteriori (MAP) estimate set
- The Here we pass in the We can run Stan’s experimental Automatic Differentiation Variational
-Inference (ADVI) using the Stan version 2.33 introduced a new variational method called
-Pathfinder, which is intended to be faster and more stable than ADVI.
-For details on how Pathfinder works see the section in the CmdStan
-User’s Guide. Pathfinder is run using the Let’s extract the draws, make the same plot we made after running the
-other algorithms, and compare them all. approximation, and compare them
-all. In this simple example the distributions are quite similar, but
-this will not always be the case for more challenging problems. For more details on the The But if your model object is large, then Storage is even faster if you discard results you do not need to
-save. The following example saves only posterior draws and discards
-sampler diagnostics, user-specified initial values, and profiling
-data. See the vignette How
-does CmdStanR work? for more information about the composition
-of CmdStanR objects. The RStan interface (rstan package) is
-an in-memory interface to Stan and relies on R packages like
-Rcpp and inline to call C++ code from
-R. On the other hand, the CmdStanR interface does not directly call any
-C++ code from R, instead relying on the CmdStan interface behind the
-scenes for compilation, running algorithms, and writing results to
-output files. Allows other developers to distribute R packages with
-pre-compiled Stan programs (like rstanarm) on
-CRAN. (Note: As of 2023, this can mostly be achieved with CmdStanR as
-well. See Developing
-using CmdStanR.) Avoids use of R6 classes, which may result in more familiar
-syntax for many R users. CRAN binaries available for Mac and Windows. Compatible with latest versions of Stan. Keeping up with Stan
-releases is complicated for RStan, often requiring non-trivial changes
-to the rstan package and new CRAN releases of both
-rstan and StanHeaders. With CmdStanR
-the latest improvements in Stan will be available from R immediately
-after updating CmdStan using
- Running Stan via external processes results in fewer unexpected
-crashes, especially in RStudio. Less memory overhead. More permissive license. RStan uses the GPL-3 license while the
-license for CmdStanR is BSD-3, which is a bit more permissive and is the
-same license used for CmdStan and the Stan C++ source code. There are additional vignettes available that discuss other aspects
-of using CmdStanR. These can be found online at the CmdStanR
-website: To ask a question please post on the Stan forums: To report a bug, suggest a feature (including additions to these
-vignettes), or to start contributing to CmdStanR development (new
-contributors welcome!) please open an issue on GitHub:Running Stan on the GPU with OpenCL
- Rok Češnovar
-and Jonah Gabry
-
-
- Source: vignettes/articles-online-only/opencl.Rmd
- opencl.RmdIntroduction
-
-_lpdf or _lpmf suffix) and when the input
-variables contain at least 20,000 elements.lpdf/lpmf functions used and whether the
-lpdf/lpmf functions are the bottlenecks of the model. The
-more computationally complex the function is, the larger the expected
-speedup. The biggest speedups are expected when using the specialized
-GLM functions.OpenCL runtime
-
-OpenCL.lib to compile the model. If you experience issue
-compiling the model with OpenCL, run the below script and set
-path_to_opencl_lib to the path to the
-OpenCL.lib file on your system. If you are using CUDA, the
-path should be similar to the one listed here.
-
path_to_opencl_lib <- "C:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v11.5/lib/x64"
-cpp_options = list(
- paste0("LDFLAGS+= -L\"",path_to_opencl_lib,"\" -lOpenCL")
-)
-
-cmdstanr::cmdstan_make_local(cpp_options = cpp_options)
-cmdstanr::rebuild_cmdstan()Compiling a model with OpenCL
-
-alpha and beta, covariates X, and
-outcome y.
-data {
- int<lower=1> k;
- int<lower=0> n;
- matrix[n, k] X;
- array[n] int y;
-}
-parameters {
- vector[k] beta;
- real alpha;
-}
-model {
- target += std_normal_lpdf(beta);
- target += std_normal_lpdf(alpha);
- target += bernoulli_logit_glm_lpmf(y | X, alpha, beta);
-}
-
library(cmdstanr)
-
-# Generate some fake data
-n <- 250000
-k <- 20
-X <- matrix(rnorm(n * k), ncol = k)
-y <- rbinom(n, size = 1, prob = plogis(3 * X[,1] - 2 * X[,2] + 1))
-mdata <- list(k = k, n = n, y = y, X = X)bernoulli_logit_glm_lpmf function. Because this is a
-supported GPU function, it should be possible to accelerate it with
-OpenCL. Check here
-for a list of functions with OpenCL support.cpp_options = list(stan_opencl = TRUE) at the compilation
-step.
-
# Compile the model with STAN_OPENCL=TRUE
-mod_cl <- cmdstan_model("opencl-files/bernoulli_logit_glm.stan",
- cpp_options = list(stan_opencl = TRUE))Running models with OpenCL
-
-0, but the
-clinfo tool can be used to figure out for sure which
-devices are available.clinfo -l outputs:
-Platform #0: AMD Accelerated Parallel Processing
- `-- Device #0: gfx906+sram-ecc
-Platform #1: Intel(R) CPU Runtime for OpenCL(TM) Applications
- `-- Device #0: Intel(R) Core(TM) i7-4790 CPU @ 3.60GHzopencl_ids argument when running a model. The
-opencl_ids is supplied as a vector of length 2, where the
-first element is the platform ID and the second argument is the device
-ID.
-
fit_cl <- mod_cl$sample(data = mdata, chains = 4, parallel_chains = 4,
- opencl_ids = c(0, 0), refresh = 0)
-
# no OpenCL version
-mod <- cmdstan_model("opencl-files/bernoulli_logit_glm.stan", force_recompile = TRUE)
-fit_cpu <- mod$sample(data = mdata, chains = 4, parallel_chains = 4, refresh = 0)
-
fit_cpu$time()$total / fit_cl$time()$totalvignettes/children/settings-knitr.Rmd
- settings-knitr.RmdHow does CmdStanR work?
- Jonah Gabry and
-Rok Češnovar
-
-
- Source: vignettes/cmdstanr-internals.Rmd
- cmdstanr-internals.RmdIntroduction
-
-$sample() method in examples, but
-all model fitting methods work in a similar way under the hood.
-
library(cmdstanr)
-check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)Compilation
-
-Immediate compilation
-
-cmdstan_model() function creates a new
-CmdStanModel object. The CmdStanModel object
-stores the path to a Stan program as well as the path to a compiled
-executable.
-
stan_file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
-mod <- cmdstan_model(stan_file)
-mod$print()
-data {
- int<lower=0> N;
- array[N] int<lower=0, upper=1> y;
-}
-parameters {
- real<lower=0, upper=1> theta;
-}
-model {
- theta ~ beta(1, 1); // uniform prior on interval 0,1
- y ~ bernoulli(theta);
-}
-
mod$stan_file()
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.36.0/examples/bernoulli/bernoulli.stan"
-
mod$exe_file()
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.36.0/examples/bernoulli/bernoulli"CmdStanModel object from
-the same Stan file then compilation will be skipped (assuming the file
-hasn’t changed).
-
mod <- cmdstan_model(stan_file)cmdstan_model() first creates the
-CmdStanModel object from just the Stan file and then calls
-its $compile()
-method. Optional arguments to the $compile() method can be
-passed via ....
-
mod <- cmdstan_model(
- stan_file,
- force_recompile = TRUE,
- include_paths = "paths/to/directories/with/included/files",
- cpp_options = list(stan_threads = TRUE, STANC2 = TRUE)
-)Delayed compilation
-
-CmdStanModel object by specifying
-compile=FALSE and then later calling the
-$compile() method directly.
-
unlink(mod$exe_file())
-mod <- cmdstan_model(stan_file, compile = FALSE)
-mod$exe_file() # not yet created
-character(0)
-
mod$compile()
-mod$exe_file()
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.36.0/examples/bernoulli/bernoulli"Pedantic check
-
-
-
-~ symbols).$check_syntax() method of a
-CmdStanModel object. Internally this corresponds to setting
-the stanc (Stan transpiler) option
-warn-pedantic. Here we demonstrate pedantic mode with a
-Stan program that is syntactically correct but is missing a lower bound
-and a prior for a parameter.
-
stan_file_pedantic <- write_stan_file("
-data {
- int N;
- array[N] int y;
-}
-parameters {
- // should have <lower=0> but omitting to demonstrate pedantic mode
- real lambda;
-}
-model {
- y ~ poisson(lambda);
-}
-")pedantic=TRUE in the call to cmdstan_model()
-(or when calling the $compile() method directly if using
-the delayed compilation approach described above).mod_pedantic <- cmdstan_model(stan_file_pedantic, pedantic = TRUE)
-Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b6650db73c.stan', line 11, column 14: A
- poisson distribution is given parameter lambda as a rate parameter
- (argument 1), but lambda was not constrained to be strictly positive.
-Warning: The parameter lambda has no priors. This means either no prior is
- provided, or the prior(s) depend on data variables. In the later case,
- this may be a false positive.pedantic argument to the $check_syntax()
-method.mod_pedantic$check_syntax(pedantic = TRUE)
-Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model_febb1e69c7387a0e64cf13583e078104.stan', line 11, column 14: A
- poisson distribution is given parameter lambda as a rate parameter
- (argument 1), but lambda was not constrained to be strictly positive.
-Warning: The parameter lambda has no priors. This means either no prior is
- provided, or the prior(s) depend on data variables. In the later case,
- this may be a false positive.
-Stan program is syntactically correctpedantic=TRUE via the $check_syntax()
-method also has the advantage that it can be used even if the model
-hasn’t been compiled yet. This can be helpful because the pedantic and
-syntax checks themselves are much faster than compilation.file.remove(mod_pedantic$exe_file()) # delete compiled executable
-[1] TRUE
-rm(mod_pedantic)
-
-mod_pedantic <- cmdstan_model(stan_file_pedantic, compile = FALSE)
-mod_pedantic$check_syntax(pedantic = TRUE)
-Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model_febb1e69c7387a0e64cf13583e078104.stan', line 11, column 14: A
- poisson distribution is given parameter lambda as a rate parameter
- (argument 1), but lambda was not constrained to be strictly positive.
-Warning: The parameter lambda has no priors. This means either no prior is
- provided, or the prior(s) depend on data variables. In the later case,
- this may be a false positive.
-Stan program is syntactically correctStan model variables
-
-$variables()
-method of the CmdStanModel object.
-
stan_file_variables <- write_stan_file("
-data {
- int<lower=1> J;
- vector<lower=0>[J] sigma;
- vector[J] y;
-}
-parameters {
- real mu;
- real<lower=0> tau;
- vector[J] theta_raw;
-}
-transformed parameters {
- vector[J] theta = mu + tau * theta_raw;
-}
-model {
- target += normal_lpdf(tau | 0, 10);
- target += normal_lpdf(mu | 0, 10);
- target += normal_lpdf(theta_raw | 0, 1);
- target += normal_lpdf(y | theta, sigma);
-}
-")
-mod_v <- cmdstan_model(stan_file_variables)
-variables <- mod_v$variables()$variables() method returns a list with
-data, parameters,
-transformed_parameters and
-generated_quantities elements, each corresponding to
-variables in their respective block of the program. Transformed data
-variables are not listed as they are not used in the model’s input or
-output.
-
names(variables)
-[1] "parameters" "included_files" "data"
-[4] "transformed_parameters" "generated_quantities"
-
names(variables$data)
-[1] "J" "sigma" "y"
-
names(variables$parameters)
-[1] "mu" "tau" "theta_raw"
-
names(variables$transformed_parameters)
-[1] "theta"
-
names(variables$generated_quantities)
-character(0)real or int)
-and the number of dimensions.
-
variables$data$J
-$type
-[1] "int"
-
-$dimensions
-[1] 0
-
variables$data$sigma
-$type
-[1] "real"
-
-$dimensions
-[1] 1
-
variables$parameters$tau
-$type
-[1] "real"
-
-$dimensions
-[1] 0
-
variables$transformed_parameters$theta
-$type
-[1] "real"
-
-$dimensions
-[1] 1Executable location
-
-dir argument.
-
mod <- cmdstan_model(stan_file, dir = "path/to/directory/for/executable")Processing data
-
-
-
-Named list of R objects
-
-N, the
-number of data points, and y an integer array of
-observations.
-
mod$print()
-data {
- int<lower=0> N;
- array[N] int<lower=0, upper=1> y;
-}
-parameters {
- real<lower=0, upper=1> theta;
-}
-model {
- theta ~ beta(1, 1); // uniform prior on interval 0,1
- y ~ bernoulli(theta);
-}
-
# data block has 'N' and 'y'
-data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
-fit <- mod$sample(data = data_list)write_stan_json(). This happens internally, but it is also
-possible to call write_stan_json() directly.
-
data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
-json_file <- tempfile(fileext = ".json")
-write_stan_json(data_list, json_file)
-cat(readLines(json_file), sep = "\n")
-{
- "N": 10,
- "y": [0, 1, 0, 0, 0, 0, 0, 0, 0, 1]
-}JSON file
-
-write_stan_json():
-
fit <- mod$sample(data = json_file)R dump file
-
-rstan::stan_rdump():Writing CmdStan output to CSV
-
-Default temporary files
-
-
-
-
fit$output_files()
-[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-1-5c7cee.csv"
-[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-2-5c7cee.csv"
-[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-3-5c7cee.csv"
-[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-4-5c7cee.csv"fit object and force (or wait for) garbage
-collection.
-
files <- fit$output_files()
-file.exists(files)
-
-[1] TRUE TRUE TRUE TRUE
- used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
-Ncells 1161565 62.1 2340494 125 NA 1595651 85.3
-Vcells 2051945 15.7 8388608 64 32768 4957304 37.9
-
file.exists(files)
-[1] FALSE FALSE FALSE FALSENon-temporary files
-
-output_dir argument to
-mod$sample() or use fit$save_output_files()
-after fitting the model.
-
# see ?save_output_files for info on optional arguments
-fit$save_output_files(dir = "path/to/directory")
-
fit <- mod$sample(
- data = data_list,
- output_dir = "path/to/directory"
-)Reading CmdStan output into R
-
-Lazy CSV reading
-
-fit$draws(),
-fit$summary(), etc.). If we examine the structure of the
-fit object, notice how the Private slot
-draws_ is NULL, indicating that the CSV files
-haven’t yet been read into R.
-
str(fit)
-Classes 'CmdStanMCMC', 'CmdStanFit', 'R6' <CmdStanMCMC>
- Inherits from: <CmdStanFit>
- Public:
- clone: function (deep = FALSE)
- cmdstan_diagnose: function ()
- cmdstan_summary: function (flags = NULL)
- code: function ()
- config_files: function (include_failed = FALSE)
- constrain_variables: function (unconstrained_variables, transformed_parameters = TRUE,
- data_file: function ()
- diagnostic_summary: function (diagnostics = c("divergences", "treedepth", "ebfmi"),
- draws: function (variables = NULL, inc_warmup = FALSE, format = getOption("cmdstanr_draws_format",
- expose_functions: function (global = FALSE, verbose = FALSE)
- functions: environment
- grad_log_prob: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- hessian: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- init: function ()
- init_model_methods: function (seed = 1, verbose = FALSE, hessian = FALSE)
- initialize: function (runset)
- inv_metric: function (matrix = TRUE)
- latent_dynamics_files: function (include_failed = FALSE)
- log_prob: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- loo: function (variables = "log_lik", r_eff = TRUE, moment_match = FALSE,
- lp: function ()
- metadata: function ()
- metric_files: function (include_failed = FALSE)
- num_chains: function ()
- num_procs: function ()
- output: function (id = NULL)
- output_files: function (include_failed = FALSE)
- print: function (variables = NULL, ..., digits = 2, max_rows = getOption("cmdstanr_max_rows",
- profile_files: function (include_failed = FALSE)
- profiles: function ()
- return_codes: function ()
- runset: CmdStanRun, R6
- sampler_diagnostics: function (inc_warmup = FALSE, format = getOption("cmdstanr_draws_format",
- save_config_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_data_file: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_latent_dynamics_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_metric_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_object: function (file, ...)
- save_output_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_profile_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- summary: function (variables = NULL, ...)
- time: function ()
- unconstrain_draws: function (files = NULL, draws = NULL, format = getOption("cmdstanr_draws_format",
- unconstrain_variables: function (variables)
- variable_skeleton: function (transformed_parameters = TRUE, generated_quantities = TRUE)
- Private:
- draws_: NULL
- init_: NULL
- inv_metric_: list
- metadata_: list
- model_methods_env_: environment
- profiles_: NULL
- read_csv_: function (variables = NULL, sampler_diagnostics = NULL, format = getOption("cmdstanr_draws_format",
- return_codes_: 0 0 0 0
- sampler_diagnostics_: 2 1 1 2 2 1 2 2 1 2 1 1 1 2 2 1 2 1 1 1 1 2 2 2 1 2 1 1 ...
- warmup_draws_: NULL
- warmup_sampler_diagnostics_: NULL draws_
-slot in Private is no longer empty.
-
draws <- fit$draws() # force CSVs to be read into R
-str(fit)
-Classes 'CmdStanMCMC', 'CmdStanFit', 'R6' <CmdStanMCMC>
- Inherits from: <CmdStanFit>
- Public:
- clone: function (deep = FALSE)
- cmdstan_diagnose: function ()
- cmdstan_summary: function (flags = NULL)
- code: function ()
- config_files: function (include_failed = FALSE)
- constrain_variables: function (unconstrained_variables, transformed_parameters = TRUE,
- data_file: function ()
- diagnostic_summary: function (diagnostics = c("divergences", "treedepth", "ebfmi"),
- draws: function (variables = NULL, inc_warmup = FALSE, format = getOption("cmdstanr_draws_format",
- expose_functions: function (global = FALSE, verbose = FALSE)
- functions: environment
- grad_log_prob: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- hessian: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- init: function ()
- init_model_methods: function (seed = 1, verbose = FALSE, hessian = FALSE)
- initialize: function (runset)
- inv_metric: function (matrix = TRUE)
- latent_dynamics_files: function (include_failed = FALSE)
- log_prob: function (unconstrained_variables, jacobian = TRUE, jacobian_adjustment = NULL)
- loo: function (variables = "log_lik", r_eff = TRUE, moment_match = FALSE,
- lp: function ()
- metadata: function ()
- metric_files: function (include_failed = FALSE)
- num_chains: function ()
- num_procs: function ()
- output: function (id = NULL)
- output_files: function (include_failed = FALSE)
- print: function (variables = NULL, ..., digits = 2, max_rows = getOption("cmdstanr_max_rows",
- profile_files: function (include_failed = FALSE)
- profiles: function ()
- return_codes: function ()
- runset: CmdStanRun, R6
- sampler_diagnostics: function (inc_warmup = FALSE, format = getOption("cmdstanr_draws_format",
- save_config_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_data_file: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_latent_dynamics_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_metric_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_object: function (file, ...)
- save_output_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- save_profile_files: function (dir = ".", basename = NULL, timestamp = TRUE, random = TRUE)
- summary: function (variables = NULL, ...)
- time: function ()
- unconstrain_draws: function (files = NULL, draws = NULL, format = getOption("cmdstanr_draws_format",
- unconstrain_variables: function (variables)
- variable_skeleton: function (transformed_parameters = TRUE, generated_quantities = TRUE)
- Private:
- draws_: -7.45006 -7.1701 -7.11437 -6.74834 -7.04353 -6.79703 -7. ...
- init_: NULL
- inv_metric_: list
- metadata_: list
- model_methods_env_: environment
- profiles_: NULL
- read_csv_: function (variables = NULL, sampler_diagnostics = NULL, format = getOption("cmdstanr_draws_format",
- return_codes_: 0 0 0 0
- sampler_diagnostics_: 2 1 1 2 2 1 2 2 1 2 1 1 1 2 2 1 2 1 1 1 1 2 2 2 1 2 1 1 ...
- warmup_draws_: NULL
- warmup_sampler_diagnostics_: NULL variables argument to fit$draws()) then
-CmdStanR will only read in the requested variables (unless they have
-already been read in).read_cmdstan_csv()
-
-read_cmdstan_csv() function is used to
-read the CmdStan CSV files into R. This function is exposed to users, so
-you can also call it directly.
-
# see ?read_cmdstan_csv for info on optional arguments controlling
-# what information is read in
-csv_contents <- read_cmdstan_csv(fit$output_files())
-str(csv_contents)
-List of 8
- $ metadata :List of 42
- ..$ stan_version_major : num 2
- ..$ stan_version_minor : num 36
- ..$ stan_version_patch : num 0
- ..$ start_datetime : chr "2025-03-31 14:50:13 UTC"
- ..$ method : chr "sample"
- ..$ save_warmup : int 0
- ..$ thin : num 1
- ..$ gamma : num 0.05
- ..$ kappa : num 0.75
- ..$ t0 : num 10
- ..$ init_buffer : num 75
- ..$ term_buffer : num 50
- ..$ window : num 25
- ..$ save_metric : int 0
- ..$ algorithm : chr "hmc"
- ..$ engine : chr "nuts"
- ..$ metric : chr "diag_e"
- ..$ stepsize_jitter : num 0
- ..$ num_chains : num 1
- ..$ id : num [1:4] 1 2 3 4
- ..$ init : num [1:4] 2 2 2 2
- ..$ seed : num 31749990
- ..$ refresh : num 100
- ..$ sig_figs : num -1
- ..$ profile_file : chr "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-profile-202503310850-1-2d08be.csv"
- ..$ save_cmdstan_config : int 0
- ..$ stanc_version : chr "stanc3 v2.36.0"
- ..$ sampler_diagnostics : chr [1:6] "accept_stat__" "stepsize__" "treedepth__" "n_leapfrog__" ...
- ..$ variables : chr [1:2] "lp__" "theta"
- ..$ step_size_adaptation: num [1:4] 0.973 0.888 0.942 1.142
- ..$ model_name : chr "bernoulli_model"
- ..$ adapt_engaged : int 1
- ..$ adapt_delta : num 0.8
- ..$ max_treedepth : num 10
- ..$ step_size : num [1:4] 1 1 1 1
- ..$ iter_warmup : num 1000
- ..$ iter_sampling : num 1000
- ..$ threads_per_chain : num 1
- ..$ time :'data.frame': 4 obs. of 4 variables:
- .. ..$ chain_id: num [1:4] 1 2 3 4
- .. ..$ warmup : num [1:4] 0.004 0.004 0.004 0.005
- .. ..$ sampling: num [1:4] 0.014 0.014 0.012 0.013
- .. ..$ total : num [1:4] 0.018 0.018 0.016 0.018
- ..$ stan_variable_sizes :List of 2
- .. ..$ lp__ : num 1
- .. ..$ theta: num 1
- ..$ stan_variables : chr [1:2] "lp__" "theta"
- ..$ model_params : chr [1:2] "lp__" "theta"
- $ time :List of 2
- ..$ total : int NA
- ..$ chains:'data.frame': 4 obs. of 4 variables:
- .. ..$ chain_id: num [1:4] 1 2 3 4
- .. ..$ warmup : num [1:4] 0.004 0.004 0.004 0.005
- .. ..$ sampling: num [1:4] 0.014 0.014 0.012 0.013
- .. ..$ total : num [1:4] 0.018 0.018 0.016 0.018
- $ inv_metric :List of 4
- ..$ 1: num 0.489
- ..$ 2: num 0.45
- ..$ 3: num 0.577
- ..$ 4: num 0.541
- $ step_size :List of 4
- ..$ 1: num 0.973
- ..$ 2: num 0.888
- ..$ 3: num 0.942
- ..$ 4: num 1.14
- $ warmup_draws : NULL
- $ post_warmup_draws : 'draws_array' num [1:1000, 1:4, 1:2] -7.45 -7.17 -7.11 -6.75 -7.04 ...
- ..- attr(*, "dimnames")=List of 3
- .. ..$ iteration: chr [1:1000] "1" "2" "3" "4" ...
- .. ..$ chain : chr [1:4] "1" "2" "3" "4"
- .. ..$ variable : chr [1:2] "lp__" "theta"
- $ warmup_sampler_diagnostics : NULL
- $ post_warmup_sampler_diagnostics: 'draws_array' num [1:1000, 1:4, 1:6] 0.9 1 1 0.911 0.801 ...
- ..- attr(*, "dimnames")=List of 3
- .. ..$ iteration: chr [1:1000] "1" "2" "3" "4" ...
- .. ..$ chain : chr [1:4] "1" "2" "3" "4"
- .. ..$ variable : chr [1:6] "accept_stat__" "stepsize__" "treedepth__" "n_leapfrog__" ...as_cmdstan_fit()
-
-as_cmdstan_fit().
-
fit2 <- as_cmdstan_fit(fit$output_files())fit object, but this function can be used to create fitted
-model objects (CmdStanMCMC, CmdStanMLE, etc.)
-from any CmdStan CSV files.Saving and accessing advanced algorithm info (latent dynamics)
-
-save_latent_dynamics is set to TRUE when
-running the $sample() method then additional CSV files are
-created (one per chain) that provide access to quantities used under the
-hood by Stan’s implementation of dynamic Hamiltonian Monte Carlo.
-
fit <- mod$sample(data = data_list, save_latent_dynamics = TRUE)
-
fit$latent_dynamics_files()
-[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-diagnostic-202503310850-1-060b43.csv"
-[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-diagnostic-202503310850-2-060b43.csv"
-[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-diagnostic-202503310850-3-060b43.csv"
-[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-diagnostic-202503310850-4-060b43.csv"
-
# read one of the files in
-x <- utils::read.csv(fit$latent_dynamics_files()[1], comment.char = "#")
-head(x)
- lp__ accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__
-1 -7.00899 0.935119 1.06923 1 3 0
-2 -7.29900 0.886772 1.06923 1 1 0
-3 -6.82892 1.000000 1.06923 1 1 0
-4 -6.82025 1.000000 1.06923 1 1 0
-5 -6.75311 0.993504 1.06923 2 3 0
-6 -6.89721 0.960772 1.06923 1 3 0
- energy__ theta p_theta g_theta
-1 7.05933 -0.633893 0.440038 1.159540
-2 7.38011 -0.432323 -0.558539 1.722860
-3 7.12220 -0.836018 -1.062060 0.628489
-4 6.85073 -0.850225 0.342380 0.592629
-5 6.85434 -1.166260 0.623978 -0.149621
-6 6.92007 -0.744453 0.296526 0.864374lp__ is also provided via
-fit$draws(), and the columns accept_stat__,
-stepsize__, treedepth__,
-n_leapfrog__, divergent__, and
-energy__ are also provided by
-fit$sampler_diagnostics(), but there are several columns
-unique to the latent dynamics file.
- theta p_theta g_theta
-1 -0.633893 0.440038 1.159540
-2 -0.432323 -0.558539 1.722860
-3 -0.836018 -1.062060 0.628489
-4 -0.850225 0.342380 0.592629
-5 -1.166260 0.623978 -0.149621
-6 -0.744453 0.296526 0.864374theta and the three
-columns above correspond to theta in the
-unconstrained space (theta on the constrained
-space is accessed via fit$draws()), the auxiliary momentum
-p_theta, and the gradient g_theta. In general,
-each of these three columns will exist for every parameter in
-the model.Developing using CmdStanR
-
-Pre-compiled Stan models in R packages
-
-configure and configure.win which run
-cmdstan_model() with compile = TRUE and save
-the compiled executables somewhere inside the inst/ folder
-of the package source. The instantiate
-package helps developers configure packages this way, and it documents
-other topics such as submitting to CRAN and administering CmdStan. Kevin
-Ushey’s configure
-package helps create and manage package configuration files in
-general.Troubleshooting and debugging
-
-"cmdstanr_verbose" to TRUE.
-
options("cmdstanr_verbose"=TRUE)
-
-mod <- cmdstan_model(stan_file, force_recompile = TRUE)
-Running make \
- /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19 \
- "STANCFLAGS += --name='bernoulli_model'"
-
---- Translating Stan model to C++ code ---
-bin/stanc --name='bernoulli_model' --o=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.hpp /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.stan
-
---- Compiling C++ code ---
-clang++ -O3 -march=native -mtune=native -Wno-deprecated-declarations -std=c++17 -Wno-unknown-warning-option -Wno-tautological-compare -Wno-sign-compare -D_REENTRANT -Wno-ignored-attributes -I stan/lib/stan_math/lib/tbb_2020.3/include -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.4.0 -I stan/lib/stan_math/lib/boost_1.84.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials -DBOOST_DISABLE_ASSERTS -c -include-pch stan/src/stan/model/model_header.hpp.gch/model_header_16_0.hpp.gch -x c++ -o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.hpp
-
---- Linking model ---
-clang++ -O3 -march=native -mtune=native -Wno-deprecated-declarations -std=c++17 -Wno-unknown-warning-option -Wno-tautological-compare -Wno-sign-compare -D_REENTRANT -Wno-ignored-attributes -I stan/lib/stan_math/lib/tbb_2020.3/include -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.4.0 -I stan/lib/stan_math/lib/boost_1.84.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials -DBOOST_DISABLE_ASSERTS -Wl,-L,"/Users/jgabry/.cmdstan/cmdstan-2.36.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/Users/jgabry/.cmdstan/cmdstan-2.36.0/stan/lib/stan_math/lib/tbb" /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.o src/cmdstan/main.o -ltbb stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a stan/lib/stan_math/lib/tbb/libtbb.dylib stan/lib/stan_math/lib/tbb/libtbbmalloc.dylib stan/lib/stan_math/lib/tbb/libtbbmalloc_proxy.dylib -o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19
-rm /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/model-d1b65860dc19.hpp
-
fit <- mod$sample(
- data = data_list,
- chains = 1,
- iter_warmup = 100,
- iter_sampling = 100
-)
-Running MCMC with 1 chain...
-
-Running ./bernoulli 'id=1' random 'seed=1376020223' data \
- 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/standata-d1b61e77387.json' \
- output \
- 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-1-1ead10.csv' \
- 'profile_file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-profile-202503310850-1-3e64eb.csv' \
- 'method=sample' 'num_samples=100' 'num_warmup=100' 'save_warmup=0' \
- 'algorithm=hmc' 'engine=nuts' adapt 'engaged=1'
-Chain 1 method = sample (Default)
-Chain 1 sample
-Chain 1 num_samples = 100
-Chain 1 num_warmup = 100
-Chain 1 save_warmup = false (Default)
-Chain 1 thin = 1 (Default)
-Chain 1 adapt
-Chain 1 engaged = true (Default)
-Chain 1 gamma = 0.05 (Default)
-Chain 1 delta = 0.8 (Default)
-Chain 1 kappa = 0.75 (Default)
-Chain 1 t0 = 10 (Default)
-Chain 1 init_buffer = 75 (Default)
-Chain 1 term_buffer = 50 (Default)
-Chain 1 window = 25 (Default)
-Chain 1 save_metric = false (Default)
-Chain 1 algorithm = hmc (Default)
-Chain 1 hmc
-Chain 1 engine = nuts (Default)
-Chain 1 nuts
-Chain 1 max_depth = 10 (Default)
-Chain 1 metric = diag_e (Default)
-Chain 1 metric_file = (Default)
-Chain 1 stepsize = 1 (Default)
-Chain 1 stepsize_jitter = 0 (Default)
-Chain 1 num_chains = 1 (Default)
-Chain 1 id = 1 (Default)
-Chain 1 data
-Chain 1 file = /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/standata-d1b61e77387.json
-Chain 1 init = 2 (Default)
-Chain 1 random
-Chain 1 seed = 1376020223
-Chain 1 output
-Chain 1 file = /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-202503310850-1-1ead10.csv
-Chain 1 diagnostic_file = (Default)
-Chain 1 refresh = 100 (Default)
-Chain 1 sig_figs = -1 (Default)
-Chain 1 profile_file = /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpI5U8HG/bernoulli-profile-202503310850-1-3e64eb.csv
-Chain 1 save_cmdstan_config = false (Default)
-Chain 1 num_threads = 1 (Default)
-Chain 1 Gradient evaluation took 1e-05 seconds
-Chain 1 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
-Chain 1 Adjust your expectations accordingly!
-Chain 1 WARNING: There aren't enough warmup iterations to fit the
-Chain 1 three stages of adaptation as currently configured.
-Chain 1 Reducing each adaptation stage to 15%/75%/10% of
-Chain 1 the given number of warmup iterations:
-Chain 1 init_buffer = 15
-Chain 1 adapt_window = 75
-Chain 1 term_buffer = 10
-Chain 1 Iteration: 1 / 200 [ 0%] (Warmup)
-Chain 1 Iteration: 100 / 200 [ 50%] (Warmup)
-Chain 1 Iteration: 101 / 200 [ 50%] (Sampling)
-Chain 1 Iteration: 200 / 200 [100%] (Sampling)
-Chain 1 Elapsed Time: 0 seconds (Warm-up)
-Chain 1 0.001 seconds (Sampling)
-Chain 1 0.001 seconds (Total)
-Chain 1 finished in 0.0 seconds.Getting started with CmdStanR
- Jonah Gabry,
-Rok Češnovar, and Andrew Johnson
-
-
- Source: vignettes/cmdstanr.Rmd
- cmdstanr.RmdIntroduction
-
-
-
# we recommend running this is a fresh R session or restarting your current session
-install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))Installing CmdStan
-
-check_cmdstan_toolchain() function:
-The C++ toolchain required for CmdStan is setup properly!install_cmdstan()
-function:
-
install_cmdstan(cores = 2)
-
-"CMDSTAN" exists at load
-time then its value will be automatically set as the default path to
-CmdStan for the R session. This is useful if your CmdStan installation
-is not located in the default directory that would have been used by
-install_cmdstan() (see #2).".cmdstan/cmdstan-[version]", for example
-".cmdstan/cmdstan-2.23.0", exists in the user’s home
-directory (Sys.getenv("HOME"), not the current
-working directory) then the path to the CmdStan with the largest version
-number will be set as the path to CmdStan for the R session. This is the
-same as the default directory that install_cmdstan() uses
-to install the latest version of CmdStan, so if that’s how you installed
-CmdStan you shouldn’t need to manually set the path to CmdStan when
-loading CmdStanR.set_cmdstan_path() function:
-
set_cmdstan_path(PATH_TO_CMDSTAN)cmdstan_path() and
-cmdstan_version():
-
cmdstan_path()
-
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.36.0"
-[1] "2.36.0"Compiling a model
-
-cmdstan_model() function creates a new CmdStanModel
-object from a file containing a Stan program. Under the hood, CmdStan is
-called to translate a Stan program to C++ and create a compiled
-executable. Here we’ll use the example Stan program that comes with the
-CmdStan installation:
-
file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
-mod <- cmdstan_model(file)mod is an R6 reference object of class CmdStanModel
-and behaves similarly to R’s reference class objects and those in object
-oriented programming languages. Methods are accessed using the
-$ operator. This design choice allows for CmdStanR and CmdStanPy to provide a
-similar user experience and share many implementation details.$print()
-method:
-
mod$print()
-data {
- int<lower=0> N;
- array[N] int<lower=0, upper=1> y;
-}
-parameters {
- real<lower=0, upper=1> theta;
-}
-model {
- theta ~ beta(1, 1); // uniform prior on interval 0,1
- y ~ bernoulli(theta);
-}$exe_file() method:
-
mod$exe_file()
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.36.0/examples/bernoulli/bernoulli"Running MCMC
-
-$sample()
-method for CmdStanModel
-objects runs Stan’s default MCMC algorithm. The data
-argument accepts a named list of R objects (like for RStan) or a path to
-a data file compatible with CmdStan (JSON or R dump).
-
# names correspond to the data block in the Stan program
-data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
-
-fit <- mod$sample(
- data = data_list,
- seed = 123,
- chains = 4,
- parallel_chains = 4,
- refresh = 500 # print update every 500 iters
-)
-Running MCMC with 4 parallel chains...
-
-Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
-Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
-Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
-Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
-Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
-Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
-Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
-Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
-Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
-Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
-Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
-Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
-Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
-Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
-Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
-Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
-Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
-Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
-Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
-Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
-Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
-Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
-Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
-Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
-Chain 1 finished in 0.0 seconds.
-Chain 2 finished in 0.0 seconds.
-Chain 3 finished in 0.0 seconds.
-Chain 4 finished in 0.0 seconds.
-
-All 4 chains finished successfully.
-Mean chain execution time: 0.0 seconds.
-Total execution time: 0.4 seconds.$sample() method. For details follow this link to its
-separate documentation page:
-
-$sample()$sample() method creates R6 CmdStanMCMC objects,
-which have many associated methods. Below we will demonstrate some of
-the most important methods. For a full list, follow this link to the
-CmdStanMCMC documentation:
-
-CmdStanMCMCPosterior summary statistics
-
-Summaries from the posterior package
-
-$summary()
-method calls summarise_draws() from the
-posterior package. The first argument specifies the
-variables to summarize and any arguments after that are passed on to
-posterior::summarise_draws() to specify which summaries to
-compute, whether to use multiple cores, etc.
-
fit$summary()
-fit$summary(variables = c("theta", "lp__"), "mean", "sd")
-
-# use a formula to summarize arbitrary functions, e.g. Pr(theta <= 0.5)
-fit$summary("theta", pr_lt_half = ~ mean(. <= 0.5))
-
-# summarise all variables with default and additional summary measures
-fit$summary(
- variables = NULL,
- posterior::default_summary_measures(),
- extra_quantiles = ~posterior::quantile2(., probs = c(.0275, .975))
-)
- variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
-1 lp__ -7.30 -6.99 0.79 0.33 -8.927 -6.75 1 1769 1938
-2 theta 0.25 0.24 0.12 0.12 0.078 0.48 1 1228 1521
- variable mean sd
-1 theta 0.25 0.12
-2 lp__ -7.30 0.79
- variable pr_lt_half
-1 theta 0.96
- variable mean median sd mad q5 q95 q2.75 q97.5
-1 lp__ -7.30 -6.99 0.79 0.33 -8.927 -6.75 -9.429 -6.75
-2 theta 0.25 0.24 0.12 0.12 0.078 0.48 0.062 0.54Posterior draws
-
-Extracting draws
-
-$draws()
-method can be used to extract the posterior draws in formats provided by
-the posterior
-package. Here we demonstrate only the draws_array and
-draws_df formats, but the posterior
-package supports other useful formats as well.
-
# default is a 3-D draws_array object from the posterior package
-# iterations x chains x variables
-draws_arr <- fit$draws() # or format="array"
-str(draws_arr)
- 'draws_array' num [1:1000, 1:4, 1:2] -7.01 -7.89 -7.41 -6.75 -6.91 ...
- - attr(*, "dimnames")=List of 3
- ..$ iteration: chr [1:1000] "1" "2" "3" "4" ...
- ..$ chain : chr [1:4] "1" "2" "3" "4"
- ..$ variable : chr [1:2] "lp__" "theta"
-
# draws x variables data frame
-draws_df <- fit$draws(format = "df")
-str(draws_df)
-draws_df [4,000 × 5] (S3: draws_df/draws/tbl_df/tbl/data.frame)
- $ lp__ : num [1:4000] -7.01 -7.89 -7.41 -6.75 -6.91 ...
- $ theta : num [1:4000] 0.168 0.461 0.409 0.249 0.185 ...
- $ .chain : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
- $ .iteration: int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
- $ .draw : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
-
print(draws_df)
-# A draws_df: 1000 iterations, 4 chains, and 2 variables
- lp__ theta
-1 -7.0 0.17
-2 -7.9 0.46
-3 -7.4 0.41
-4 -6.7 0.25
-5 -6.9 0.18
-6 -6.9 0.33
-7 -7.2 0.15
-8 -6.8 0.29
-9 -6.8 0.24
-10 -6.8 0.24
-# ... with 3990 more draws
-# ... hidden reserved variables {'.chain', '.iteration', '.draw'}posterior::as_draws_*() functions.
-
# this should be identical to draws_df created via draws(format = "df")
-draws_df_2 <- as_draws_df(draws_arr)
-identical(draws_df, draws_df_2)
-[1] TRUE$draws() method, but in most cases the speed
-difference will be minor.rstan::extract().Sampler diagnostics
-
-Extracting diagnostic values for each iteration and chain
-
-$sampler_diagnostics()
-method extracts the values of the sampler parameters
-(treedepth__, divergent__, etc.) in formats
-supported by the posterior package. The default is as a
-3-D array (iteration x chain x variable).
-
# this is a draws_array object from the posterior package
-str(fit$sampler_diagnostics())
- 'draws_array' num [1:1000, 1:4, 1:6] 2 1 2 2 2 1 2 1 2 1 ...
- - attr(*, "dimnames")=List of 3
- ..$ iteration: chr [1:1000] "1" "2" "3" "4" ...
- ..$ chain : chr [1:4] "1" "2" "3" "4"
- ..$ variable : chr [1:6] "treedepth__" "divergent__" "energy__" "accept_stat__" ...
-
# this is a draws_df object from the posterior package
-str(fit$sampler_diagnostics(format = "df"))
-draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
- $ treedepth__ : num [1:4000] 2 1 2 2 2 1 2 1 2 1 ...
- $ divergent__ : num [1:4000] 0 0 0 0 0 0 0 0 0 0 ...
- $ energy__ : num [1:4000] 8.95 8.77 7.87 7.64 6.93 ...
- $ accept_stat__: num [1:4000] 0.688 0.811 1 0.966 0.976 ...
- $ stepsize__ : num [1:4000] 0.905 0.905 0.905 0.905 0.905 ...
- $ n_leapfrog__ : num [1:4000] 3 3 3 3 3 3 3 3 3 3 ...
- $ .chain : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
- $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
- $ .draw : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...Sampler diagnostic warnings and summaries
-
-$diagnostic_summary() method will display any
-sampler diagnostic warnings and return a summary of diagnostics for each
-chain.
-
fit$diagnostic_summary()
-$num_divergent
-[1] 0 0 0 0
-
-$num_max_treedepth
-[1] 0 0 0 0
-
-$ebfmi
-[1] 1.11 0.76 1.19 1.08
-
fit_with_warning <- cmdstanr_example("schools")
-Warning: 143 of 4000 (4.0%) transitions ended with a divergence.
-See https://mc-stan.org/misc/warnings for details.
-Warning: 1 of 4 chains had an E-BFMI less than 0.3.
-See https://mc-stan.org/misc/warnings for details.fit$diagnostic_summary().
-
diagnostics <- fit_with_warning$diagnostic_summary()
-Warning: 143 of 4000 (4.0%) transitions ended with a divergence.
-See https://mc-stan.org/misc/warnings for details.
-Warning: 1 of 4 chains had an E-BFMI less than 0.3.
-See https://mc-stan.org/misc/warnings for details.
-
print(diagnostics)
-$num_divergent
-[1] 1 37 75 30
-
-$num_max_treedepth
-[1] 0 0 0 0
-
-$ebfmi
-[1] 0.17 0.35 0.34 0.42
-
# number of divergences reported in warning is the sum of the per chain values
-sum(diagnostics$num_divergent)
-[1] 143Running optimization and variational inference
-
-$optimize(), $laplace(),
-$variational(), and $pathfinder() methods,
-which are called in a similar way to the $sample() method
-demonstrated above.Optimization
-
-$optimize().
-
fit_mle <- mod$optimize(data = data_list, seed = 123)
-Initial log joint probability = -16.144
- Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
- 6 -5.00402 0.000246518 8.73164e-07 1 1 9
-Optimization terminated normally:
- Convergence detected: relative gradient magnitude is below tolerance
-Finished in 0.2 seconds.
-
fit_mle$print() # includes lp__ (log prob calculated by Stan program)
- variable estimate
- lp__ -5.00
- theta 0.20
-
fit_mle$mle("theta")
-theta
- 0.2 theta.
jacobian=TRUE. See the Maximum
-Likelihood Estimation section of the CmdStan User’s Guide for more
-details.
-
fit_map <- mod$optimize(
- data = data_list,
- jacobian = TRUE,
- seed = 123
-)
-Initial log joint probability = -18.2733
- Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
- 5 -6.74802 0.000708195 1.43227e-05 1 1 8
-Optimization terminated normally:
- Convergence detected: relative gradient magnitude is below tolerance
-Finished in 0.1 seconds.Laplace Approximation
-
-$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 MAP
-estimate, the samples provide an estimate of the mean and standard
-deviation of the posterior distribution. If the mode is the MLE, the
-sample provides an estimate of the standard error of the likelihood.
-Whether the mode is the MAP or MLE depends on the value of the
-jacobian argument when running optimization. See the Laplace
-Sampling chapter of the CmdStan User’s Guide for more details.fit_map object from above as the
-mode argument. If mode is omitted then
-optimization will be run internally before taking draws from the normal
-approximation.
-
fit_laplace <- mod$laplace(
- mode = fit_map,
- draws = 4000,
- data = data_list,
- seed = 123,
- refresh = 1000
- )
-Calculating Hessian
-Calculating inverse of Cholesky factor
-Generating draws
-iteration: 0
-iteration: 1000
-iteration: 2000
-iteration: 3000
-Finished in 0.1 seconds.
-
fit_laplace$print("theta")
- variable mean median sd mad q5 q95
- theta 0.27 0.25 0.12 0.12 0.10 0.51
-
mcmc_hist(fit_laplace$draws("theta"), binwidth = 0.025)
Variational (ADVI)
-
-$variational()
-method. For details on the ADVI algorithm see the CmdStan
-User’s Guide.
-
fit_vb <- mod$variational(
- data = data_list,
- seed = 123,
- draws = 4000
-)
-------------------------------------------------------------
-EXPERIMENTAL ALGORITHM:
- This procedure has not been thoroughly tested and may be unstable
- or buggy. The interface is subject to change.
-------------------------------------------------------------
-Gradient evaluation took 9e-06 seconds
-1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
-Adjust your expectations accordingly!
-Begin eta adaptation.
-Iteration: 1 / 250 [ 0%] (Adaptation)
-Iteration: 50 / 250 [ 20%] (Adaptation)
-Iteration: 100 / 250 [ 40%] (Adaptation)
-Iteration: 150 / 250 [ 60%] (Adaptation)
-Iteration: 200 / 250 [ 80%] (Adaptation)
-Success! Found best value [eta = 1] earlier than expected.
-Begin stochastic gradient ascent.
- iter ELBO delta_ELBO_mean delta_ELBO_med notes
- 100 -6.164 1.000 1.000
- 200 -6.225 0.505 1.000
- 300 -6.186 0.339 0.010 MEDIAN ELBO CONVERGED
-Drawing a sample of size 4000 from the approximate posterior...
-COMPLETED.
-Finished in 0.1 seconds.
-
fit_vb$print("theta")
- variable mean median sd mad q5 q95
- theta 0.26 0.24 0.11 0.11 0.11 0.46
-
mcmc_hist(fit_vb$draws("theta"), binwidth = 0.025)
Variational (Pathfinder)
-
-$pathfinder()
-method.
-
fit_pf <- mod$pathfinder(
- data = data_list,
- seed = 123,
- draws = 4000
-)
-Path [1] :Initial log joint density = -18.273334
-Path [1] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes
- 5 -6.748e+00 7.082e-04 1.432e-05 1.000e+00 1.000e+00 126 -6.145e+00 -6.145e+00
-Path [1] :Best Iter: [5] ELBO (-6.145070) evaluations: (126)
-Path [2] :Initial log joint density = -19.192715
-Path [2] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes
- 5 -6.748e+00 2.015e-04 2.228e-06 1.000e+00 1.000e+00 126 -6.223e+00 -6.223e+00
-Path [2] :Best Iter: [2] ELBO (-6.170358) evaluations: (126)
-Path [3] :Initial log joint density = -6.774820
-Path [3] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes
- 4 -6.748e+00 1.137e-04 2.596e-07 1.000e+00 1.000e+00 101 -6.178e+00 -6.178e+00
-Path [3] :Best Iter: [4] ELBO (-6.177909) evaluations: (101)
-Path [4] :Initial log joint density = -7.949193
-Path [4] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes
- 5 -6.748e+00 2.145e-04 1.301e-06 1.000e+00 1.000e+00 126 -6.197e+00 -6.197e+00
-Path [4] :Best Iter: [5] ELBO (-6.197118) evaluations: (126)
-Total log probability function evaluations:4379
-Finished in 0.1 seconds.
-
fit_pf$print("theta")
- variable mean median sd mad q5 q95
- theta 0.25 0.24 0.12 0.12 0.08 0.47
-
mcmc_hist(fit_pf$draws("theta"), binwidth = 0.025) +
- ggplot2::labs(subtitle = "Approximate posterior from pathfinder") +
- ggplot2::xlim(0, 1)
-
mcmc_hist(fit_vb$draws("theta"), binwidth = 0.025) +
- ggplot2::labs(subtitle = "Approximate posterior from variational") +
- ggplot2::xlim(0, 1)
-
mcmc_hist(fit_laplace$draws("theta"), binwidth = 0.025) +
- ggplot2::labs(subtitle = "Approximate posterior from Laplace") +
- ggplot2::xlim(0, 1)
-
mcmc_hist(fit$draws("theta"), binwidth = 0.025) +
- ggplot2::labs(subtitle = "Posterior from MCMC") +
- ggplot2::xlim(0, 1)
$optimize(),
-$laplace(), $variational(), and
-pathfinder() methods, follow these links to their
-documentation pages.Saving fitted model objects
-
-$save_object()
-method provided by CmdStanR is the most convenient way to save a fitted
-model object to disk and ensure that all of the contents are available
-when reading the object back into R.
-
fit$save_object(file = "fit.RDS")
-
-# can be read back in using readRDS
-fit2 <- readRDS("fit.RDS")$save_object()
-could take a long time. $save_object()
-reads the CmdStan results files into memory, stores them in the model
-object, and saves the object with saveRDS(). To speed up
-the process, you can emulate $save_object()
-and replace saveRDS with the much faster
-qsave() function from the qs package.
-
# Load CmdStan output files into the fitted model object.
-fit$draws() # Load posterior draws into the object.
-try(fit$sampler_diagnostics(), silent = TRUE) # Load sampler diagnostics.
-try(fit$init(), silent = TRUE) # Load user-defined initial values.
-try(fit$profiles(), silent = TRUE) # Load profiling samples.
-
-# Save the object to a file.
-qs::qsave(x = fit, file = "fit.qs")
-
-# Read the object.
-fit2 <- qs::qread("fit.qs")
-
# Load posterior draws into the fitted model object and omit other output.
-fit$draws()
-
-# Save the object to a file.
-qs::qsave(x = fit, file = "fit.qs")
-
-# Read the object.
-fit2 <- qs::qread("fit.qs")Comparison with RStan
-
-Different ways of interfacing with Stan’s C++
-
-Advantages of RStan
-
-
-
-Advantages of CmdStanR
-
-
-
-cmdstanr::install_cmdstan().Additional resources
-
-
