From 92ce480101ee9ca0e522ce6eef4f341a76950094 Mon Sep 17 00:00:00 2001
From: jgabry Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson. Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson, Steve Bronder.Page not found (404)
diff --git a/docs/CONTRIBUTING.html b/docs/CONTRIBUTING.html
index ce8b14b16..af3c28f4c 100644
--- a/docs/CONTRIBUTING.html
+++ b/docs/CONTRIBUTING.html
@@ -1,5 +1,5 @@
-Code of Conduct
License
diff --git a/docs/LICENSE.html b/docs/LICENSE.html
index 9cfe0b946..8cb6242d8 100644
--- a/docs/LICENSE.html
+++ b/docs/LICENSE.html
@@ -1,5 +1,5 @@
-BSD 3-Clause License
diff --git a/docs/articles/articles-online-only/opencl.html b/docs/articles/articles-online-only/opencl.html
index af47c06b3..df160a0ad 100644
--- a/docs/articles/articles-online-only/opencl.html
+++ b/docs/articles/articles-online-only/opencl.html
@@ -12,7 +12,7 @@
-
+
@@ -40,7 +40,7 @@
@@ -226,8 +226,8 @@ Compiling a model with OpenCLhere for a
-list of functions with OpenCL support.
To build the model with OpenCL support, add
cpp_options = list(stan_opencl = TRUE)
at the compilation
step.
Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson.
+Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson, Steve Bronder.
data {
int<lower=0> N;
- array[N] int<lower=0,upper=1> y;
+ array[N] int<lower=0, upper=1> y;
}
parameters {
- real<lower=0,upper=1> theta;
+ real<lower=0, upper=1> theta;
}
model {
- theta ~ beta(1,1); // uniform prior on interval 0,1
+ theta ~ beta(1, 1); // uniform prior on interval 0,1
y ~ bernoulli(theta);
}
mod$stan_file()
[1] "/Users/jgabry/.cmdstan/cmdstan-2.33.1/examples/bernoulli/bernoulli.stan"
+[1] "/Users/jgabry/.cmdstan/cmdstan-2.35.0/examples/bernoulli/bernoulli.stan"
mod$exe_file()
[1] "/Users/jgabry/.cmdstan/cmdstan-2.33.1/examples/bernoulli/bernoulli"
+[1] "/Users/jgabry/.cmdstan/cmdstan-2.35.0/examples/bernoulli/bernoulli"
Subsequently, if you create a CmdStanModel
object from
the same Stan file then compilation will be skipped (assuming the file
hasn’t changed).
mod$compile()
mod$exe_file()
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.33.1/examples/bernoulli/bernoulli"
+[1] "/Users/jgabry/.cmdstan/cmdstan-2.35.0/examples/bernoulli/bernoulli"
mod_pedantic <- cmdstan_model(stan_file_pedantic, pedantic = TRUE) -Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-74281eb75594.stan', line 11, column 14: A +Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a1418059d6f.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 @@ -288,7 +288,7 @@
argument to thePedantic checkpedantic
$check_syntax()
method.+mod_pedantic$check_syntax(pedantic = TRUE) -Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model_febb1e69c7387a0e64cf13583e078104.stan', line 11, column 14: A +Warning in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/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 @@ -300,18 +300,18 @@
Pedantic check
+[1] TRUEfile.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/RtmpXOnqh8/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 correct
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/Rtmp2hQnhs/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 correct
+stan_file_variables <- write_stan_file(" data { int<lower=1> J; @@ -351,47 +351,47 @@
Stan model variables
+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)
Each variable is represented as a list containing the type information (currently limited to
-real
orint
) and the number of dimensions.+variables$data$J
-$type [1] "int" $dimensions [1] 0
+variables$data$sigma
-$type [1] "real" $dimensions [1] 1
@@ -427,20 +427,20 @@+variables$parameters$tau
-$type [1] "real" $dimensions [1] 0
+variables$transformed_parameters$theta
$type [1] "real" @@ -405,7 +405,7 @@
Executable location
+mod <- cmdstan_model(stan_file, dir = "path/to/directory/for/executable")
Named list of R objectsN, 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; + array[N] int<lower=0, upper=1> y; } parameters { - real<lower=0,upper=1> theta; + real<lower=0, upper=1> theta; } model { - theta ~ beta(1,1); // uniform prior on interval 0,1 + theta ~ beta(1, 1); // uniform prior on interval 0,1 y ~ bernoulli(theta); }
+@@ -448,7 +448,7 @@# 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)
Named list of R objectswrite_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) @@ -465,7 +465,7 @@
: -JSON filewrite_stan_json()
+fit <- mod$sample(data = json_file)
@@ -475,7 +475,7 @@R dump filerstan::stan_rdump(): -
+rdump_file <- tempfile(fileext = ".data.R") rstan::stan_rdump(names(data_list), file = rdump_file, envir = list2env(data_list)) cat(readLines(rdump_file), sep = "\n") @@ -488,31 +488,31 @@
Writing CmdStan output to CSV
Default temporary files
-+When fitting a model, the default behavior is to write the output from CmdStan to CSV files in a temporary directory.
-+-fit$output_files()
+[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-202312131009-1-37e810.csv" -[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-202312131009-2-37e810.csv" -[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-202312131009-3-37e810.csv" -[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-202312131009-4-37e810.csv"
[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-202407021546-1-2808db.csv" +[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-202407021546-2-2808db.csv" +[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-202407021546-3-2808db.csv" +[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-202407021546-4-2808db.csv"
These files will be lost if you end your R session or if you remove the
-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 1260600 67.4 2436292 130.2 NA 2436292 130.2 -Vcells 2212085 16.9 8388608 64.0 32768 3459642 26.4
@@ -523,10 +523,10 @@++used (Mb) gc trigger (Mb) limit (Mb) max used (Mb) +Ncells 1155803 61.8 2350794 125.6 NA 1542121 82.4 +Vcells 2048173 15.7 8388608 64.0 32768 2851724 21.8
file.exists(files)
[1] FALSE FALSE FALSE FALSE
Non-temporary files
+-# see ?save_output_files for info on optional arguments fit$save_output_files(dir = "path/to/directory")
++[1] 1.11 0.80 1.05 0.95fit <- mod$sample( data = data_list, output_dir = "path/to/directory" @@ -546,7 +546,7 @@
object, notice how theLazy CSV readingfit
Private
slotdraws_
isNULL
, indicating that the CSV files haven’t yet been read into R. -+str(fit)
Classes 'CmdStanMCMC', 'CmdStanFit', 'R6' <CmdStanMCMC> Inherits from: <CmdStanFit> @@ -555,6 +555,7 @@
Lazy CSV readingLazy CSV readingLazy CSV readingLazy CSV readingLazy CSV reading
After we call a method that requires the draws then if we reexamine the structure of the object we will see that the
-draws_
slot inPrivate
is no longer empty.+draws <- fit$draws() # force CSVs to be read into R str(fit)
Classes 'CmdStanMCMC', 'CmdStanFit', 'R6' <CmdStanMCMC> @@ -617,6 +621,7 @@
Lazy CSV readingLazy CSV readingLazy CSV readingLazy CSV readingLazy CSV reading
For models with many parameters, transformed parameters, or generated @@ -678,19 +686,19 @@
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 40 + $ metadata :List of 42 ..$ stan_version_major : num 2 - ..$ stan_version_minor : num 33 + ..$ stan_version_minor : num 35 ..$ stan_version_patch : num 0 - ..$ start_datetime : chr "2023-12-13 17:09:13 UTC" + ..$ start_datetime : chr "2024-07-02 21:46:33 UTC" ..$ method : chr "sample" - ..$ save_warmup : num 0 + ..$ save_warmup : int 0 ..$ thin : num 1 ..$ gamma : num 0.05 ..$ kappa : num 0.75 @@ -698,6 +706,7 @@
read_cmdstan_csv()read_cmdstan_csv()read_cmdstan_csv()read_cmdstan_csv()as_cmdstan_fit()
If you need to manually create fitted model objects from CmdStan CSV files use
-as_cmdstan_fit()
.+fit2 <- as_cmdstan_fit(fit$output_files())
This is pointless in our case since we have the original
fit
object, but this function can be used to create fitted @@ -783,32 +793,32 @@Saving and
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.
--fit <- mod$sample(data = data_list, save_latent_dynamics = TRUE)
++fit <- mod$sample(data = data_list, save_latent_dynamics = TRUE)
-fit$latent_dynamics_files()
-[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-diagnostic-202312131009-1-425883.csv" -[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-diagnostic-202312131009-2-425883.csv" -[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-diagnostic-202312131009-3-425883.csv" -[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-diagnostic-202312131009-4-425883.csv"
-++[1] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-diagnostic-202407021546-1-54a823.csv" +[2] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-diagnostic-202407021546-2-54a823.csv" +[3] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-diagnostic-202407021546-3-54a823.csv" +[4] "/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-diagnostic-202407021546-4-54a823.csv"
# read one of the files in x <- utils::read.csv(fit$latent_dynamics_files()[1], comment.char = "#") head(x)
+1 -6.77025 0.987691 0.980824 2 3 0 +2 -7.20991 0.932461 0.980824 2 3 0 +3 -7.66830 0.933065 0.980824 1 1 0 +4 -7.37078 1.000000 0.980824 1 1 0 +5 -7.59961 0.948892 0.980824 2 3 0 +6 -7.95801 0.940559 0.980824 2 3 0 + energy__ theta p_theta g_theta +1 7.08007 -0.959636 -1.160790 0.323414 +2 7.21726 -1.779880 0.178790 -1.268180 +3 7.67261 -2.089350 -0.136971 -1.678370 +4 7.66302 -1.898700 -1.127370 -1.436940 +5 8.20617 -2.047820 -1.624180 -1.628720 +6 8.67205 -2.252980 1.762220 -1.858880lp__ accept_stat__ stepsize__ treedepth__ n_leapfrog__ divergent__ -1 -9.73111 1.000000 0.932232 2 3 0 -2 -6.75346 1.000000 0.932232 1 3 0 -3 -8.48107 0.665681 0.932232 1 3 0 -4 -7.74234 1.000000 0.932232 1 1 0 -5 -6.82333 1.000000 0.932232 2 3 0 -6 -6.82333 0.784555 0.932232 1 3 0 - energy__ theta p_theta g_theta -1 10.97150 0.3939970 -2.155610 4.166930 -2 9.05119 -1.0295100 -2.933880 0.158161 -3 8.61284 0.0530292 0.702591 3.159050 -4 8.45793 -0.2148440 1.637290 2.357940 -5 7.76204 -1.3633100 1.875250 -0.555575 -6 8.04989 -1.3633100 -2.143570 -0.555575
The column
lp__
is also provided viafit$draws()
, and the columnsaccept_stat__
,stepsize__
,treedepth__
, @@ -816,15 +826,15 @@Saving and
energy__
are also provided byfit$sampler_diagnostics()
, but there are several columns unique to the latent dynamics file. -+ -+theta p_theta g_theta -1 0.3939970 -2.155610 4.166930 -2 -1.0295100 -2.933880 0.158161 -3 0.0530292 0.702591 3.159050 -4 -0.2148440 1.637290 2.357940 -5 -1.3633100 1.875250 -0.555575 -6 -1.3633100 -2.143570 -0.555575
theta p_theta g_theta +1 -0.959636 -1.160790 0.323414 +2 -1.779880 0.178790 -1.268180 +3 -2.089350 -0.136971 -1.678370 +4 -1.898700 -1.127370 -1.436940 +5 -2.047820 -1.624180 -1.628720 +6 -2.252980 1.762220 -1.858880
Our model has a single parameter
theta
and the three columns above correspond totheta
in the unconstrained space (theta
on the constrained @@ -867,22 +877,24 @@Troubleshooting and debugging to
TRUE
. -+options("cmdstanr_verbose"=TRUE) mod <- cmdstan_model(stan_file, force_recompile = TRUE)
-Running make \ - /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467 \ + /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326 \ "STANCFLAGS += --name='bernoulli_model'" --- Translating Stan model to C++ code --- -bin/stanc --name='bernoulli_model' --o=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.hpp /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.stan +bin/stanc --name='bernoulli_model' --o=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.hpp /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.stan + +--- Compiling C++ code --- +clang++ -Wno-deprecated-declarations -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_15_0.hpp.gch -x c++ -o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.hpp ---- Compiling, linking C++ code --- -clang++ -Wno-deprecated-declarations -Wno-deprecated-declarations -std=c++1y -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.78.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_13_0.hpp.gch -x c++ -o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.o /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.hpp -clang++ -Wno-deprecated-declarations -Wno-deprecated-declarations -std=c++1y -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.78.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.33.1/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/Users/jgabry/.cmdstan/cmdstan-2.33.1/stan/lib/stan_math/lib/tbb" /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.o src/cmdstan/main.o -Wl,-L,"/Users/jgabry/.cmdstan/cmdstan-2.33.1/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/Users/jgabry/.cmdstan/cmdstan-2.33.1/stan/lib/stan_math/lib/tbb" 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/RtmpXOnqh8/model-742871ce8467 -rm -f /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/model-742871ce8467.o
diff --git a/docs/articles/cmdstanr.html b/docs/articles/cmdstanr.html index cac1ca1bb..06f9ad3ca 100644 --- a/docs/articles/cmdstanr.html +++ b/docs/articles/cmdstanr.html @@ -12,7 +12,7 @@ - + @@ -26,7 +26,7 @@ - ++--- Linking model --- +clang++ -Wno-deprecated-declarations -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.35.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/Users/jgabry/.cmdstan/cmdstan-2.35.0/stan/lib/stan_math/lib/tbb" /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.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/Rtmp2hQnhs/model-9a14560da326 +rm /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.hpp /var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/model-9a14560da326.o+fit <- mod$sample( data = data_list, chains = 1, @@ -891,28 +903,29 @@
Troubleshooting and debugging)
Running MCMC with 1 chain... -Running ./bernoulli 'id=1' random 'seed=90958316' data \ - 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/standata-742874baeb02.json' \ +Running ./bernoulli 'id=1' random 'seed=707272264' data \ + 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/standata-9a144d3dc244.json' \ output \ - 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-202312131009-1-4db120.csv' \ - 'profile_file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/RtmpXOnqh8/bernoulli-profile-202312131009-1-1c1f83.csv' \ + 'file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-202407021546-1-38ac27.csv' \ + 'profile_file=/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmp2hQnhs/bernoulli-profile-202407021546-1-910909.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 = 0 (Default) +Chain 1 save_warmup = false (Default) Chain 1 thin = 1 (Default) Chain 1 adapt -Chain 1 engaged = 1 (Default) -Chain 1 gamma = 0.050000000000000003 (Default) -Chain 1 delta = 0.80000000000000004 (Default) +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) @@ -925,19 +938,20 @@
Troubleshooting and debuggingTroubleshooting and debuggingTroubleshooting and debugging -
Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson.
+Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson, Steve Bronder.
- +@@ -49,7 +49,7 @@ - + @@ -64,7 +64,7 @@ Other Packages - + @@ -125,7 +125,7 @@ @@ -133,15 +133,15 @@- - + +Getting started with CmdStanR
Jonah Gabry, Rok Češnovar, and Andrew Johnson
- - + + Source:vignettes/cmdstanr.Rmd
cmdstanr.Rmd
+install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))Introduction
@@ -154,7 +154,7 @@Introduction
# 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")))
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.
@@ -214,10 +214,10 @@Installing CmdStancmdstan_version(): -
+[1] "/Users/jgabry/.cmdstan/cmdstan-2.33.1"
-[1] "/Users/jgabry/.cmdstan/cmdstan-2.35.0"
+[1] "2.33.1"
[1] "2.35.0"
Compiling a model @@ -241,20 +241,20 @@
Compiling a modelmod$print()
data { int<lower=0> N; - array[N] int<lower=0,upper=1> y; + array[N] int<lower=0, upper=1> y; } parameters { - real<lower=0,upper=1> theta; + real<lower=0, upper=1> theta; } model { - theta ~ beta(1,1); // uniform prior on interval 0,1 + theta ~ beta(1, 1); // uniform prior on interval 0,1 y ~ bernoulli(theta); }
The path to the compiled executable is returned by the
$exe_file()
method:-mod$exe_file()
+[1] "/Users/jgabry/.cmdstan/cmdstan-2.33.1/examples/bernoulli/bernoulli"
[1] "/Users/jgabry/.cmdstan/cmdstan-2.35.0/examples/bernoulli/bernoulli"
Running MCMC @@ -277,30 +277,30 @@
Running MCMC)
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 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. @@ -347,17 +347,17 @@
Summaries from the posterior packa 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.27 -7.00 0.71 0.34 -8.70 -6.75 1 1852 2114 -2 theta 0.25 0.23 0.12 0.12 0.08 0.47 1 1611 1678
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail +1 lp__ -7.31 -7.01 0.81 0.34 -8.871 -6.75 1 1404 1631 +2 theta 0.26 0.24 0.13 0.13 0.079 0.49 1 1130 1314
+1 theta 0.26 0.13 +2 lp__ -7.31 0.81variable mean sd -1 theta 0.25 0.12 -2 lp__ -7.27 0.71
-variable pr_lt_half -1 theta 0.97
+1 theta 0.96 +variable mean median sd mad q5 q95 q2.75 q97.5 -1 lp__ -7.27 -7.00 0.71 0.34 -8.70 -6.75 -9.165 -6.75 -2 theta 0.25 0.23 0.12 0.12 0.08 0.47 0.065 0.52
variable mean median sd mad q5 q95 q2.75 q97.5 +1 lp__ -7.31 -7.01 0.81 0.34 -8.871 -6.75 -9.443 -6.75 +2 theta 0.26 0.24 0.13 0.13 0.079 0.49 0.063 0.54
-CmdStan’s stansummary utility @@ -384,7 +384,7 @@
Extracting draws# iterations x chains x variables draws_arr <- fit$draws() # or format="array" str(draws_arr)
'draws_array' num [1:1000, 1:4, 1:2] -6.78 -6.9 -7.05 -6.85 -6.75 ... +
'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" @@ -394,8 +394,8 @@
Extracting drawsdraws_df <- fit$draws(format = "df") str(draws_df)
@@ -403,16 +403,16 @@draws_df [4,000 × 5] (S3: draws_df/draws/tbl_df/tbl/data.frame) - $ lp__ : num [1:4000] -6.78 -6.9 -7.05 -6.85 -6.75 ... - $ theta : num [1:4000] 0.284 0.186 0.162 0.196 0.252 ... + $ 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 ...
Extracting drawsprint(draws_df)
# A draws_df: 1000 iterations, 4 chains, and 2 variables lp__ theta -1 -6.8 0.28 -2 -6.9 0.19 -3 -7.0 0.16 -4 -6.9 0.20 -5 -6.7 0.25 -6 -7.1 0.36 -7 -9.0 0.55 -8 -7.2 0.15 -9 -6.8 0.23 -10 -7.5 0.42 +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'}
To convert an existing draws object to a different format use the @@ -426,6 +426,10 @@
Extracting draws$draws() method, but in most cases the speed difference will be minor. +
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
rstan::extract()
.Plotting draws @@ -453,7 +457,7 @@
Extracting di
-# this is a draws_array object from the posterior package str(fit$sampler_diagnostics())
'draws_array' num [1:1000, 1:4, 1:6] 1 2 2 2 2 1 1 1 1 2 ... +
'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" @@ -462,12 +466,12 @@
Extracting di
# this is a draws_df object from the posterior package str(fit$sampler_diagnostics(format = "df"))
@@ -487,7 +491,7 @@draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame) - $ treedepth__ : num [1:4000] 1 2 2 2 2 1 1 1 1 2 ... + $ 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] 6.79 6.9 7.06 7 7.9 ... - $ accept_stat__: num [1:4000] 0.996 0.984 0.988 1 0.835 ... - $ stepsize__ : num [1:4000] 1.03 1.03 1.03 1.03 1.03 ... - $ n_leapfrog__ : num [1:4000] 1 3 3 3 3 3 1 3 3 3 ... + $ 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 a [1] 0 0 0 0 $ebfmi -[1] 1.0 1.3 1.1 1.2
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.
@@ -496,34 +500,33 @@Sampler diagnostic warnings a suffers from divergences.
-fit_with_warning <- cmdstanr_example("schools")
-Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
-Chain 4 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/s0/zfzm55px2nd2v__zlw5xfj2h0000gn/T/Rtmpl3ZPjU/model-76ac71f3c1fe.stan', line 14, column 2 to column 41)
-Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
-Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
-Chain 4
Warning: 149 of 4000 (4.0%) transitions ended with a divergence. +
+Warning: 374 of 4000 (9.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.
After fitting there is a warning about divergences. We can also regenerate this warning message later using
-fit$diagnostic_summary()
.+-diagnostics <- fit_with_warning$diagnostic_summary()
Warning: 149 of 4000 (4.0%) transitions ended with a divergence. +
+Warning: 374 of 4000 (9.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.
@@ -562,26 +565,26 @@+print(diagnostics)
-$num_divergent -[1] 11 54 32 52 +[1] 269 29 11 65 $num_max_treedepth [1] 0 0 0 0 $ebfmi -[1] 0.40 0.36 0.37 0.30
+[1] 0.15 0.34 0.38 0.37+-# number of divergences reported in warning is the sum of the per chain values sum(diagnostics$num_divergent)
+[1] 149
[1] 374
CmdStan’s diagnose utility @@ -544,7 +547,7 @@
Create a
stanfit
object+stanfit <- rstan::read_stan_csv(fit$output_files())
Optimization
$optimize()
. -+-fit_mle <- mod$optimize(data = data_list, seed = 123)
Initial log joint probability = -9.51104 - Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes - 6 -5.00402 0.000103557 2.55661e-07 1 1 9 -Optimization terminated normally: - Convergence detected: relative gradient magnitude is below tolerance +
-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
Here’s a plot comparing the penalized MLE to the posterior distribution of
-theta
.+ @@ -592,18 +595,18 @@OptimizationMaximum 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 = -11.0088 - Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes - 5 -6.74802 0.000938344 1.39149e-05 1 1 8 -Optimization terminated normally: - Convergence detected: relative gradient magnitude is below tolerance -Finished in 0.1 seconds.
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.2 seconds.
Laplace Approximation @@ -621,7 +624,7 @@
Laplace Approximationmode 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, @@ -629,19 +632,19 @@
Laplace Approximation seed = 123, refresh = 1000 )
Calculating Hessian -Calculating inverse of Cholesky factor -Generating draws -iteration: 0 -iteration: 1000 -iteration: 2000 -iteration: 3000 +
-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.50
@@ -652,40 +655,40 @@+ 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 5e-06 seconds -1000 transitions using 10 leapfrog steps per transition would take 0.05 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.262 1.000 1.000 - 200 -6.263 0.500 1.000 - 300 -6.307 0.336 0.007 MEDIAN ELBO CONVERGED -Drawing a sample of size 4000 from the approximate posterior... -COMPLETED. +
------------------------------------------------------------- +EXPERIMENTAL ALGORITHM: + This procedure has not been thoroughly tested and may be unstable + or buggy. The interface is subject to change. +------------------------------------------------------------ +Gradient evaluation took 1.3e-05 seconds +1000 transitions using 10 leapfrog steps per transition would take 0.13 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.27 0.25 0.12 0.12 0.10 0.49
@@ -697,31 +700,31 @@+ theta 0.26 0.24 0.11 0.11 0.11 0.46+mcmc_hist(fit_vb$draws("theta"), binwidth = 0.025)
Variational (Pathfinder)CmdStan User’s Guide. Pathfinder is run using the
$pathfinder()
method. -+-fit_pf <- mod$pathfinder( data = data_list, seed = 123, draws = 4000 )
Path [1] :Initial log joint density = -11.008832 -Path [1] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes - 5 -6.748e+00 9.383e-04 1.391e-05 1.000e+00 1.000e+00 126 -6.264e+00 -6.264e+00 -Path [1] :Best Iter: [3] ELBO (-6.195408) evaluations: (126) -Path [2] :Initial log joint density = -7.318450 -Path [2] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes - 4 -6.748e+00 5.414e-03 1.618e-04 1.000e+00 1.000e+00 101 -6.251e+00 -6.251e+00 -Path [2] :Best Iter: [3] ELBO (-6.229174) evaluations: (101) -Path [3] :Initial log joint density = -12.374612 -Path [3] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes - 5 -6.748e+00 1.419e-03 2.837e-05 1.000e+00 1.000e+00 126 -6.199e+00 -6.199e+00 -Path [3] :Best Iter: [5] ELBO (-6.199185) evaluations: (126) -Path [4] :Initial log joint density = -13.009824 -Path [4] : Iter log prob ||dx|| ||grad|| alpha alpha0 # evals ELBO Best ELBO Notes - 5 -6.748e+00 1.677e-03 3.885e-05 1.000e+00 1.000e+00 126 -6.173e+00 -6.173e+00 -Path [4] :Best Iter: [5] ELBO (-6.172860) evaluations: (126) -Total log probability function evaluations:4379 +
-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")
@@ -729,22 +732,22 @@variable mean median sd mad q5 q95 theta 0.25 0.24 0.12 0.12 0.08 0.47
Variational (Pathfinder) -
+-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)
+@@ -768,7 +771,7 @@mcmc_hist(fit$draws("theta"), binwidth = 0.025) + ggplot2::labs(subtitle = "Posterior from MCMC") + ggplot2::xlim(0, 1)
Saving fitted model objects
+fit$save_object(file = "fit.RDS") # can be read back in using readRDS @@ -780,7 +783,7 @@
Saving fitted model objects
$save_object()
and replacesaveRDS
with the much fasterqsave()
function from theqs
package. -++ - - + diff --git a/docs/articles/cmdstanr_files/figure-html/laplace-1.png b/docs/articles/cmdstanr_files/figure-html/laplace-1.png index ed19be52582cb48cb15203f27c1df33b10014478..975e63a2e8cd24745093931bd5d884a9933debd8 100644 GIT binary patch literal 38787 zcmeIb2{@Ed*f%bzv?6# 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. @@ -796,7 +799,7 @@
Saving fitted model objects
+# Load posterior draws into the fitted model object and omit other output. fit$draws() @@ -891,21 +894,21 @@
Additional resources -
Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson.
+Developed by Jonah Gabry, Rok Češnovar, Andrew Johnson, Steve Bronder.
4yb?EAj2lZvt>3Ry OEWH<6jQD}VJF}$jqWr|=6GReOD;>&nfq~ST9 z{^eLhW>K=s#2?@H2K8c8KQf&*N}ljxf0upx9f9Zh#jhx-w@s7>#|n6|ih}ZuOKK2H zgAr41L*3Ss-G^V~v#GxBK8w$GB9(ohxi70~u)JKz=A-M^qpn}sS^ZzJygIk%WdA2v z$HWwUWp~>`9ulKJmU>%ik9XIT_74Kp41KzItK=1#>K+7_(Vq{tJzP9ke#s<)zVCQh z4xTbq6w$BDhOod3-kR6g&L{b 7n9$}GB%fZ zqTO?^bs3j;dddZ4ak-N9?9E>YVX=Ds!;2Vt_b0-}JB`CItAZG#TVQx_nh{UEsQ8K1 zp=r@_?Yqv!r-p43llU9qtO+e)`ouk+4z?eKt1BlZmUFl@BM#ZRDMTC@B%;OZ?PFn2 zcS17Krpp+zc-QzT8(gfK@m|ayq|n26on+TH6l(h_7^nCAGTWq2{|t#zf@ryNVpPbL z6;1X$vdXg}Pc Y=^p(%<$~jE2#G1J_?g4#> z&oZ4s)q9{Y=!t>bL3BYAYNF-jp?+V)KR$LGC+E^4m3To Wd~F ub`dDTBL=#Yfy64U%oY9@Mj7Y~jAndHtp5 z;wgrppIp?ZW9VNpeNoi$)oB#Eag49(HG{RzY(Zm-;~qu@D8Em!?XF%2j=x(BABZ ziW^XV;`WC22IJ~iZDtxmu4c+UMUqO;qE_n>n=|?+jM!~HceF%WHDEqz4z;Cg4)LoT zSv}yk*L|#0DUWZAw^r$pYoWPnxzwzPtd3|L`(2X 1yy?3ml3F_4lhT z|3rQfdujXy`I7FOdw=&&(I55JP%EfqkH=^AUJffwD@{u}Yp9e_i{lip@nl*|#-TyO z5P1ajib^}X@OhmlBA?1W>D-hXCV%?)-~&!4JN{%+KT7qY{Z;AaI?1jdJSY>b$IU)u zwf+#hBl&tilS_R66wh8+v;N17*TtwAo|qOnJ)W44w{yEG(vb$Kw4OE*-uKW&(}h=* zOze||Cv}fcsr>kKi`7VG&9OdT!^rHMH>1zXzWnTecqr`x3&oKEc=Y+P!r`J*jq%5F zLJS7*(w#2|d@{ejx 8V0N|}uid3Y1A|_1a-H4V(u0>DjoFJi>u|V2yR`-Y zgn;h7>v}RWcr4(W5O0r{Aszv^+5`Tc0DtlD@T2k9e%XhJ-n(|a7mB^f8}c>)4^JHL zy2Rx>4tvJCWPH_!?jQJttoNm uRef0r>{G3#;>2(QrNoSqte4VLFA$^4-zq3JhY }M6yxH!-K^1w?p8Ao>xk9mu7jMC9h7# zUxYaK?PvDetTQcV=q=I`RgKPydx*#3OL@heGF#Pgj5HDTN0$cttZzDIyV=+v6@-Kc z%H|uCfGTkfh`#yjnU);k`4a1dVE*oU#ef@I$bsk*E`Dv%ooi&@a^tW0YU%BPplgp> z2ybuBY*)p5OMOKP_icZ>mW|#5TCm5TX Fvds zegur{0&19I``cVu9&p#e{Q^;&`MJie*Qe<%G$rY}9&LYHkAgVA0fARVy6w>57YGO@ zYldS~*IL`mJ9sTVJ~HfNufz5Tc<(6~8MUj;;q7mS-WTr`wcdM5ZhHu=Atuh}j0rHi zy!~yS!-NQZGmn0tU)%6piyDmVzooF@jkQCU{>XY s&nyi#{iR4kGoqW* zePOUo8T0!6nYgWt2_PZP7u)|?8!L@_QSg-@FME|&`pR+38r3ILt9#5wkLM?t_`lmM zVsV%Uyrie%JXY;_=`GTiDvYY+20@TuZ~ZOi#1_G76$616zY4=DYC0F>rdM3?uUn}% z4 E8FF)9`Sg)= zAVlGH5sR{7Yux#*VuBL@jX!wX^g1vs*F&e&m|^)UuiV>*!s;3Rk@xT+LPXIyNnpXZ zQ{WFWHOcI(RM_J_{iAdbfJTG@{3js%Nwz3pf>|E2`nFIqHIyW^#nM5a%==!ta{dff z&f(w8gmv2wl?ARKG9~h#UVbZc^l8x7)jGSoTj}~^C;v@fhQvIb98CfY$bSJp%$Wwu z;TTZn$o|6X0}a(vE&3NX*U`a1MT<07tcYcx>W!p@<3BuuqyS||;iV;}=kkId*1V5t zL8UK`oYD^vm>m+0vR*M_t_l~suQAolJ}?=Oc_-%B))vHzKm&HpQ-cHwha!TlX48Wy z6)BME?`uM6A5Y#lTK-ig(;F2>wIyByl0N}r4Ss}+GSdt7-ofekPjyokUuIB0A{8#e zlq!Ax1wHXJ;Cn(U==b1L6GXWd275|C&_z?G`1a47pStt5{G+ iagEs zuipV)clTvn&Ovq=hRY7D+UzRHR@2{Ak}VLmt0e#8>;F(Ckt(QB{Nrs56nlH$yG*p! zjO8^U+|)~7#8W)kJn0AdI(%fTpF9PYd)e&iECX!{2TcymYz=ar0^eOG8udG z0XGC4=Pv!44URv)3CKgGL3W4agXy2_4rJw!%ffdBj _FBPU`WIv zU)CGjd=afF-rGI(F!*2Q{{H1LNVuKm2--fW>w&wn^e?}O3BmgDRkd?FnXnt63g#of z169PwiSxbkZ|$_2u>z=aqm0>1%sL0?!AdldW!O%uM<^h;6%NUrr~+6=+jk*s2bRo% zRZFzGFD{<0_Y1O&FTnE9uAp&aJIJn}H|I8iV(nV=KP=g`=zkz|7rOsp$^Qzvd(ZQB zt~*MIJpS0cr(rk-X0=i?rZ36fF&lfyy#xgnp&Jpv^*d}TVD)==q43lnx6Bm+KY*l- zaIv_pQD_%e%&Saub*OtOWnO-%byUZ|Uw@O=;nIMo1yh^rp4;rw*Zp%lOYbE_UX#Nb zB r_k>siUu`L?xBV0s4{{lhN+000iK%PMq8HbS0JR(uLHBpM zeH% Wm4IR@gSV6RpP}?jU%{Okj8VSFN9d|1tZFX0ru{*+>OhIs2>-b*zJl zriGxHpsh}0Ow~QXon!L}jDn>dcjoNPi!)Jn{zi9j>z*Ogp!;bhz2l%348(tGr8aNb z?Ftv&d&u*r$XfJ3c=KMO_h%r?L_MW?N{`KXmLIRDv*j#UH9Lj7 1ym@US5;!H|wMC!)SV(kQz$fn>J zVqE-tyzkg2RyG_tZq|4qK*X_UxS=7gqa_Autg5;XJH1T<5 ooU9P`SC8K0s7{82JtnW;7#6 zuf3q!i6~1z6r%PJE)Q|~{t8Ns?WK?4{DM0ONIpmwnd$AYA%S=p N4v zd#hUz-v#Y0O0w&laUCxfQg(gwR*w9qhG1}cmuEP!pB<-#)y$VDnA=$L7%Q~6W~SnD z!7NJi!z#;7LENfydMJq9H0AC{{Ea0I6na-jJdNo-RGS9}V9q$F++3QHxsH$2ag{%| zP7mZ32sLkY-9Ehbi+Dly5@!ooIwocchY{u@<)fVLT-k_D#Wg^5+1>r#_ixbrH5TsI zuOK>mCn~we5f }OW6C2w{?Qg@4y{e2x2;=@bc=`p|BJ_erw`-c8YU%v4$ zvT^M!fd5NL*d6xz1m{1fZE*Q(IV>BR3h&-G={oo8v#lcvRLh!AmUxgHNR_KjHZLwj z)q5w# |uIaxKNE>c%WSE^j?BvRLHhf|}m6{x>P!#@r-2Sm9GwN>;u!-COF z1x5--*GwG`@_L_(fl*DDz5Ujx;Ey@Bhk;@12VdM;+50EV1aA{evfKA`ZC|yZr4NuL zQoh4`d(&0{`?m5BuJ$`@ZjS`C+e;cUZa**su* tCSWOaZ2^JPRyH!++C_O;a5qk8Q8wu4G&@nk~+=# z9kv*JF&N80#R}aemp&uKV9Irj!MQL6*`v;3Yo0t|A+~` mNb? z9wB%e?rO %ezs&2+ssktQk${ zTz6T#=PkgV-L|xWg2m3MLx@?Ms~ejy&}6^fCZb#Hc38P*Qv#I>?X87W>-}%K7zbNo zOXW=Uw_kgQKLkF-$27%fhufmy_HNhfHuC#om*Lx(-G5s(>ov&Wx8~lkb^<*T0TEY+ z@gz0Z$psXp9&;bJSLu?*QrnE*H*`SQAPx-tTL>CJy4c(8Nu8`vYLU;< B`EQzV`8>#-|=%JF=wg5ZX(YTJHy zOoQ6|`=@JKt92isrS5#Q>4w;1kAvOR8q}N{WTg-;zPgZN(;>10;wpd=o!H`wJ4BC4 z0MB+`UKbCz(As6#hPwS`3}jbHwi@$)L`i +L7oPZCB%rFC14m_+w_Q~6V_>!i^~ z&*(rNzu>ibGpBTSR_vQ?kvXi`X?}Fxb-vzp${sc~*Htf4GgT^zY=xP1Nohsn<4^7# zxGS)gnm O|+ h+*Pfi*lk0ZsI(#d0UXcCY#MY6`zPPEmokR58XQ zpZ^1f9ku&+OC@W?~Cn)e5rAtO%nI4|mO()UKKFgJ;Bs~^PRbq{CK`j;2i z6KiWC_$?~&1;EztP==#BAjUah)cm`hZ{pc+z1n5_uEMPW%wG)LHMRfWAajfFZD7ad z?GH7aTs+q4^@~5IshLV`7A8GaXYRc*Bv+ksxoUiqHl^`z&6gHp9|59jg)|z2gNelV z#I^Jmdi?>qSEYC+e`l+Jn!-+>VC;>mqi+0DaL3h+4R%~t=Qo}DDsu3a48=9p7;kjx z$h**# nW(?{EZ>(N zoH8_VS&5vuW7`PHxOjnMdfS-1!*P|8*MW}=u8%MM7k#k5z`qg@AiN~g9&Qi(RA)is z@ ?ze$Y>cE@FnL5;$wwFl65Vx zhto)<;;G~#Bi`-&6eA-D?=B=*ISL&QTwr~!lws#+q_} 3io6>q-X6nvQOXkJiq-fn{0gysgWpI;xfI6q*Ana~3J`L9M5J_yDz zoJa+SGg>g 3!Z%LhDZo5y>?9f+duo#a_h zgKoJn0u86_ktf++i^Rp@UHY54!g;QRp1xOO?>mO=srB&iC~8<|(Uyz%mqAFkU*g&h zyx6nH9=Zb?=x=+1pyR0b^PPBsKPkFXX=L<%Ag2?= 9Hh#mEMf_7jtG8!q>_lJJWb#;3kZB 1Ao$<5~t!MO-p(^#Ktw-Ln**4w(d>&qkNS60yV0tMDl z4@^rY>fg0z+26}a{&Q}xuNMF$hY0-=+~5r^edtzYk7=4(pj=!0IilRN33hQxIVuM_ z_5KFLRK7sPk5t}G#8DF*8FT3YVAQtzWV^%J!27Y5Or;#ReLM*dJqTuj-w(E_Z?ET| za}TtZ+N*6!zOgE-eGW9P$E0Sz11fQWv3XkFWx?Si4i&y(@Rrbd $0t*4qg;&hb zz0G|2$_2U0>eKf*XNlJ-*JVx=dd~3Ap&G-g@R0ZgnV}t`EP9LE;_#K|CucV|w+ lqSUWDvr9%4FrEaijZKOYD)4ipKw-9u45|rOfut zK@9;RTLS#2j(kMn13_c8;?$F7%+DPneg9bI99!>AMd$dJfw+RgQ7|?z_R@4e5uBfJ z>{dKlA0TH`qRvhtKMZ8bX`)TEz3mVW-v@B7&qBSKG3!Hte**e@xU+q2!oL;TohX1o zgx!$bK_Ao8;Vqwv#@90pYp?sJ!N AJ98YRoxD@15ZjvmJ$fp`{kw<5Bo%3nnVuoh(Xvb0XY=#|=rh4EuMx+90JD_2U^d^gQi}%M(NT?YuJELjq(U zh 7k=jLlA^7gp6UQwYrhgic zLP$tJ;2KDpcyg;9K-37Du0BrjEhRAPDSa2}jOU&!@DO0Z6x2sChb9k}n9tp!m( 3p|R1~QdD*Q`FrHKvIJgv^3LsDFZw&k|bw2d`x~>7Mg@riWNe z0dAQ8(|6JeMQy&g@pwj!I}r{Is?Vo%+PIF6-q$jEw|H+kY9vsEF1?monXDrB{mkVW zwTbq-l;#Kdw_of9I_bswqLB_+YdbKrh)j7U8%bTmf2XzTp=FtQeJQ 6A${; z=ecFNi*cMbbt7epxvV+NIqV_Kt`BN!mqRLlI&m~Z7iDy8E`+hKG>4074}_2opagG1 zGa^ 7X4yI<_IYQ0S- z^;x;pDTqR;^^4Ba6ejapo7=s=nU}LyFU&7Xw@9JW8buwAnpI1+^fViCZbyHbS}?e3 zmKOZ|xT%nR`$d{j7<;bic)DwQxPJiDwg2*raKfUXL+6CKf#WAvctRDc?*1b#`4!5X z@4Qr%BC8VU-dRN`;IS)bsdc_G^Nz3itgXlN>U)!%bgG^;%wa}Qwdo*VJ;yU+4(n8? zd;fS-LmYPuQP!=gIt4ka6_{a!I% X|WaZtCGOdAvRt2$M^Mmjpl7a+B@uwapoT uy=D+)7?_<1qDwECf3cb*}U!U(1z3m%foHiuXxuZCN&dl;FcjPFj(w) ik%M)W-Rq-iL zdp=&RgM8v1B{ajUE=NDNJwm=nL9zaME)6v_BfKx`H6k;#FU!`%WtEO 86YD|6B=WWpGH|=;+u1SxZhDzV(p}ntun@;2VB}oZqT5_?y9}HwD zHGP)43bNTn^_i-M# ( OOf6QOka+Djb 3qsL z-uox9&L(<4#0CR(J0{(Hlpk^XmYiJM?IbWKymt-~^ezfgs>Ghhgt*z3R`F@5`qjBw z)f6SCv>zRzQ@)N$QxV3GpxC?)h;4&{f<_*$c(klOd-lv&`*j^sHTL4~c{Es-3W;sf zUWW1LeRm1?D(?Kj0viWN_ph#aKA1dwiSgLn(!KgYG DSnn zP{WGn 9TvYNWkJaCQ#SF~_V^7ZZXg-+HdS`})~N(R hYn^rZI7bZRO?VM I zOuLxt-^ZGpo9UisxKDPwBPY8q2@Pd#nMv&svuqcVY(GexC Ko29eQ6 zhdBg{&j`hkg7HvU?$KWGhH7 sx zC%K*Ptx)+8J(E37&-D_&M$k;d>dRRe334&{^&x8HL{b0HBe6+4wXU`xGu!V{9`IfQ z6;D 29O?}D|eAL%aGMbc}I`wksg}VG?a!f78 z=5ysSR?Gg?uzXARTVCoHfl+Ch5InU#dLA-%BFn2QAX97*FxB8W+$tcPQp4fLSw%%f z=tT1ns=;j`OC56|Cw1&Pbix|c-SHjty#sd@9iAp4iEgJ7r7IeEt?E{I;tp*1osWNH zM8vYhxbq!?^P(lL=mjRzuEP1KvYZdt5nn>A3SyJL1&5CqBO+iQVL~pl=l%Ln7f9RQ zKC>#h(34|7)osmh-n~3fo$K#2ty9zSEm_%}r&=&o`)PIM`zd{f-c?BTIUlU_Ss@S9 z_ulwfSULUx88aRr2z2xf1D7ulql-NvMl69YCf^*xY1y;-{E2_cm{aEMrCd__bQgTa z*&AUz^#a Y(9PWO0 zV$A+G%g6FG{eiZLD-t~AGmCOv_&(OBGSg^r 7gh5*`$p+3t zGBA|Q^&=xAM7FiOunhmP>v=O1jN)2L1d@ZmCszN;Awf6{2zR;ri8U`pKB4O8e@W-R zXDTyWY=1GGkN7eZWpz>eLg83SQE^3wxfqs{99V-e%h#l{gOQ}MIm^57e}>vw`R1yV zwS|Mua}jj>T@5}=s@?Zf?L*}xzIQV8AQ^gMgvRrJVXc0nK_J;h|DFegq{y^70zi z$T+dmd3b%~!JD1{^c>V1J=ZYzju*pIR>L^O&me1y;izh{5*bk_RR@#y81bc+t_w={ z-S389G! jx~ zY1)E%?w@auQq%5v#f@$-EQshkxAM}epjp_>!9m4Dj=vw)V+b?FqFTD9mUo-S-h^m> z?n_>7r-onb17qbC*YfTkOe-(D*f7G^9(BHd=!t%;<0*UkW3O5}6NFNkCUjMp#I+po z8A&WOf&Wu{WHB$?pV$9veonM&Ujkt!My7xVpDrK(sS~M|%BibiUt$}Q7hL)oe9H=Q z?hsIa-vm_-zwn_6%dC~b3Isy~omAA+of(mLnJ& h2R&j(2NaZA^gUf?vie z`?1{U^P^N**27Xca<4&*@uCRMRarMD6lUQlCw#kkK8NdW(IV|fYlgu&k0n#}faHEB zg_KS(doLrWdgkF=9JCTfUO9IxWnT~U)l3C+8bOm~9ybtrPB}$o-Wp9AvS9uj6PuhE zvEMv;!}MDPBHS-Fg qBj<}j~lFYQNfM$d_7XVs7BbXE_T%^H;wnk-+*D~??Frf%LfU~80x z6&=-C$OHfIX{NVnrRd3`p3Y=_HT_bRFX(~lnbSUjX)Yoi@)`XHGorw^7|Sf{3VBfI zFZtB&3LgS$k!h!z99cmV_a4Ey()Q@`jQwWocRK8$!n~KGAM_8Z!9~=XHB&hjwG4nV zor}0+dWK6D-;XF?>MNt`TeT&ju(|NEDEvvwT8ImGzQ{~1PlH-cFsbGiL0cnyneI%s z7ut?RX>|2<1mh1|E2>$ivWFzkE(SGwm7{`^0fhp33{Cxg%wYi)MAYH^&yeU4bS)4t zmp^`4@%}SYzcEhxl-MQIEf4XmeQz$9ReD|s@I@nehH@xfB5$d;1?fTWrZ90i{T|@2 zAF-o`&-P3#j5>u36eXkmMHbrqzRg#>9!t%eoS`s=tpLgm3dMse>+&x~DmqnBI=d&S zrPedbp?p&8B`j9WwcSi>FYTXW>~YL+zln$t%XKQ&1)T8>f~TkLvGx02Hd!U$syFuv z57J0%P@PXb5rNKZxs8!6em>ZIc%;o% Fp{>ac>}*e$ILV@0Jw!1`me#3JS*Rg2^t!F7vf=`P&O_Z$f`?ARf@AC zhhao6Bla$Dk@Ci^9~>1<&1MTMBxa;J)rkt3_4QE76^faxI$sweH#`Y+S5ir ioY4@c}rGy^>S0dSbb4y;b>CbO2o6eaDH <^1DfF+}(qxYlqP$_m$XI(N^F8jN|KLV{uMxJ5m%k2b399 z(M3~ZN$L;5x1KT&$aOOrYfYhng%~n-l&HzQR){JlAUvWGJ-r7>kJ5X^>=u-tTiC5D z=qjjgi;}gC^ZKJtw1xE8h0YV~FXn;4-zYq7ma9qfDX*OMf!k>1Npr+qFS}yBK}PvW zOk%LDneKN6@Nnfd`Y~g1$Ky^^StlF+21B34>BlT*`+PXv?je%J&Ind$=DWS;NTW?< zvb`I33ZKA;IQ0En-?tPbFz(`l4iBoYt#{js*yxU11oTgizeI`_jbaUMvRUgSR#p%3 z!vew;K_DP47ehJSy#VzUJ5UlQyPB_860gu*cG0vWJ2TS^t7nXyFP8Xvi7{r{c0njh z>4Xo!w`ob^7X*Zpm~d(`W$}26YA`A0ja`qNKktJFDWjx%;^E8=>+-#feqlusW7^O& z6*R9aI8ginX;H&k!wr^a;qMRq6gI02Vs|H-q^y*gZcZ(psWu)HJ@7jmmLKd&MC2!G zX4p-1oOLC--Zddl3wZKcm43)ryy#jUchU6~gQ*Z(fB{FY(jzUij$2?l7X^D?OdS(e zOEYT=Dr4^iUxIC61&V?$bwwXL7%XMczo_B;@P+8F$Ifg$#Cw0Fl$=r?OvEZLYo{UB zNyw9bwng}Aq0=gnkBYfc*&CjR_qm^l#^f$#=4ralf6EcQEC#>G`UWs0hw!J<(cz{7 ze{H+H=*Bo7J8bm&3tuoGeEP Fbt9Lbb~3)V zRQ`@xH0?w5Dy7 O9T}X~^xrgZ+%~up(AfhwTqtSZ>I~l OL*+uai{w}CeSj`t z2SrVgNyFb?>cY`WtG&w9_amIfGkw@fgj&KDeE%z24k8a$OPm+R9dFLLBa1&W){7t! z SfC(}A*%CfI4%`@=u2j^3+v%-}9!ULRlqfz>H~ z&f;I6yee0g;vGGSj(&QWG4_Ri^bXjwv~ti_0xq%|JKQ-f@k6to^i+w}KsA$1zqi<| zAS@FjmuvEU# {rBxn$(zggVOoUW8wP^$}Hr4~3h_Tct|8#fZ Z*B!JFmxvX%xjX9$HvVfQr)tSBXkUa% z^m8Wp2ZWEgdyjn)tDOHq1Tes+{7MtRJXLZ>Id9`?QpSd{KFSKpSB^0eF-mC~E7tP$ z9qkK^4mVbxT4|z%S6(Tpm1&SU!l-OZGhJa=0|ic+D|L@2Rn~!T0hTDXVOD_lXSv>6 zW`#ia2wksQaa>wmaQ7d2FgY2V;Mz$KLRs5;Ae0rrXPkcQ{I0XuaW+KT!}-37|B_;V z0~-4Jb!UP@(l1JFOn_AML%CY0TOD7UiFt#fdI0=$7Gsb3H7D`d9s$&FjC*1s?6*#L zDAiy!=o;j*U$~DSW@Nu(-5WkSqUNkV)UTSVQ`WWqX3VRo(RxBKy>RRyYl}8!d>>Uh zZeymBWArUq8;$W6YUT_`3uvobZrn$ymR^d`tk*`Dz4Hx_ic;0?Y^dP$KnAY|NcsDZ zJglX?BUU-h^R@wy?*)&ptgCR_KO7sb!AadoCkjg&YWs4eWVHJ5}weHlNfMsuLo zXnBwy<~U%O5ZU@%yi%KiYWk<5tl&T=e~-(onmFTrIHUH;NXANC;JLmppN1sUwLd2G z)+k9j!pN35&$nqWwZ0!ptC{Cb!vb1z?l z7q=*XIq&;xI!SQpjD8z$MZzifyRJ9;U7X^krm z<+lGGq3Gl{cfJJ?j>%?yY<6N9qz!C`t~e%l9 dB1aVKTXLF@_7K{@E=4J*1cYr8FG} zUO}g9v?Ax@XMoRF3^&Ev{rQ6}f@XFnG1@%`KBX*YHJ&A@%+Zl5KBT8TKd& $aT?`FV{5}kAu%ks z+wu!qrZE;G4T3`vqVr2Z|2^G!bWBpyhngz&iXIOFjSi`YA4eVFPwSoy=kDcHe z<}ONg^KB~lMhVQ$hV;ELLana>Pf>m &itTl;=txQV=IU7*b!RMdr}Ozi zc<_Hzx@~tnz|j5^m|)-mm)P6!{uxf&N5rG~JT8!f&m8~lQVT&K>wip4ut~1Ye))R; z{%9D^)3|yPO$V7HLq7@BARW XDF^i}STDg}Z}0EL@a5i$z#-N^I|-J(69a1S zk^|!Rcc4)Gu9ohH<(sh=+M@>XC-0udqWNF_tdxdR3i&mDX 9bv7R?YG7$Y(1v%4^dtxa$u8x$kbu^Rn(Zuq|=dm zna}tXTKN3Uf%OJ`Ge98}71Tw2y;?Cp0N($w#&|NzHk;p98!CnN9bNw>>OZxfkXiz_ z^BuFR(m`Q4#xLrqt3-*K6u+;tMw 2aT_vOD#NjN#g zl9h9MPG|-Ua7v=C#ESXYNS)cXG6DF`_wbl^nwJ>s7D5$X6uY7l+i#1k2F+48X$F-t zG$A_|lxMri(X^zAZOr?|dZ-LFWAD^YN+atzm1{~Ea!sGN7$(IuBqFpKOL!EFYKRzC zkkTR7Zd^t(NutaH6~2n)jIjpF`swUY??6H>Ee92|6OBtqPrf|H{>oK1YL3hv-NM{k z%5q_)M_Si$ =lRAEA93c_CB @k)B94JE#0X{su)uIvKbVd9*XSitGUoD zyszJ#T>*VZ3aTvi7V2}SH3WY=G9E*9vPPJ{awaVD!gtB77V<}jV$QyZlH~L(GhBT} zThvBtfe8_qj0ucqEMaY84y YFYO`#7}V;e3e0oxdr- zCgVn$CbQ>PLpRruVnb!8Hs@+&m^8C@mubhTOUb(K%&p2C7UE*A8l<{$%AbR6V>M}G zO%KaJ{b+Q%$s(p++rPblv;F2PA*8~r^8gzl{|uE<6a?nMwT{X;zLVvF?*rzcjNa zoTWcN+jCRYTn*TYwcG|oxaIatnl+<4?4oGKU|>1$?Kcfngodf&(KiVRokb kOF gV6GE&)v zv$_X&-3PkqLB8+KjZNsylM*U?xsqRPRp21cChX{Tc3&=?^5pOZv81w68Jb+%juvAm z-N^%Dmzq`>RvZ%^6`GJ*Og++GKUo6g*$5EfuB+MW)#_Rn_$B+HXw#4W0_i6#G&|o$ zlUZ;k9H_Vg(-d-2y8oGBRs-!~9voCEkke^%Yy8&bB_le8)o&wo8`cJqCj$$E6 o8zTICsJ&H750rWCE(GsH fsNf>Q|zJBepMnN-!mgCts_4i?f8MT`Od3$_Z?u zjC->V7=*r00HocE(2HD`>hRTQ|6XglZ?sm@^M@q^d_|y_bj(*+R~8(v@eT|zmIx>~ zy;+S6toI4-*2R33Ci2IzMy?r9aK`BV+;;7-UF_`9&9lK(sckyl2wsB9cCM=BGf^{w zGR4#nNjY7MRb+5Iui%uqax(DNIza@L*mc`+*@9#1sLU&^Wg3? vMAo4b+hKYryctE!TCx>c97iDE3W z$3~@@+oR%#t{pdl#{60@vV#kNlCFT+?``X+4nfSYf-8DunUBuRA^#PeJ0JHUz7589 zalUDY^xca92_d`EbXdSr(NIVo4sUa?yoI_&t^@io!`GX?3Tl4jjH6UknPIeJCp&ZJ zGo1`|I-0IlYM_=o 3Xa)V~!l?DyhiT$q4 Xa$jnG?zs6T)ZIhdhX{I%*xzltH z(X$L9Yzj`Hdakij*PPU2of+a>()_Z&w#Y%JygyZ=YL=28s`9-Mr5 kRn@P?B; zLAf-7BC*-V$z(Y}dZH?}cb4=L4Hs 8|tR9)UC!Kba zd6tbAbK|FD@A73Lzomze(L1tELdzMmRoM-)en&@c-d0R^dMhMGj; pV^W;LB&*!U@dG4)8g!Ej8@Vi!vZOBjwL>oGSA`##? z=A;%&c3XDO_>^R|;f;x#a%E6jW_EXxCXA1bY3u^$6+amfnU%Sv$L9r0XqiN=ySv7Q zGn~zmRdlkw<@8O~DFK};eR#k}{oI<#ttt1K`E2j{HwJ;LsJN=K+s-aGc6H yFGu|PO& literal 39193 zcmeIbcT`i&*DsEU1sgp2C@2CpKoq1Yodhc)Rir6ZkY1&?fO+%*m8O7n2o^x;y%VA! zAYJJ-AcPK~2S`ZnM5QVG{`lSX-ur#;i)*>o;hZxwdv^KkJ$v??Jd>A|*uUr49tsMI z{gOAYDNs=C=%Aq3j@Z2u-1)}7Eu4aa>VVPJtMVpSC9YbST3Fq&ysxV#p=YjVWu&Vh zafO0H(A!%{!|=G$-;l7}^z&zG9yi3q_#Gf#3dk^3`1+)6kf@h#*1_R#nsNEXm(k7$ zBUZj%!aE~&$lsS~Kfdqw? yhX=heAatlJoR zKNYo&VNaBHHP2<^NCqQG@)FygcK+LvAXstgyB-Dc@aHs*M^u;b<pI8TI4@#J`_o+mr|0fb?XsR|_YuAHK9T)Y<@GQ8@u9M( z>kA*_-LkHCnv}FZmh(;LhW(A+o--G8&ML6SgBGFG*J1NC`GQgUXIhe8t`YU*X$p<} zf`6^{OhQU@m7EKX588Z)6ZnccAKMUOK#P9tVEa+DqICQkA(KZde2=ZGLioOZ8YYVF z@8>vEoN%X^Qe=$IJ+1UsG2}0sdODxee?*@WHqL2aB-~sl6s`aKGRL@Q@6=z)v5*p# z50ODvmb5st|1M<_bJq;tXH(8?OPF1>Q;7)K{oPN9ZtpQK7B1z?$t&e {h9s8)FCB@gT23qe|+LNdYD`L(6tvQXpCIcwBN?=IW}nc zdNA7X^?YOZDEbro0+AwH(U@31GWeDIN(|}rL|xn`tfA8Bf_<)zhViy_)cTgE8+atS zB>f)GA3x>a$9?qV+Y^3lU+z5cdhk{F)`3f9fv2n=Oy_=WaNKrU0U_X7;5lf2xs+(| zAh bX<;QW|b-FmzFAU!qK2rI-x^rLc(?h960R$334tq~H&4UX80CgaI( zQ}^^61>NpIT)Ef6DSGz7zha+?KRu9=8$9g(g!(?0lbt~PA#Xa3{GDYf77yZK36JPL zSdW?~q&M|I6>kLYWOIq>#qw>JHSc|LS`vEn)W2r=PEW>WW9(d|#M+YKrPh b{NQ{+UEwv$*K_I}oF zsssQbsHZBaFC#;74qWf1*ydqGK?Sb1f&a(Ae+ml9D9YtKJ36AaFJEs*kRM(QdhJU= zAx aP!Hx*Y`>taG20XbInssh8pxho z?5L=ig^UoIEMTtjuV%-pAF?BN)1DB2M6r#MYU^T)3I9hl0X4M4y}rmc&^{8aSV&OC zWDB6tZlgFP%Vwj?o|s?@yVR>& reELb8wT=?2 _su0lY-?L(Qp znjpGMzN=$~-O}Q|s?(YEjUuRbcZe|L^d3ARt{sL-k?iz6ke%4lX~K*?kC}I{$G3<4 zK8cBo+u->6Ck3*PDBwO+HH@#ODjLIg5hvTD3rM$Xx@@lhJ}-w8-CQ!ZG@Dqh_@hTC zF^b#No0# zP aZ=yXhVFNFW Wh*KpijA50;c{x!e|o@>!g92(c6m{bGyP{oo$ z-BlTsT&=WZYU@U@W6#D$Cc4}I;f0(4;W*dr;Rp(N>yB~epoh~hdG%Td;rRRCNR6<% zFs|QiivpfbJ%2MKeZf&+oflmtX8xLZ_xP#Rzt>R-rvMI*;?2Q7ssC{Yv^=oYa*IPJ z#J4noV#`*xFySZ2Z)w7RrY682KRR~nm%sAf?bTa+16Mjl*PP_*7MwdqV-W;(AXWG1 z_cP7jSik0E0$vxI)Z_faxz>52 mJh ziSBRzDjdK~N8lW$56xp?FOs`3$AT4HQ93~ETJC)Y}Yen G}yqt>;X8BIPNH3#neg&5@?m$O?rO1<*j^oIF4O z+2ww~8HVGbfz-4PX{g+{3UddE&LY^o?FdosABmQA!<;{KlU^oKy5-^W%-8o9wHL&U zbh>soH|Arons_%PJ7sZwqvZkW{(Xch^yCpavtBgI6RqL|$s9FR*v4HCUJ>u^`W(#Y z_a{@xclwXe=G1Z~KsN5$=Hv-WtWMPd$z@X3^VWLr?v8;Y%7~4-&I2n1dmep{-*WN4 z<_OOPh+w@pUT)*AHGtCR)2CnEM4TPKWODy)WLr }b6Y^P bCt axk>ryv6TOFwL7uyq%xHjrd6L8MG?oc%(b0_(`@!Dd(f zOcIdqn1B>#^46r+F=tECe`FT^Yt 9wT z7N{eeAE_^F%wG!$FQmk~BTbm{mg{H3?{*U>iw_pdU$8*p`$P8DsA?at6yXb030<2j z@&rjwI!xA4c ysU#3e?g@SlzEth83u6CanA4M!I0$y-p-Ym4 z&-|Hxt+j4_a$0d;f+#$8H=(9Kmd0&?-K0SJhbGCgSf*%Ic-A6@6Vr=^602)7_jVW< zWc^AD{jmGJAg@TbGn*04My)QWQgQneiJ!>=Lckd?Lq6!%uASsEOw7L#|MV2=$+;*7 zpWUmDRIPZd)r#~HKy2f%tmQwvv)0O!vXz7MHZxKT89CO` dg4E@$2Cbc57+25NGC$VWGKBD%Ru zC9tc=xX8EC8B5TfzuV}{!MZGr6XH=)Al)6LYj60ALn!kd=4e z=}S|gWv)zUoKw &ivj_Uf}| zk5?6Xu1%yTNCNK@V OZ zolRSG%{ yuA0qGWy{)&OO+{Z@l1NE`Q+IKZ_v*ENH7%%R{j=jRF7ebJc zM|DLn6p@-3ilGn2Z?D-ka{?s}IBq?Il@f>dN9dUgdpqPSawUqYah)+m#V#5{Pgh11 zcLOmZQGdK-C7E%UvR3mfX6o~E0gjxR(rSfqX0DvqiEDtX1(e93D|)XlD=@rkhYvqg z#C6g|_TT%{Z^zai7%o2zHp{B{_D1}Y^~xweFT{(%= fh-vanyE&s#jzcuF#0&DR{Tfp!i1q_xiVh^wNdGe<%N%K9! z56>oRqVo`JY;_kyTN3Op__fWgrh#bBa{JJD?!UU=q}_V= z79(Q-wnjnq_NvEx25JkRxxYTVdMEa;9EYTU;^O)O>wcUmkp^xL4R)l#tVR=O4$|hZ z9dlP+_3Gg_x~D!XM0OG?=4=ZV${FC!Ssw;Hmcs;X*==q{!{K?xn^Y(Y05@)aV=U@7 z_x`IncPvOW$!cV-j=we J@r@1b9M%Di}$QXZ!=ij(%8xWgd^&1Y{THk1y zqt+m6J TsOQ#$9)Uq+kBg%*UG)qJ{%m2-!An zoD+nb17+#-YP-C#EP$Nh_$@K7(X}lxuj=JbL-;SXX4GW?HrusGaO5z3nMC{kXu`L> z&1Lz%k$w2_eq2RS5U*&HEt1K0r2_B;s8^44hUBf5{nsU+1}Jflr+htH<)&gW>I;{1 z{R?GQ4}tU@o_NXYXz-SAV{kCMiIU3^xX6*XYG#H?-LHbAgoG=amm;b?ItW;`g~e|x zn@r+i!0`A7VoX;1lc63pBX+heJt4V{Jx02$ZuvevI8j!4R?>DEE!HR;fW{~jfQZe# z_71B(t$pDs0B;G&xT_nZg9g~C%4aGn-UO)wz*k>xDKKu_^WZ%Ih_myIk-y9Qr*C}& zu;aGh@LKR*r9r$IOh2g%JF%A3WgG_IbYMAQ>FcMs;!#!`lE=yen{)PJt{$6I&z}YS z2zV#FanE&`zy*v V_{W!zerVrGP=3$) js zma2Ui(4c=<*Jk%dgYTwo0h$PFK9sV_7HtAJ_kw)RBTTiiI*kJ3Lyl0lY;tx*8614k z`b^Kb$sVp9P)77D)#xV57!T4y$vR#KF0Z${W&Ftp_J?(%oZH0ylt9Fx!gIr8)tuJ` z+i)`E*o$q$ZzN){6LJo~Mm+u7tLsf~nF%Pe3Fb#WTH4qI0naIMN?u=2E|a=mjM`+> z(U3zC*WX^AWy&59&zvhAq}}8(09t$!o3`7=@>B(