Skip to content

Commit

Permalink
Merge branch 'develop' into release/v2.36.0
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicusor Serban committed Dec 3, 2024
2 parents 13b0c5c + f4fb7cb commit 43f3fea
Show file tree
Hide file tree
Showing 6 changed files with 334 additions and 12 deletions.
10 changes: 9 additions & 1 deletion src/stan/io/stan_csv_reader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,6 +250,10 @@ class stan_csv_reader {

int rows = lines - 3;
int cols = std::count(line.begin(), line.end(), ',') + 1;
if (cols == 1) {
// model has no parameters
return;
}
adaptation.metric.resize(rows, cols);
char comment; // Buffer for comment indicator, #

Expand Down Expand Up @@ -330,7 +334,11 @@ class stan_csv_reader {
for (int col = 0; col < cols; col++) {
std::getline(ls, line, ',');
boost::trim(line);
std::stringstream(line) >> samples(row, col);
try {
samples(row, col) = static_cast<double>(std::stold(line));
} catch (const std::out_of_range& e) {
samples(row, col) = std::numeric_limits<double>::quiet_NaN();
}
}
}
}
Expand Down
42 changes: 32 additions & 10 deletions src/stan/mcmc/chainset.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,7 @@ class chainset {
double median(const std::string& name) const { return median(index(name)); }

/**
* Compute maximum absolute deviation (mad) for specified parameter.
* Compute median absolute deviation (mad) for specified parameter.
*
* Follows R implementation: constant * median(abs(x - center))
* where the value of center is median(x) and the constant is 1.4826,
Expand All @@ -263,16 +263,20 @@ class chainset {
* @param index parameter index
* @return sample mad
*/
double max_abs_deviation(const int index) const {
double med_abs_deviation(const int index) const {
Eigen::MatrixXd draws = samples(index);
auto center = median(index);
Eigen::MatrixXd abs_dev = (draws.array() - center).abs();
Eigen::Map<Eigen::VectorXd> map(abs_dev.data(), abs_dev.size());
return 1.4826 * stan::math::quantile(map, 0.5);
try {
return 1.4826 * stan::math::quantile(map, 0.5);
} catch (const std::logic_error& e) {
return std::numeric_limits<double>::quiet_NaN();
}
}

/**
* Compute maximum absolute deviation (mad) for specified parameter.
* Compute median absolute deviation (mad) for specified parameter.
*
* Follows R implementation: constant * median(abs(x - center))
* where the value of center is median(x) and the constant is 1.4826,
Expand All @@ -282,27 +286,33 @@ class chainset {
* @param name parameter name
* @return sample mad
*/
double max_abs_deviation(const std::string& name) const {
return max_abs_deviation(index(name));
double med_abs_deviation(const std::string& name) const {
return med_abs_deviation(index(name));
}

/**
* Compute the quantile value of the specified parameter
* at the specified probability.
* at the specified probability via a call to stan::math::quantile.
*
* Calls stan::math::quantile which throws
* std::invalid_argument If any element of samples_vec is NaN, or size 0.
* and std::domain_error If `p<0` or `p>1`.
* If this happens, error will be caught, quantile value is NaN.
*
* @param index parameter index
* @param prob probability
* @return parameter value at quantile
*/
double quantile(const int index, const double prob) const {
// Ensure the probability is within [0, 1]
Eigen::MatrixXd draws = samples(index);
Eigen::Map<Eigen::VectorXd> map(draws.data(), draws.size());
return stan::math::quantile(map, prob);
double result;
try {
result = stan::math::quantile(map, prob);
} catch (const std::logic_error& e) {
return std::numeric_limits<double>::quiet_NaN();
}
return result;
}

/**
Expand All @@ -321,6 +331,11 @@ class chainset {
* Compute the quantile values of the specified parameter
* for a set of specified probabilities.
*
* Calls stan::math::quantile which throws
* std::invalid_argument If any element of samples_vec is NaN, or size 0.
* and std::domain_error If `p<0` or `p>1`.
* If this happens, error will be caught, quantile value is NaN.
*
* @param index parameter index
* @param probs vector of probabilities
* @return vector of parameter values for quantiles
Expand All @@ -332,7 +347,14 @@ class chainset {
Eigen::MatrixXd draws = samples(index);
Eigen::Map<Eigen::VectorXd> map(draws.data(), draws.size());
std::vector<double> probs_vec(probs.data(), probs.data() + probs.size());
std::vector<double> quantiles = stan::math::quantile(map, probs_vec);
std::vector<double> quantiles;
try {
quantiles = stan::math::quantile(map, probs_vec);
} catch (const std::logic_error& e) {
Eigen::VectorXd nans(probs.size());
nans.setConstant(std::numeric_limits<double>::quiet_NaN());
return nans;
}
return Eigen::Map<Eigen::VectorXd>(quantiles.data(), quantiles.size());
}

Expand Down
24 changes: 24 additions & 0 deletions src/test/unit/io/stan_csv_reader_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,8 @@ class StanIoStanCsvReader : public testing::Test {
"src/test/unit/io/test_csv_files/bernoulli_warmup.csv");
fixed_param_stream.open(
"src/test/unit/io/test_csv_files/fixed_param_output.csv");
no_params_stream.open(
"src/test/unit/io/test_csv_files/no_parameters_hmc.csv");
}

void TearDown() {
Expand All @@ -46,6 +48,7 @@ class StanIoStanCsvReader : public testing::Test {
bernoulli_thin_stream.close();
bernoulli_warmup_stream.close();
fixed_param_stream.close();
no_params_stream.close();
}

std::ifstream blocker0_stream, epil0_stream;
Expand All @@ -58,6 +61,7 @@ class StanIoStanCsvReader : public testing::Test {
std::ifstream bernoulli_thin_stream;
std::ifstream bernoulli_warmup_stream;
std::ifstream fixed_param_stream;
std::ifstream no_params_stream;
};

TEST_F(StanIoStanCsvReader, read_metadata1) {
Expand Down Expand Up @@ -570,6 +574,14 @@ TEST_F(StanIoStanCsvReader, fixed_param) {
ASSERT_EQ(10, fixed_param.samples.rows());
}

TEST_F(StanIoStanCsvReader, no_parameters) {
stan::io::stan_csv no_parameters_hmc;
std::stringstream out;
no_parameters_hmc = stan::io::stan_csv_reader::parse(no_params_stream, &out);
ASSERT_EQ(100, no_parameters_hmc.samples.rows());
ASSERT_EQ(0, no_parameters_hmc.adaptation.metric.size());
}

TEST_F(StanIoStanCsvReader, no_samples) {
std::ifstream no_samples_stream;
no_samples_stream.open(
Expand All @@ -590,4 +602,16 @@ TEST_F(StanIoStanCsvReader, variational) {
= stan::io::stan_csv_reader::parse(variational_stream, &out);
variational_stream.close();
ASSERT_EQ(1000, variational.metadata.num_samples);
ASSERT_EQ(0, variational.adaptation.metric.size());
}

TEST_F(StanIoStanCsvReader, read_nans) {
std::ifstream datagen_stream;
datagen_stream.open("src/test/unit/mcmc/test_csv_files/datagen_output.csv",
std::ifstream::in);
std::stringstream out;
stan::io::stan_csv datagen
= stan::io::stan_csv_reader::parse(datagen_stream, &out);
datagen_stream.close();
ASSERT_TRUE(std::isnan(datagen.samples(0, 2)));
}
157 changes: 157 additions & 0 deletions src/test/unit/io/test_csv_files/no_parameters_hmc.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
# stan_version_major = 2
# stan_version_minor = 35
# stan_version_patch = 0
# model = sim_model
# start_datetime = 2024-11-25 16:05:51 UTC
# method = sample (Default)
# sample
# num_samples = 100
# num_warmup = 100
# save_warmup = false (Default)
# thin = 1 (Default)
# adapt
# engaged = true (Default)
# gamma = 0.05 (Default)
# delta = 0.8 (Default)
# kappa = 0.75 (Default)
# t0 = 10 (Default)
# init_buffer = 75 (Default)
# term_buffer = 50 (Default)
# window = 25 (Default)
# save_metric = false (Default)
# algorithm = hmc (Default)
# hmc
# engine = nuts (Default)
# nuts
# max_depth = 10 (Default)
# metric = diag_e (Default)
# metric_file = (Default)
# stepsize = 1 (Default)
# stepsize_jitter = 0 (Default)
# num_chains = 1 (Default)
# id = 1 (Default)
# data
# file = (Default)
# init = 2 (Default)
# random
# seed = 1799096017 (Default)
# output
# file = output.csv (Default)
# diagnostic_file = (Default)
# refresh = 100 (Default)
# sig_figs = -1 (Default)
# profile_file = profile.csv (Default)
# save_cmdstan_config = false (Default)
# num_threads = 1 (Default)
# stanc_version = %%NAME%%3 %%VERSION%%
# stancflags =
lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,x
# Adaptation terminated
# Step size = 1.35996e+14
# Diagonal elements of inverse mass matrix:
#
0,1,1.35996e+14,1,1,0,0,8.57152
0,1,1.35996e+14,1,1,0,0,-4.56101
0,1,1.35996e+14,1,1,0,0,9.30019
0,1,1.35996e+14,1,1,0,0,-14.6254
0,1,1.35996e+14,1,1,0,0,8.30123
0,1,1.35996e+14,1,1,0,0,19.2067
0,1,1.35996e+14,1,1,0,0,12.6027
0,1,1.35996e+14,1,1,0,0,-5.70173
0,1,1.35996e+14,1,1,0,0,3.04002
0,1,1.35996e+14,1,1,0,0,-8.28215
0,1,1.35996e+14,1,1,0,0,3.20386
0,1,1.35996e+14,1,1,0,0,-6.45707
0,1,1.35996e+14,1,1,0,0,-5.29269
0,1,1.35996e+14,1,1,0,0,5.1527
0,1,1.35996e+14,1,1,0,0,-3.69775
0,1,1.35996e+14,1,1,0,0,-1.77925
0,1,1.35996e+14,1,1,0,0,-13.7562
0,1,1.35996e+14,1,1,0,0,-0.789554
0,1,1.35996e+14,1,1,0,0,-14.0911
0,1,1.35996e+14,1,1,0,0,-4.71911
0,1,1.35996e+14,1,1,0,0,-12.7874
0,1,1.35996e+14,1,1,0,0,15.029
0,1,1.35996e+14,1,1,0,0,-6.30709
0,1,1.35996e+14,1,1,0,0,3.50805
0,1,1.35996e+14,1,1,0,0,9.01495
0,1,1.35996e+14,1,1,0,0,5.9809
0,1,1.35996e+14,1,1,0,0,11.02
0,1,1.35996e+14,1,1,0,0,9.68621
0,1,1.35996e+14,1,1,0,0,-4.20659
0,1,1.35996e+14,1,1,0,0,-7.71128
0,1,1.35996e+14,1,1,0,0,-13.6125
0,1,1.35996e+14,1,1,0,0,0.541567
0,1,1.35996e+14,1,1,0,0,-6.96969
0,1,1.35996e+14,1,1,0,0,2.49902
0,1,1.35996e+14,1,1,0,0,8.31074
0,1,1.35996e+14,1,1,0,0,13.4856
0,1,1.35996e+14,1,1,0,0,-4.60591
0,1,1.35996e+14,1,1,0,0,6.84733
0,1,1.35996e+14,1,1,0,0,-16.8676
0,1,1.35996e+14,1,1,0,0,4.43181
0,1,1.35996e+14,1,1,0,0,-10.6785
0,1,1.35996e+14,1,1,0,0,-5.56113
0,1,1.35996e+14,1,1,0,0,-0.95401
0,1,1.35996e+14,1,1,0,0,11.2198
0,1,1.35996e+14,1,1,0,0,3.43417
0,1,1.35996e+14,1,1,0,0,-11.2942
0,1,1.35996e+14,1,1,0,0,-14.3029
0,1,1.35996e+14,1,1,0,0,3.69492
0,1,1.35996e+14,1,1,0,0,0.319324
0,1,1.35996e+14,1,1,0,0,-5.95097
0,1,1.35996e+14,1,1,0,0,5.99333
0,1,1.35996e+14,1,1,0,0,-6.59629
0,1,1.35996e+14,1,1,0,0,14.1795
0,1,1.35996e+14,1,1,0,0,-7.58818
0,1,1.35996e+14,1,1,0,0,4.89377
0,1,1.35996e+14,1,1,0,0,4.63195
0,1,1.35996e+14,1,1,0,0,-4.62905
0,1,1.35996e+14,1,1,0,0,-11.4145
0,1,1.35996e+14,1,1,0,0,4.03017
0,1,1.35996e+14,1,1,0,0,-10.0459
0,1,1.35996e+14,1,1,0,0,-11.8674
0,1,1.35996e+14,1,1,0,0,-0.161997
0,1,1.35996e+14,1,1,0,0,-5.75037
0,1,1.35996e+14,1,1,0,0,-13.3027
0,1,1.35996e+14,1,1,0,0,4.86817
0,1,1.35996e+14,1,1,0,0,11.1937
0,1,1.35996e+14,1,1,0,0,13.918
0,1,1.35996e+14,1,1,0,0,-12.2423
0,1,1.35996e+14,1,1,0,0,22.3588
0,1,1.35996e+14,1,1,0,0,-8.30628
0,1,1.35996e+14,1,1,0,0,-6.87127
0,1,1.35996e+14,1,1,0,0,12.721
0,1,1.35996e+14,1,1,0,0,-7.86135
0,1,1.35996e+14,1,1,0,0,-8.56196
0,1,1.35996e+14,1,1,0,0,-4.04709
0,1,1.35996e+14,1,1,0,0,-21.6766
0,1,1.35996e+14,1,1,0,0,-19.6485
0,1,1.35996e+14,1,1,0,0,1.99421
0,1,1.35996e+14,1,1,0,0,11.2645
0,1,1.35996e+14,1,1,0,0,-9.35154
0,1,1.35996e+14,1,1,0,0,-3.37081
0,1,1.35996e+14,1,1,0,0,2.46874
0,1,1.35996e+14,1,1,0,0,7.28248
0,1,1.35996e+14,1,1,0,0,19.3846
0,1,1.35996e+14,1,1,0,0,-4.92502
0,1,1.35996e+14,1,1,0,0,-11.4687
0,1,1.35996e+14,1,1,0,0,3.12818
0,1,1.35996e+14,1,1,0,0,-5.79636
0,1,1.35996e+14,1,1,0,0,12.2196
0,1,1.35996e+14,1,1,0,0,5.59427
0,1,1.35996e+14,1,1,0,0,-22.9084
0,1,1.35996e+14,1,1,0,0,2.70951
0,1,1.35996e+14,1,1,0,0,-0.604509
0,1,1.35996e+14,1,1,0,0,1.06201
0,1,1.35996e+14,1,1,0,0,-3.10986
0,1,1.35996e+14,1,1,0,0,6.57804
0,1,1.35996e+14,1,1,0,0,9.0426
0,1,1.35996e+14,1,1,0,0,3.30854
0,1,1.35996e+14,1,1,0,0,8.85378
0,1,1.35996e+14,1,1,0,0,-1.71584
#
# Elapsed Time: 0 seconds (Warm-up)
# 0 seconds (Sampling)
# 0 seconds (Total)
#
31 changes: 30 additions & 1 deletion src/test/unit/mcmc/chainset_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ TEST_F(McmcChains, summary_stats) {
EXPECT_NEAR(theta_sd_expect, bern_chains.sd("theta"), 1e-5);

double theta_mad_expect = 0.12309;
EXPECT_NEAR(theta_mad_expect, bern_chains.max_abs_deviation("theta"), 1e-5);
EXPECT_NEAR(theta_mad_expect, bern_chains.med_abs_deviation("theta"), 1e-5);

double theta_mcse_mean_expect = 0.003234;
EXPECT_NEAR(theta_mcse_mean_expect, bern_chains.mcse_mean("theta"), 1e-4);
Expand Down Expand Up @@ -207,3 +207,32 @@ TEST_F(McmcChains, summary_stats) {
EXPECT_NEAR(theta_ac(i), theta_ac_expect(i), 0.0005);
}
}

TEST_F(McmcChains, quantile_tests) {
std::ifstream datagen_stream;
datagen_stream.open("src/test/unit/mcmc/test_csv_files/datagen_output.csv",
std::ifstream::in);
stan::io::stan_csv datagen_csv
= stan::io::stan_csv_reader::parse(datagen_stream, &out);
datagen_stream.close();
stan::mcmc::chainset datagen_chains(datagen_csv);

Eigen::VectorXd probs(6);
probs << 0.0, 0.01, 0.05, 0.95, 0.99, 1.0;
Eigen::VectorXd stepsize_quantiles
= datagen_chains.quantiles("stepsize__", probs);
for (size_t i = 0; i < probs.size(); ++i) {
EXPECT_TRUE(std::isnan(stepsize_quantiles(i)));
}

double stepsize_MAD = datagen_chains.med_abs_deviation("stepsize__");
EXPECT_TRUE(std::isnan(stepsize_MAD));

Eigen::VectorXd bad_probs(3);
bad_probs << 5, 50, 95;
Eigen::VectorXd y_sim_quantiles
= datagen_chains.quantiles("y_sim[1]", bad_probs);
for (size_t i = 0; i < bad_probs.size(); ++i) {
EXPECT_TRUE(std::isnan(y_sim_quantiles(i)));
}
}
Loading

0 comments on commit 43f3fea

Please sign in to comment.