diff --git a/src/cmdstan/diagnose.cpp b/src/cmdstan/diagnose.cpp index 76d8fdbbf8..9a451e9d57 100644 --- a/src/cmdstan/diagnose.cpp +++ b/src/cmdstan/diagnose.cpp @@ -1,12 +1,15 @@ +#include #include -#include +#include #include #include #include #include #include -double RHAT_MAX = 1.05; +using cmdstan::return_codes; + +double RHAT_MAX = 1.01499; // round to 1.01 void diagnose_usage() { std::cout << "USAGE: diagnose [ ... ]" @@ -26,7 +29,7 @@ void diagnose_usage() { int main(int argc, const char *argv[]) { if (argc == 1) { diagnose_usage(); - return 0; + return return_codes::OK; } // Parse any arguments specifying filenames @@ -45,49 +48,47 @@ int main(int argc, const char *argv[]) { if (!filenames.size()) { std::cout << "No valid input files, exiting." << std::endl; - return 0; + return return_codes::NOT_OK; } std::cout << std::fixed << std::setprecision(2); - // Parse specified files - std::cout << "Processing csv files: " << filenames[0]; - ifstream.open(filenames[0].c_str()); - - stan::io::stan_csv stan_csv - = stan::io::stan_csv_reader::parse(ifstream, &std::cout); - stan::mcmc::chains<> chains(stan_csv); - ifstream.close(); - - if (filenames.size() > 1) - std::cout << ", "; - else - std::cout << std::endl << std::endl; - - for (std::vector::size_type chain = 1; chain < filenames.size(); - ++chain) { - std::cout << filenames[chain]; - ifstream.open(filenames[chain].c_str()); - stan_csv = stan::io::stan_csv_reader::parse(ifstream, &std::cout); - chains.add(stan_csv); - ifstream.close(); - if (chain < filenames.size() - 1) - std::cout << ", "; - else - std::cout << std::endl << std::endl; + std::vector csv_parsed; + for (int i = 0; i < filenames.size(); ++i) { + std::ifstream infile; + std::stringstream out; + stan::io::stan_csv sample; + infile.open(filenames[i].c_str()); + try { + sample = stan::io::stan_csv_reader::parse(infile, &out); + // csv_reader warnings are errors - fail fast. + if (!out.str().empty()) { + throw std::invalid_argument(out.str()); + } + csv_parsed.push_back(sample); + } catch (const std::invalid_argument &e) { + std::cout << "Cannot parse input csv file: " << filenames[i] << e.what() + << "." << std::endl; + return return_codes::NOT_OK; + } } - + stan::mcmc::chainset chains(csv_parsed); + stan::io::stan_csv_metadata metadata = csv_parsed[0].metadata; + std::vector param_names = csv_parsed[0].header; + size_t num_params = param_names.size(); int num_samples = chains.num_samples(); std::vector bad_n_eff_names; std::vector bad_rhat_names; bool has_errors = false; - for (int i = 0; i < chains.num_params(); ++i) { - if (chains.param_name(i) == std::string("treedepth__")) { + for (int i = 0; i < num_params; ++i) { + if (param_names[i] == std::string("treedepth__")) { std::cout << "Checking sampler transitions treedepth." << std::endl; - int max_limit = stan_csv.metadata.max_depth; + int max_limit = metadata.max_depth; long n_max = 0; - Eigen::VectorXd t_samples = chains.samples(i); + Eigen::MatrixXd draws = chains.samples(i); + Eigen::VectorXd t_samples + = Eigen::Map(draws.data(), draws.size()); for (long n = 0; n < t_samples.size(); ++n) { if (t_samples(n) >= max_limit) { ++n_max; @@ -109,7 +110,7 @@ int main(int argc, const char *argv[]) { std::cout << "Treedepth satisfactory for all transitions." << std::endl << std::endl; } - } else if (chains.param_name(i) == std::string("divergent__")) { + } else if (param_names[i] == std::string("divergent__")) { std::cout << "Checking sampler transitions for divergences." << std::endl; int n_divergent = chains.samples(i).sum(); if (n_divergent > 0) { @@ -129,26 +130,22 @@ int main(int argc, const char *argv[]) { std::cout << "No divergent transitions found." << std::endl << std::endl; } - } else if (chains.param_name(i) == std::string("energy__")) { + } else if (param_names[i] == std::string("energy__")) { std::cout << "Checking E-BFMI - sampler transitions HMC potential energy." << std::endl; - Eigen::VectorXd e_samples = chains.samples(i); + Eigen::MatrixXd draws = chains.samples(i); + Eigen::VectorXd e_samples + = Eigen::Map(draws.data(), draws.size()); double delta_e_sq_mean = 0; - double e_mean = 0; - double e_var = 0; - e_mean += e_samples(0); - e_var += e_samples(0) * (e_samples(0) - e_mean); + double e_mean = chains.mean(i); + double e_var = chains.variance(i); for (long n = 1; n < e_samples.size(); ++n) { double e = e_samples(n); double delta_e_sq = (e - e_samples(n - 1)) * (e - e_samples(n - 1)); double d = delta_e_sq - delta_e_sq_mean; delta_e_sq_mean += d / n; d = e - e_mean; - e_mean += d / (n + 1); - e_var += d * (e - e_mean); } - - e_var /= static_cast(e_samples.size() - 1); double e_bfmi = delta_e_sq_mean / e_var; double e_bfmi_threshold = 0.3; if (e_bfmi < e_bfmi_threshold) { @@ -163,14 +160,16 @@ int main(int argc, const char *argv[]) { } else { std::cout << "E-BFMI satisfactory." << std::endl << std::endl; } - } else if (chains.param_name(i).find("__") == std::string::npos) { - double n_eff = chains.effective_sample_size(i); + } else if (param_names[i].find("__") == std::string::npos) { + auto [ess_bulk, ess_tail] = chains.split_rank_normalized_ess(i); + double n_eff = ess_bulk < ess_tail ? ess_bulk : ess_tail; if (n_eff / num_samples < 0.001) - bad_n_eff_names.push_back(chains.param_name(i)); + bad_n_eff_names.push_back(param_names[i]); - double split_rhat = chains.split_potential_scale_reduction(i); + auto [rhat_bulk, rhat_tail] = chains.split_rank_normalized_rhat(i); + double split_rhat = rhat_bulk > rhat_tail ? rhat_bulk : rhat_tail; if (split_rhat > RHAT_MAX) - bad_rhat_names.push_back(chains.param_name(i)); + bad_rhat_names.push_back(param_names[i]); } } if (bad_n_eff_names.size() > 0) { @@ -187,13 +186,15 @@ int main(int argc, const char *argv[]) { << " may be substantially lower than quoted." << std::endl << std::endl; } else { - std::cout << "Effective sample size satisfactory." << std::endl + std::cout << "Rank-normalized split effective sample size satisfactory " + << "for all parameters." << std::endl << std::endl; } if (bad_rhat_names.size() > 0) { has_errors = true; - std::cout << "The following parameters had split R-hat greater than " + std::cout << "The following parameters had rank-normalized split R-hat " + "greater than " << RHAT_MAX << ":" << std::endl; std::cout << " "; for (size_t n = 0; n < bad_rhat_names.size() - 1; ++n) @@ -207,7 +208,8 @@ int main(int argc, const char *argv[]) { << " effective parameterization." << std::endl << std::endl; } else { - std::cout << "Split R-hat values satisfactory all parameters." << std::endl + std::cout << "Rank-normalized split R-hat values satisfactory " + << "for all parameters." << std::endl << std::endl; } if (!has_errors) @@ -215,5 +217,5 @@ int main(int argc, const char *argv[]) { else std::cout << "Processing complete." << std::endl; - return 0; + return return_codes::OK; } diff --git a/src/cmdstan/stansummary.cpp b/src/cmdstan/stansummary.cpp index fe4483d4b3..518a2f75bf 100644 --- a/src/cmdstan/stansummary.cpp +++ b/src/cmdstan/stansummary.cpp @@ -1,6 +1,5 @@ #include #include -#include #include #include #include @@ -34,7 +33,7 @@ Example: stansummary model_chain_1.csv model_chain_2.csv -c, --csv_filename [file] Write statistics to a csv file. -h, --help Produce help message, then exit. -p, --percentiles [values] Percentiles to report as ordered set of - comma-separated numbers from (0.1,99.9), inclusive. + comma-separated numbers from (0.0,100.0), inclusive. Default is 5,50,95. -s, --sig_figs [n] Significant figures reported. Default is 2. Must be an integer from (1, 18), inclusive. @@ -140,8 +139,8 @@ Example: stansummary model_chain_1.csv model_chain_2.csv // check for stan csv file parse errors written to output stream std::stringstream cout_ss; - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); // Get column headers for sampler, model params size_t max_name_length = 0; diff --git a/src/cmdstan/stansummary_helper.hpp b/src/cmdstan/stansummary_helper.hpp index 4fb482852e..6ac7496325 100644 --- a/src/cmdstan/stansummary_helper.hpp +++ b/src/cmdstan/stansummary_helper.hpp @@ -1,7 +1,7 @@ #ifndef CMDSTAN_STANSUMMARY_HELPER_HPP #define CMDSTAN_STANSUMMARY_HELPER_HPP -#include +#include #include #include #include @@ -163,10 +163,12 @@ bool is_container(const std::string ¶meter_name) { /** * Return parameter name corresponding to column label. * + * @tparam Tchains - either a mcmc::chains or mcmc::chainset object * @param in column index * @return variable name */ -std::string base_param_name(const stan::mcmc::chains<> &chains, int index) { +template +std::string base_param_name(const Tchains &chains, int index) { std::string name = chains.param_name(index); return name.substr(0, name.find("[")); } @@ -174,11 +176,13 @@ std::string base_param_name(const stan::mcmc::chains<> &chains, int index) { /** * Return parameter name corresponding to column label. * + * @tparam Tchains - either a mcmc::chains or mcmc::chainset object * @param in set of samples from one or more chains * @param in column index * @return parameter name */ -std::string matrix_index(const stan::mcmc::chains<> &chains, int index) { +template +std::string matrix_index(const Tchains &chains, int index) { std::string name = chains.param_name(index); return name.substr(name.find("[")); } @@ -186,12 +190,13 @@ std::string matrix_index(const stan::mcmc::chains<> &chains, int index) { /** * Return vector of dimensions for container variable. * + * @tparam Tchains - either a mcmc::chains or mcmc::chainset object * @param in set of samples from one or more chains * @param in column index of first container element * @return vector of dimensions */ -std::vector dimensions(const stan::mcmc::chains<> &chains, - int start_index) { +template +std::vector dimensions(const Tchains &chains, int start_index) { std::vector dims; int dim; @@ -284,10 +289,10 @@ int matrix_index(std::vector &index, const std::vector &dims) { } /** - * Convert percentiles - int values in range (1,99) + * Convert percentiles - string-encoded doubles in range (0,100) * to probabilities - double values in range (0, 1). * - *

Input values must be in strictly increasing order. + * Input values must be in strictly increasing order. * * @param vector of strings * @return vector of doubles @@ -300,12 +305,12 @@ Eigen::VectorXd percentiles_to_probs( for (size_t i = 0; i < percentiles.size(); ++i) { try { pct = std::stod(percentiles[i]); - if (!std::isfinite(pct) || pct < 0.1 || pct > 99.9 || pct < cur_pct) + if (!std::isfinite(pct) || pct < 0.0 || pct > 100.0 || pct < cur_pct) throw std::exception(); cur_pct = pct; } catch (const std::exception &e) { throw std::invalid_argument( - "values must be in range (0.1,99.9)" + "values must be in range (0, 100)" ", inclusive, and strictly increasing."); } probs[i] = pct / 100.0; @@ -324,75 +329,62 @@ Eigen::VectorXd percentiles_to_probs( * @param out output stream * @return stan::mcmc::chains object */ -stan::mcmc::chains<> parse_csv_files(const std::vector &filenames, +stan::mcmc::chainset parse_csv_files(const std::vector &filenames, stan::io::stan_csv_metadata &metadata, Eigen::VectorXd &warmup_times, Eigen::VectorXd &sampling_times, Eigen::VectorXi &thin, std::ostream *out) { - // instantiate stan::mcmc::chains object by parsing first file + std::vector csvs; + csvs.reserve(filenames.size()); std::ifstream ifstream; - ifstream.open(filenames[0].c_str()); - stan::io::stan_csv stan_csv = stan::io::stan_csv_reader::parse(ifstream, out); - ifstream.close(); - if (stan_csv.samples.rows() < 1) { - std::stringstream message_stream(""); - message_stream << "No sampling draws found in Stan CSV file: " - << filenames[0] << "."; - throw std::invalid_argument(message_stream.str()); - } - warmup_times(0) = stan_csv.timing.warmup; - sampling_times(0) = stan_csv.timing.sampling; - stan::mcmc::chains<> chains(stan_csv); - thin(0) = stan_csv.metadata.thin; - metadata = stan_csv.metadata; - - // parse rest of input files, add to chains - for (std::vector::size_type chain = 1; chain < filenames.size(); - chain++) { + + for (size_t chain = 0; chain < filenames.size(); chain++) { ifstream.open(filenames[chain].c_str()); - stan_csv = stan::io::stan_csv_reader::parse(ifstream, out); + csvs.push_back(stan::io::stan_csv_reader::parse(ifstream, out)); ifstream.close(); + auto &stan_csv = csvs[chain]; if (stan_csv.samples.rows() < 1) { std::stringstream message_stream(""); message_stream << "No sampling draws found in Stan CSV file: " << filenames[chain] << "."; throw std::invalid_argument(message_stream.str()); } - chains.add(stan_csv); thin(chain) = stan_csv.metadata.thin; warmup_times(chain) = stan_csv.timing.warmup; sampling_times(chain) = stan_csv.timing.sampling; } - return chains; + metadata = csvs[0].metadata; + return stan::mcmc::chainset(csvs); } /** * Assemble vector of output column labels as follows: - * Mean, MCSE, StdDev, specified quantile labels, N_eff, N_eff/S, R-hat + * Mean, MCSE, StdDev, MAD, ... percentiles ..., ESS_bulk, ESS_tail, R_hat * * @param in vector of percentile values as strings * @return vector column labels */ std::vector get_header( const std::vector &percentiles) { - // Mean, MCSE, StdDev, ... percentiles ..., N_eff, N_eff/s, R_hat - std::vector header(percentiles.size() + 6); + // Mean, MCSE, StdDev, MAD, ... percentiles ..., ESS_bulk, ESS_tail, R_hat + std::vector header(percentiles.size() + 7); header.at(0) = "Mean"; header.at(1) = "MCSE"; header.at(2) = "StdDev"; + header.at(3) = "MAD"; for (size_t i = 0; i < percentiles.size(); ++i) { - header[i + 3] = percentiles[i] + '%'; + header[i + 4] = percentiles[i] + '%'; } - size_t offset = 3 + percentiles.size(); - header.at(offset) = "N_Eff"; - header.at(offset + 1) = "N_Eff/s"; + size_t offset = 4 + percentiles.size(); + header.at(offset) = "ESS_bulk"; + header.at(offset + 1) = "ESS_tail"; header.at(offset + 2) = "R_hat"; return header; } /** * Compute statistics for span of output columns - * Mean, MCSE, StdDev, specified quantile, N_eff, N_eff/S, R-hat + * Mean, MCSE, StdDev, MAD, ... percentiles ..., ESS_bulk, ESS_tail, R_hat * * @param in set of samples from one or more chains * @param in vector of warmup times (required for N_eff/S) @@ -402,12 +394,11 @@ std::vector get_header( * @param in span length * @param in out matrix of model param statistics */ -void get_stats(const stan::mcmc::chains<> &chains, +void get_stats(const stan::mcmc::chainset &chains, const Eigen::VectorXd &sampling_times, const Eigen::VectorXd &probs, std::vector cols, Eigen::MatrixXd ¶ms) { params.setZero(); - double total_sampling_time = sampling_times.sum(); if (params.rows() != cols.size()) { throw std::domain_error("get_stats: size mismatch"); @@ -416,18 +407,22 @@ void get_stats(const stan::mcmc::chains<> &chains, // Model parameters int i = 0; for (int i_chains : cols) { - double sd = chains.sd(i_chains); - double n_eff = chains.effective_sample_size(i_chains); params(i, 0) = chains.mean(i_chains); - params(i, 1) = sd / sqrt(n_eff); - params(i, 2) = sd; + params(i, 1) = chains.mcse_mean(i_chains); + params(i, 2) = chains.sd(i_chains); + params(i, 3) = chains.max_abs_deviation(i_chains); Eigen::VectorXd quantiles = chains.quantiles(i_chains, probs); for (int j = 0; j < quantiles.size(); j++) - params(i, 3 + j) = quantiles(j); - params(i, quantiles.size() + 3) = n_eff; - params(i, quantiles.size() + 4) = n_eff / total_sampling_time; - params(i, quantiles.size() + 5) - = chains.split_potential_scale_reduction(i_chains); + params(i, 4 + j) = quantiles(j); + + auto [ess_bulk, ess_tail] = chains.split_rank_normalized_ess(i_chains); + + params(i, quantiles.size() + 4) = ess_bulk; + params(i, quantiles.size() + 5) = ess_tail; + + auto [rhat_bulk, rhat_tail] = chains.split_rank_normalized_rhat(i_chains); + params(i, quantiles.size() + 6) + = rhat_bulk > rhat_tail ? rhat_bulk : rhat_tail; i++; } } @@ -473,7 +468,7 @@ void write_header(const std::vector &header, * @param in output format flag: true for csv; false for plain text * @param in output stream */ -void write_params(const stan::mcmc::chains<> &chains, +void write_params(const stan::mcmc::chainset &chains, const Eigen::MatrixXd ¶ms, const Eigen::VectorXi &col_widths, const Eigen::Matrix &chains, * @param in output format flag: true for csv; false for plain text * @param in output stream */ -void write_all_model_params(const stan::mcmc::chains<> &chains, +void write_all_model_params(const stan::mcmc::chainset &chains, const Eigen::MatrixXd ¶ms, const Eigen::VectorXi &col_widths, const Eigen::Matrix &chains, if (as_csv) { *out << "\"" << chains.param_name(row_maj_index_chains) << "\""; for (int j = 0; j < params.cols(); j++) { - *out << "," << std::fixed - << std::setprecision(compute_precision( - params(row_maj_index, j), sig_figs, false)) - << params(row_maj_index, j); + *out << "," << params(row_maj_index, j); } } else { *out << std::setw(max_name_length + 1) << std::left @@ -600,26 +592,17 @@ void write_all_model_params(const stan::mcmc::chains<> &chains, * @param in prefix string - used to output as comments in csv file * @param out output stream */ -void write_timing(const stan::mcmc::chains<> &chains, +void write_timing(const stan::mcmc::chainset &chains, const stan::io::stan_csv_metadata &metadata, const Eigen::VectorXd &warmup_times, const Eigen::VectorXd &sampling_times, const Eigen::VectorXi &thin, const std::string &prefix, std::ostream *out) { *out << prefix << "Inference for Stan model: " << metadata.model << std::endl - << prefix << chains.num_chains() << " chains: each with iter=(" - << chains.num_kept_samples(0); - for (int chain = 1; chain < chains.num_chains(); chain++) - *out << "," << chains.num_kept_samples(chain); - *out << ")"; - *out << "; warmup=(" << chains.warmup(0); - for (int chain = 1; chain < chains.num_chains(); chain++) - *out << "," << chains.warmup(chain); - *out << ")"; - *out << "; thin=(" << thin(0); - for (int chain = 1; chain < chains.num_chains(); chain++) - *out << "," << thin(chain); - *out << ")"; + << prefix << chains.num_chains() + << " chains: each with iter=" << metadata.num_samples; + *out << "; warmup=" << metadata.num_warmup; + *out << "; thin=" << metadata.thin; *out << "; " << chains.num_samples() << " iterations saved." << std::endl << prefix << std::endl; @@ -702,14 +685,16 @@ void write_sampler_info(const stan::io::stan_csv_metadata &metadata, *out << prefix << "Samples were drawn using " << metadata.algorithm << " with " << metadata.engine << "." << std::endl; *out << prefix - << "For each parameter, N_Eff is a crude measure of effective " - "sample size," + << "For each parameter, ESS_bulk and ESS_tail measure the " + "effective sample size " + << std::endl + << "for the entire sample (bulk) and for the " + "the .05 and .95 tails (tail), " << std::endl; *out << prefix - << "and R_hat is the potential scale reduction factor on split " - "chains (at " - << std::endl; - *out << prefix << "convergence, R_hat=1)." << std::endl; + << "and R_hat measures the potential scale reduction on split chains." + << std::endl + << "At convergence R_hat will be very close to 1.00." << std::endl; } /** @@ -725,11 +710,11 @@ void write_sampler_info(const stan::io::stan_csv_metadata &metadata, * @param in size of longest sampler param name */ // autocorrelation report prints to std::out -void autocorrelation(const stan::mcmc::chains<> &chains, +void autocorrelation(const stan::mcmc::chainset &chains, const stan::io::stan_csv_metadata &metadata, int autocorr_idx, int max_name_length) { int c = autocorr_idx - 1; - Eigen::MatrixXd autocorr(chains.num_params(), chains.num_samples(c)); + Eigen::MatrixXd autocorr(chains.num_params(), chains.num_samples()); for (int i = 0; i < chains.num_params(); ++i) { autocorr.row(i) = chains.autocorrelation(c, i); } diff --git a/src/test/interface/example_output/corr_gauss.nom b/src/test/interface/example_output/corr_gauss.nom index a6cd709481..64d07f8279 100644 --- a/src/test/interface/example_output/corr_gauss.nom +++ b/src/test/interface/example_output/corr_gauss.nom @@ -9,8 +9,8 @@ No divergent transitions found. Checking E-BFMI - sampler transitions HMC potential energy. E-BFMI satisfactory. -Effective sample size satisfactory. +Rank-normalized split effective sample size satisfactory for all parameters. -Split R-hat values satisfactory all parameters. +Rank-normalized split R-hat values satisfactory for all parameters. Processing complete. diff --git a/src/test/interface/example_output/corr_gauss_depth15.nom b/src/test/interface/example_output/corr_gauss_depth15.nom index f348c9e339..aa9095d6c7 100644 --- a/src/test/interface/example_output/corr_gauss_depth15.nom +++ b/src/test/interface/example_output/corr_gauss_depth15.nom @@ -7,8 +7,8 @@ No divergent transitions found. Checking E-BFMI - sampler transitions HMC potential energy. E-BFMI satisfactory. -Effective sample size satisfactory. +Rank-normalized split effective sample size satisfactory for all parameters. -Split R-hat values satisfactory all parameters. +Rank-normalized split R-hat values satisfactory for all parameters. Processing complete, no problems detected. diff --git a/src/test/interface/example_output/corr_gauss_depth8.nom b/src/test/interface/example_output/corr_gauss_depth8.nom index 32c345dc92..5a476ee328 100644 --- a/src/test/interface/example_output/corr_gauss_depth8.nom +++ b/src/test/interface/example_output/corr_gauss_depth8.nom @@ -9,8 +9,8 @@ No divergent transitions found. Checking E-BFMI - sampler transitions HMC potential energy. E-BFMI satisfactory. -Effective sample size satisfactory. +Rank-normalized split effective sample size satisfactory for all parameters. -Split R-hat values satisfactory all parameters. +Rank-normalized split R-hat values satisfactory for all parameters. Processing complete. diff --git a/src/test/interface/example_output/eight_schools.nom b/src/test/interface/example_output/eight_schools.nom index b306a5bb6b..ae4fc866dc 100644 --- a/src/test/interface/example_output/eight_schools.nom +++ b/src/test/interface/example_output/eight_schools.nom @@ -11,8 +11,8 @@ Checking E-BFMI - sampler transitions HMC potential energy. The E-BFMI, 0.26, is below the nominal threshold of 0.30 which suggests that HMC may have trouble exploring the target distribution. If possible, try to reparameterize the model. -Effective sample size satisfactory. +Rank-normalized split effective sample size satisfactory for all parameters. -Split R-hat values satisfactory all parameters. +Rank-normalized split R-hat values satisfactory for all parameters. Processing complete. diff --git a/src/test/interface/example_output/mix.nom b/src/test/interface/example_output/mix.nom index fbe76363e9..11693838f2 100644 --- a/src/test/interface/example_output/mix.nom +++ b/src/test/interface/example_output/mix.nom @@ -7,12 +7,10 @@ No divergent transitions found. Checking E-BFMI - sampler transitions HMC potential energy. E-BFMI satisfactory. -The following parameters had fewer than 0.001 effective draws per transition: - mu[1], mu[2], theta -Such low values indicate that the effective sample size estimators may be biased high and actual performance may be substantially lower than quoted. +Rank-normalized split effective sample size satisfactory for all parameters. -The following parameters had split R-hat greater than 1.05: - mu[1], mu[2], theta +The following parameters had rank-normalized split R-hat greater than 1.01: + mu[1], mu[2], sigma[1], theta Such high values indicate incomplete mixing and biased estimation. You should consider regularizating your model with additional prior information or a more effective parameterization. diff --git a/src/test/interface/multi_chain_test.cpp b/src/test/interface/multi_chain_test.cpp index 8564cf4e8e..40ba10d879 100644 --- a/src/test/interface/multi_chain_test.cpp +++ b/src/test/interface/multi_chain_test.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include @@ -30,8 +29,8 @@ TEST(interface, output_multi) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); constexpr std::array names{ "lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__", "divergent__", @@ -55,8 +54,8 @@ TEST(interface, output_multi) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); constexpr std::array names{ "lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__", "divergent__", diff --git a/src/test/interface/output_sig_figs_test.cpp b/src/test/interface/output_sig_figs_test.cpp index ea1f91b0f7..a96d22f0d8 100644 --- a/src/test/interface/output_sig_figs_test.cpp +++ b/src/test/interface/output_sig_figs_test.cpp @@ -1,7 +1,6 @@ #include #include #include -#include #include #include @@ -28,8 +27,8 @@ TEST(interface, output_sig_figs_1) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); EXPECT_NEAR(chains.samples(8)(0, 0), 0.1, 1E-16); } @@ -56,8 +55,8 @@ TEST(interface, output_sig_figs_2) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); EXPECT_NEAR(chains.samples(8)(0, 0), 0.12, 1E-16); } @@ -84,7 +83,7 @@ TEST(interface, output_sig_figs_9) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); EXPECT_NEAR(chains.samples(8)(0, 0), 0.123456789, 1E-16); } diff --git a/src/test/interface/stansummary_test.cpp b/src/test/interface/stansummary_test.cpp index 034cff84ba..4876333df1 100644 --- a/src/test/interface/stansummary_test.cpp +++ b/src/test/interface/stansummary_test.cpp @@ -136,10 +136,10 @@ TEST(CommandStansummary, matrix_index_2d) { TEST(CommandStansummary, header_tests) { std::string expect - = " Mean MCSE StdDev 10% 50% 90% N_Eff " - "N_Eff/s R_hat\n"; + = " Mean MCSE StdDev MAD 10% 50% 90%" + " ESS_bulk ESS_tail R_hat\n"; std::string expect_csv - = "name,Mean,MCSE,StdDev,10%,50%,90%,N_Eff,N_Eff/s,R_hat\n"; + = "name,Mean,MCSE,StdDev,MAD,10%,50%,90%,ESS_bulk,ESS_tail,R_hat\n"; std::vector pcts; pcts.push_back("10"); pcts.push_back("50"); @@ -150,12 +150,14 @@ TEST(CommandStansummary, header_tests) { EXPECT_FLOAT_EQ(probs[2], 0.9); std::vector header = get_header(pcts); - EXPECT_EQ(header.size(), pcts.size() + 6); + EXPECT_EQ(header.size(), pcts.size() + 7); EXPECT_EQ(header[0], "Mean"); EXPECT_EQ(header[2], "StdDev"); - EXPECT_EQ(header[4], "50%"); - EXPECT_EQ(header[6], "N_Eff"); - EXPECT_EQ(header[8], "R_hat"); + EXPECT_EQ(header[3], "MAD"); + EXPECT_EQ(header[5], "50%"); + EXPECT_EQ(header[7], "ESS_bulk"); + EXPECT_EQ(header[8], "ESS_tail"); + EXPECT_EQ(header[9], "R_hat"); Eigen::VectorXi column_widths(header.size()); for (size_t i = 0, w = 5; i < header.size(); ++i, ++w) { @@ -218,8 +220,8 @@ TEST(CommandStansummary, param_tests) { Eigen::VectorXd warmup_times(filenames.size()); Eigen::VectorXd sampling_times(filenames.size()); Eigen::VectorXi thin(filenames.size()); - stan::mcmc::chains<> chains = parse_csv_files( - filenames, metadata, warmup_times, sampling_times, thin, &std::cout); + auto chains = parse_csv_files(filenames, metadata, warmup_times, + sampling_times, thin, &std::cout); EXPECT_EQ(chains.num_chains(), 1); EXPECT_EQ(chains.num_params(), 8); @@ -237,7 +239,7 @@ TEST(CommandStansummary, param_tests) { size_t num_model_params = chains.num_params() - model_params_offset; EXPECT_EQ(num_model_params, 1); - Eigen::MatrixXd model_params(num_model_params, 9); + Eigen::MatrixXd model_params(num_model_params, 10); std::vector model_param_idxes(num_model_params); std::iota(model_param_idxes.begin(), model_param_idxes.end(), model_params_offset); @@ -246,7 +248,7 @@ TEST(CommandStansummary, param_tests) { double mean_theta = model_params(0, 0); EXPECT_TRUE(mean_theta > 0.25); EXPECT_TRUE(mean_theta < 0.27); - double rhat_theta = model_params(0, 8); + double rhat_theta = model_params(0, 9); EXPECT_TRUE(rhat_theta > 0.999); EXPECT_TRUE(rhat_theta < 1.01); } @@ -386,7 +388,7 @@ TEST(CommandStansummary, bad_percentiles_arg) { ASSERT_TRUE(out.hasError) << "\"" << out.command << "\" failed to quit with an error"; - arg_percentiles = "--percentiles \"0,100\""; + arg_percentiles = "--percentiles \"101\""; out = run_command(command + " " + arg_percentiles + " " + csv_file); EXPECT_TRUE(boost::algorithm::contains(out.output, expected_message)); ASSERT_TRUE(out.hasError) @@ -423,17 +425,17 @@ TEST(CommandStansummary, bad_include_param_args) { TEST(CommandStansummary, check_console_output) { std::string lp - = "lp__ -7.3 3.7e-02 0.77 -9.1 -7.0 -6.8 443 " - "19275 1.0"; + = "lp__ -7.3 3.7e-02 0.77 0.30 -9.0 -7.0 -6.8 " + " 519 503 1.0"; std::string theta - = "theta 0.26 6.1e-03 0.12 0.079 0.25 0.47 384 " - "16683 1.00"; + = "theta 0.26 6.1e-03 0.12 0.12 0.080 0.25 0.47 " + " 362 396 1.0"; std::string accept_stat - = "accept_stat__ 0.90 4.6e-03 1.5e-01 0.57 0.96 1.0 1026 " - "44597 1.00"; + = "accept_stat__ 0.90 4.6e-03 1.5e-01 0.064 0.57 0.96 1.0 " + "1284 941 1.00"; std::string energy - = "energy__ 7.8 5.1e-02 1.0e+00 6.8 7.5 9.9 411 " - "17865 1.0"; + = "energy__ 7.8 5.1e-02 1.0e+00 0.75 6.8 7.5 9.9 " + " 490 486 1.0"; std::string path_separator; path_separator.push_back(get_path_separator()); @@ -478,17 +480,16 @@ TEST(CommandStansummary, check_console_output) { TEST(CommandStansummary, check_csv_output) { std::string csv_header - = "name,Mean,MCSE,StdDev,5%,50%,95%,N_Eff,N_Eff/s,R_hat"; + = "name,Mean,MCSE,StdDev,MAD,5%,50%,95%,ESS_bulk,ESS_tail,R_hat"; std::string lp - = "\"lp__\",-7.2719,0.0365168,0.768874,-9.05757,-6.96978,-6.75008,443." - "328,19275.1,1.00037"; + = "\"lp__\",-7.2719,0.0365168,0.768874,0.303688,-8.98426,-6.97009,-6." + "75007,519.29,503.309,1.00141"; std::string energy - = "\"energy__\"" - ",7.78428,0.0508815,1.0314,6.80383,7.46839,9.88601,410." - "898,17865.1,1.00075"; + = "\"energy__\",7.78428,0.0508815,1.0314,0.745859,6.80565,7.46758,9.8864," + "489.874,486.438,1.00495"; std::string theta - = "\"theta\",0.256552,0.00610844,0.119654,0.0786292,0.24996,0.470263,383." - "704,16682.8,0.999309"; + = "\"theta\",0.256552,0.00610844,0.119654,0.120965,0.0802982,0.24996,0." + "47034,361.506,395.736,1.00186"; std::string path_separator; path_separator.push_back(get_path_separator()); @@ -531,9 +532,9 @@ TEST(CommandStansummary, check_csv_output) { } TEST(CommandStansummary, check_csv_output_no_percentiles) { - std::string csv_header = "name,Mean,MCSE,StdDev,N_Eff,N_Eff/s,R_hat"; + std::string csv_header = "name,Mean,MCSE,StdDev,MAD,ESS_bulk,ESS_tail,R_hat"; std::string lp - = "\"lp__\",-7.2719,0.0365168,0.768874,443.328,19275.1,1.00037"; + = "\"lp__\",-7.2719,0.0365168,0.768874,0.303688,519.29,503.309,1.00141"; std::string path_separator; path_separator.push_back(get_path_separator()); @@ -570,11 +571,12 @@ TEST(CommandStansummary, check_csv_output_no_percentiles) { TEST(CommandStansummary, check_csv_output_sig_figs) { std::string csv_header - = "name,Mean,MCSE,StdDev,5%,50%,95%,N_Eff,N_Eff/s,R_hat"; - std::string lp = "\"lp__\",-7.3,0.037,0.77,-9.1,-7,-6.8,4.4e+02,1.9e+04,1"; - std::string energy = "\"energy__\",7.8,0.051,1,6.8,7.5,9.9,4.1e+02,1.8e+04,1"; + = "name,Mean,MCSE,StdDev,MAD,5%,50%,95%,ESS_bulk,ESS_tail,R_hat"; + std::string lp = "\"lp__\",-7.3,0.037,0.77,0.3,-9,-7,-6.8,5.2e+02,5e+02,1"; + std::string energy + = "\"energy__\",7.8,0.051,1,0.75,6.8,7.5,9.9,4.9e+02,4.9e+02,1"; std::string theta - = "\"theta\",0.26,0.0061,0.12,0.079,0.25,0.47,3.8e+02,1.7e+04,1"; + = "\"theta\",0.26,0.0061,0.12,0.12,0.08,0.25,0.47,3.6e+02,4e+02,1"; std::string path_separator; path_separator.push_back(get_path_separator()); @@ -619,20 +621,20 @@ TEST(CommandStansummary, check_csv_output_sig_figs) { TEST(CommandStansummary, check_csv_output_include_param) { std::string csv_header - = "name,Mean,MCSE,StdDev,5%,50%,95%,N_Eff,N_Eff/s,R_hat"; + = "name,Mean,MCSE,StdDev,MAD,5%,50%,95%,ESS_bulk,ESS_tail,R_hat"; std::string lp - = "\"lp__\",-15.5617,0.97319,6.05585,-25.3432,-15.7562,-5.48405,38.7217," - "372.539,1.00208"; + = "\"lp__\",-15.5617,0.97319,6.05585,6.3817,-25.3182,-15.7598,-5.47732," + "41.1897,113.537,1.00153"; std::string energy - = "\"energy__\",20.5888,1.01449,6.43127,10.2634,20.8487,30.9906,40.1879," - "386.645,1.00109"; + = "\"energy__\",20.5888,1.01449,6.43127,6.6161,10.2809,20.8278,30.9921," + "42.5605,140.171,1.00069"; // note: skipping theta 1-5 std::string theta6 - = "\"theta[6]\",5.001,0.365016,5.76072,-4.99688,5.23474,14.1597,249.074," - "2396.33,1.00109"; + = "\"theta[6]\",5.001,0.365016,5.76072,5.37947,-4.95375,5.22746,14.1688," + "230.645,464.978,1.00054"; std::string theta7 - = "\"theta[7]\",8.54125,0.650098,6.22195,-0.841225,8.09613,19.256,91." - "6001,881.279,0.999012"; + = "\"theta[7]\",8.54125,0.650098,6.22195,5.35785,-0.814388,8.09342,19." + "2622,92.3075,241.177,1.00244"; // note: skipping theta 8 std::string message = "# Inference for Stan model: eight_schools_cp_model";