diff --git a/src/stan/analyze/mcmc/mcse.hpp b/src/stan/analyze/mcmc/mcse.hpp index 6481796b72..f8789fb670 100644 --- a/src/stan/analyze/mcmc/mcse.hpp +++ b/src/stan/analyze/mcmc/mcse.hpp @@ -2,6 +2,7 @@ #define STAN_ANALYZE_MCMC_MCSE_HPP #include +#include #include #include #include @@ -42,11 +43,22 @@ inline double mcse_sd(const Eigen::MatrixXd& chains) { if (chains.rows() < 4 || !is_finite_and_varies(chains)) return std::numeric_limits::quiet_NaN(); - Eigen::MatrixXd diffs = (chains.array() - chains.mean()).matrix(); - double Evar = diffs.array().square().mean(); - double varvar = (math::mean(diffs.array().pow(4) - Evar * Evar)) - / ess(diffs.array().abs().matrix()); - return std::sqrt(varvar / Evar / 4); + // center the data, take abs value + Eigen::MatrixXd draws_ctr = (chains.array() - chains.mean()).abs().matrix(); + + // posterior pkg fn `ess_mean` computes on split chains + double ess_mean = ess(split_chains(draws_ctr)); + + // estimated variance (2nd moment) + double Evar = draws_ctr.array().square().mean(); + + // variance of variance, adjusted for ESS + double fourth_moment = draws_ctr.array().pow(4).mean(); + double varvar = (fourth_moment - std::pow(Evar, 2)) / ess_mean; + + // variance of standard deviation - use Taylor series approximation + double varsd = varvar / Evar / 4.0; + return std::sqrt(varsd); } } // namespace analyze diff --git a/src/test/unit/analyze/mcmc/ess_basic_test.cpp b/src/test/unit/analyze/mcmc/ess_basic_test.cpp index 8d43ce9217..d02c91d41d 100644 --- a/src/test/unit/analyze/mcmc/ess_basic_test.cpp +++ b/src/test/unit/analyze/mcmc/ess_basic_test.cpp @@ -7,29 +7,42 @@ #include #include -TEST(RankNormalizedEss, test_basic_ess) { - std::stringstream out; - Eigen::MatrixXd chains_lp(1000, 4); - Eigen::MatrixXd chains_theta(1000, 4); - - std::vector draws_theta(4); - std::vector draws_lp(4); - std::vector sizes(4); - - for (size_t i = 0; i < 4; ++i) { - std::stringstream fname; - fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) - << ".csv"; - std::ifstream bern_stream(fname.str(), std::ifstream::in); - stan::io::stan_csv bern_csv +class EssBasic : public testing::Test { +public: + void SetUp() { + chains_lp.resize(1000, 4); + chains_theta.resize(1000, 4); + draws_theta.resize(4); + draws_lp.resize(4); + sizes.resize(4); + for (size_t i = 0; i < 4; ++i) { + std::stringstream fname; + fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) + << ".csv"; + std::ifstream bern_stream(fname.str(), std::ifstream::in); + stan::io::stan_csv bern_csv = stan::io::stan_csv_reader::parse(bern_stream, &out); - bern_stream.close(); - chains_lp.col(i) = bern_csv.samples.col(0); - chains_theta.col(i) = bern_csv.samples.col(7); - draws_lp[i] = chains_lp.col(i).data(); - draws_theta[i] = chains_theta.col(i).data(); - sizes[i] = 1000; + bern_stream.close(); + chains_lp.col(i) = bern_csv.samples.col(0); + chains_theta.col(i) = bern_csv.samples.col(7); + draws_lp[i] = chains_lp.col(i).data(); + draws_theta[i] = chains_theta.col(i).data(); + sizes[i] = 1000; + } + } + + void TearDown() { } + + std::stringstream out; + Eigen::MatrixXd chains_lp; + Eigen::MatrixXd chains_theta; + std::vector draws_theta; + std::vector draws_lp; + std::vector sizes; +}; + +TEST_F(EssBasic, test_basic_ess) { double ess_lp_expect = 1335.4137; double ess_theta_expect = 1377.503; diff --git a/src/test/unit/analyze/mcmc/mcse_test.cpp b/src/test/unit/analyze/mcmc/mcse_test.cpp new file mode 100644 index 0000000000..1fffd81fe9 --- /dev/null +++ b/src/test/unit/analyze/mcmc/mcse_test.cpp @@ -0,0 +1,72 @@ +#include +#include +#include +#include +#include +#include +#include + +class MonteCarloStandardError : public testing::Test { +public: + void SetUp() { + chains_lp.resize(1000, 4); + chains_theta.resize(1000, 4); + chains_divergent.resize(1000, 4); + for (size_t i = 0; i < 4; ++i) { + std::stringstream fname; + fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) + << ".csv"; + std::ifstream bern_stream(fname.str(), std::ifstream::in); + stan::io::stan_csv bern_csv + = stan::io::stan_csv_reader::parse(bern_stream, &out); + bern_stream.close(); + chains_lp.col(i) = bern_csv.samples.col(0); + chains_theta.col(i) = bern_csv.samples.col(7); + chains_divergent.col(i) = bern_csv.samples.col(5); + } + } + + void TearDown() { + } + + std::stringstream out; + Eigen::MatrixXd chains_lp; + Eigen::MatrixXd chains_theta; + Eigen::MatrixXd chains_divergent; +}; + +TEST_F(MonteCarloStandardError, test_mcse) { + double mcse_mean_lp_expect = 0.020164778; + double mcse_mean_theta_expect = 0.0032339916; + + double mcse_sd_lp_expect = 0.0355305; + double mcse_sd_theta_expect = 0.0021642137; + EXPECT_NEAR(mcse_mean_lp_expect, stan::analyze::mcse_mean(chains_lp), 0.0001); + EXPECT_NEAR(mcse_mean_theta_expect, stan::analyze::mcse_mean(chains_theta), 0.0001); + + EXPECT_NEAR(mcse_sd_lp_expect, stan::analyze::mcse_sd(chains_lp), 0.0001); + EXPECT_NEAR(mcse_sd_theta_expect, stan::analyze::mcse_sd(chains_theta), 0.0001); +} + +TEST_F(MonteCarloStandardError, const_fail) { + auto mcse_mean = stan::analyze::mcse_mean(chains_divergent); + auto mcse_sd = stan::analyze::mcse_sd(chains_divergent); + EXPECT_TRUE(std::isnan(mcse_mean)); + EXPECT_TRUE(std::isnan(mcse_sd)); +} + +TEST_F(MonteCarloStandardError, inf_fail) { + chains_theta(0,0) = std::numeric_limits::infinity(); + auto mcse_mean = stan::analyze::mcse_mean(chains_theta); + auto mcse_sd = stan::analyze::mcse_sd(chains_theta); + EXPECT_TRUE(std::isnan(mcse_mean)); + EXPECT_TRUE(std::isnan(mcse_sd)); +} + +TEST_F(MonteCarloStandardError, short_chains_fail) { + chains_theta.resize(3, 4); + auto mcse_mean = stan::analyze::mcse_mean(chains_theta); + auto mcse_sd = stan::analyze::mcse_sd(chains_theta); + EXPECT_TRUE(std::isnan(mcse_mean)); + EXPECT_TRUE(std::isnan(mcse_sd)); +} diff --git a/src/test/unit/analyze/mcmc/rhat_basic_test.cpp b/src/test/unit/analyze/mcmc/rhat_basic_test.cpp new file mode 100644 index 0000000000..15db985030 --- /dev/null +++ b/src/test/unit/analyze/mcmc/rhat_basic_test.cpp @@ -0,0 +1,61 @@ +#include +#include +#include +#include +#include +#include + +class RhatBasic : public testing::Test { +public: + void SetUp() { + chains_lp.resize(1000, 4); + chains_theta.resize(1000, 4); + draws_theta.resize(4); + draws_lp.resize(4); + sizes.resize(4); + for (size_t i = 0; i < 4; ++i) { + std::stringstream fname; + fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) + << ".csv"; + std::ifstream bern_stream(fname.str(), std::ifstream::in); + stan::io::stan_csv bern_csv + = stan::io::stan_csv_reader::parse(bern_stream, &out); + bern_stream.close(); + chains_lp.col(i) = bern_csv.samples.col(0); + chains_theta.col(i) = bern_csv.samples.col(7); + draws_lp[i] = chains_lp.col(i).data(); + draws_theta[i] = chains_theta.col(i).data(); + sizes[i] = 1000; + } + } + + void TearDown() { + } + + std::stringstream out; + Eigen::MatrixXd chains_lp; + Eigen::MatrixXd chains_theta; + std::vector draws_theta; + std::vector draws_lp; + std::vector sizes; +}; + +TEST_F(RhatBasic, test_basic_rhat) { + double rhat_lp_basic_expect = 1.0001296; + double rhat_theta_basic_expect = 1.0029197; + + auto rhat_basic_lp = stan::analyze::rhat(chains_lp); + auto old_rhat_basic_lp + = stan::analyze::compute_potential_scale_reduction(draws_lp, sizes); + + auto rhat_basic_theta = stan::analyze::rhat(chains_theta); + auto old_rhat_basic_theta + = stan::analyze::compute_potential_scale_reduction(draws_theta, sizes); + + EXPECT_NEAR(rhat_lp_basic_expect, rhat_basic_lp, 0.00001); + EXPECT_NEAR(rhat_theta_basic_expect, rhat_basic_theta, 0.00001); + + EXPECT_NEAR(old_rhat_basic_lp, rhat_basic_lp, 0.00001); + EXPECT_NEAR(old_rhat_basic_theta, rhat_basic_theta, 0.00001); +} + diff --git a/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp b/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp index ac52cf99db..0b4ecb6d25 100644 --- a/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp +++ b/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp @@ -6,29 +6,36 @@ #include #include -TEST(RankNormalizedEss, test_bulk_tail_ess) { - std::stringstream out; - Eigen::MatrixXd chains_lp(1000, 4); - Eigen::MatrixXd chains_theta(1000, 4); - - std::vector draws_theta(4); - std::vector draws_lp(4); - std::vector sizes(4); - - for (size_t i = 0; i < 4; ++i) { - std::stringstream fname; - fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) - << ".csv"; - std::ifstream bern_stream(fname.str(), std::ifstream::in); - stan::io::stan_csv bern_csv +class RankNormalizedEss : public testing::Test { +public: + void SetUp() { + chains_lp.resize(1000, 4); + chains_theta.resize(1000, 4); + chains_divergent.resize(1000, 4); + for (size_t i = 0; i < 4; ++i) { + std::stringstream fname; + fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) + << ".csv"; + std::ifstream bern_stream(fname.str(), std::ifstream::in); + stan::io::stan_csv bern_csv = stan::io::stan_csv_reader::parse(bern_stream, &out); - bern_stream.close(); - chains_lp.col(i) = bern_csv.samples.col(0); - chains_theta.col(i) = bern_csv.samples.col(7); - draws_lp[i] = chains_lp.col(i).data(); - draws_theta[i] = chains_theta.col(i).data(); - sizes[i] = 1000; + bern_stream.close(); + chains_lp.col(i) = bern_csv.samples.col(0); + chains_theta.col(i) = bern_csv.samples.col(7); + chains_divergent.col(i) = bern_csv.samples.col(5); + } } + + void TearDown() { + } + + std::stringstream out; + Eigen::MatrixXd chains_lp; + Eigen::MatrixXd chains_theta; + Eigen::MatrixXd chains_divergent; +}; + +TEST_F(RankNormalizedEss, test_bulk_tail_ess) { double ess_lp_bulk_expect = 1512.7684; double ess_lp_tail_expect = 1591.9707; @@ -45,7 +52,20 @@ TEST(RankNormalizedEss, test_bulk_tail_ess) { EXPECT_NEAR(ess_theta_tail_expect, ess_theta.second, 0.001); } -TEST(RankNormalizedEss, short_chains_fail) { +TEST_F(RankNormalizedEss, const_fail) { + auto ess = stan::analyze::split_rank_normalized_ess(chains_divergent); + EXPECT_TRUE(std::isnan(ess.first)); + EXPECT_TRUE(std::isnan(ess.second)); +} + +TEST_F(RankNormalizedEss, inf_fail) { + chains_theta(0,0) = std::numeric_limits::infinity(); + auto ess = stan::analyze::split_rank_normalized_ess(chains_theta); + EXPECT_TRUE(std::isnan(ess.first)); + EXPECT_TRUE(std::isnan(ess.second)); +} + +TEST_F(RankNormalizedEss, short_chains_fail) { std::stringstream out; std::ifstream eight_schools_5iters_1_stream, eight_schools_5iters_2_stream; stan::io::stan_csv eight_schools_5iters_1, eight_schools_5iters_2; @@ -72,29 +92,3 @@ TEST(RankNormalizedEss, short_chains_fail) { } } -TEST(RankNormalizedEss, const_fail) { - std::stringstream out; - std::ifstream bernoulli_const_1_stream, bernoulli_const_2_stream; - stan::io::stan_csv bernoulli_const_1, bernoulli_const_2; - bernoulli_const_1_stream.open( - "src/test/unit/mcmc/test_csv_files/bernoulli_const_1.csv", - std::ifstream::in); - bernoulli_const_1 - = stan::io::stan_csv_reader::parse(bernoulli_const_1_stream, &out); - bernoulli_const_1_stream.close(); - bernoulli_const_2_stream.open( - "src/test/unit/mcmc/test_csv_files/bernoulli_const_2.csv", - std::ifstream::in); - bernoulli_const_2 - = stan::io::stan_csv_reader::parse(bernoulli_const_2_stream, &out); - bernoulli_const_2_stream.close(); - - Eigen::MatrixXd chains(bernoulli_const_1.samples.rows(), 2); - chains.col(0) - = bernoulli_const_1.samples.col(bernoulli_const_1.samples.cols() - 1); - chains.col(1) - = bernoulli_const_2.samples.col(bernoulli_const_2.samples.cols() - 1); - auto ess = stan::analyze::split_rank_normalized_ess(chains); - EXPECT_TRUE(std::isnan(ess.first)); - EXPECT_TRUE(std::isnan(ess.second)); -} diff --git a/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp b/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp index a1ca5629c4..6cdc16c2e1 100644 --- a/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp +++ b/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp @@ -1,4 +1,3 @@ -#include #include #include #include @@ -7,81 +6,59 @@ #include #include #include +#include -TEST(RankNormalizedRhat, test_basic_bulk_tail_rhat) { - std::stringstream out; - Eigen::MatrixXd chains_lp(1000, 4); - Eigen::MatrixXd chains_theta(1000, 4); - - std::vector draws_theta(4); - std::vector draws_lp(4); - std::vector sizes(4); - - for (size_t i = 0; i < 4; ++i) { - std::stringstream fname; - fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) - << ".csv"; - std::ifstream bern_stream(fname.str(), std::ifstream::in); - stan::io::stan_csv bern_csv +class RankNormalizedRhat : public testing::Test { +public: + void SetUp() { + chains_lp.resize(1000, 4); + chains_theta.resize(1000, 4); + chains_divergent.resize(1000, 4); + for (size_t i = 0; i < 4; ++i) { + std::stringstream fname; + fname << "src/test/unit/analyze/mcmc/test_csv_files/bern" << (i + 1) + << ".csv"; + std::ifstream bern_stream(fname.str(), std::ifstream::in); + stan::io::stan_csv bern_csv = stan::io::stan_csv_reader::parse(bern_stream, &out); - bern_stream.close(); - chains_lp.col(i) = bern_csv.samples.col(0); - chains_theta.col(i) = bern_csv.samples.col(7); - draws_lp[i] = chains_lp.col(i).data(); - draws_theta[i] = chains_theta.col(i).data(); - sizes[i] = 1000; + bern_stream.close(); + chains_lp.col(i) = bern_csv.samples.col(0); + chains_theta.col(i) = bern_csv.samples.col(7); + chains_divergent.col(i) = bern_csv.samples.col(5); + } } - double rhat_lp_basic_expect = 1.0001296; - double rhat_lp_new_expect = 1.0007301; - - double rhat_theta_basic_expect = 1.0029197; - double rhat_theta_new_expect = 1.0067897; - auto rhat_basic_lp = stan::analyze::rhat(chains_lp); - auto old_rhat_basic_lp - = stan::analyze::compute_potential_scale_reduction(draws_lp, sizes); - auto rhat_lp = stan::analyze::split_rank_normalized_rhat(chains_lp); + void TearDown() { + } - auto rhat_basic_theta = stan::analyze::rhat(chains_theta); - auto old_rhat_basic_theta - = stan::analyze::compute_potential_scale_reduction(draws_theta, sizes); - auto rhat_theta = stan::analyze::split_rank_normalized_rhat(chains_theta); + std::stringstream out; + Eigen::MatrixXd chains_lp; + Eigen::MatrixXd chains_theta; + Eigen::MatrixXd chains_divergent; +}; - EXPECT_NEAR(rhat_lp_basic_expect, rhat_basic_lp, 0.00001); - EXPECT_NEAR(rhat_theta_basic_expect, rhat_basic_theta, 0.00001); +TEST_F(RankNormalizedRhat, test_bulk_tail_rhat) { + double rhat_lp_expect = 1.0007301; + double rhat_theta_expect = 1.0067897; - EXPECT_NEAR(old_rhat_basic_lp, rhat_basic_lp, 0.00001); - EXPECT_NEAR(old_rhat_basic_theta, rhat_basic_theta, 0.00001); + auto rhat_lp = stan::analyze::split_rank_normalized_rhat(chains_lp); + auto rhat_theta = stan::analyze::split_rank_normalized_rhat(chains_theta); - EXPECT_NEAR(rhat_lp_new_expect, std::max(rhat_lp.first, rhat_lp.second), + EXPECT_NEAR(rhat_lp_expect, std::max(rhat_lp.first, rhat_lp.second), 0.00001); - EXPECT_NEAR(rhat_theta_new_expect, + EXPECT_NEAR(rhat_theta_expect, std::max(rhat_theta.first, rhat_theta.second), 0.00001); } -TEST(RankNormalizedRhat, const_fail) { - std::stringstream out; - std::ifstream bernoulli_const_1_stream, bernoulli_const_2_stream; - stan::io::stan_csv bernoulli_const_1, bernoulli_const_2; - bernoulli_const_1_stream.open( - "src/test/unit/mcmc/test_csv_files/bernoulli_const_1.csv", - std::ifstream::in); - bernoulli_const_1 - = stan::io::stan_csv_reader::parse(bernoulli_const_1_stream, &out); - bernoulli_const_1_stream.close(); - bernoulli_const_2_stream.open( - "src/test/unit/mcmc/test_csv_files/bernoulli_const_2.csv", - std::ifstream::in); - bernoulli_const_2 - = stan::io::stan_csv_reader::parse(bernoulli_const_2_stream, &out); - bernoulli_const_2_stream.close(); +TEST_F(RankNormalizedRhat, const_fail) { + auto rhat = stan::analyze::split_rank_normalized_rhat(chains_divergent); + EXPECT_TRUE(std::isnan(rhat.first)); + EXPECT_TRUE(std::isnan(rhat.second)); +} - Eigen::MatrixXd chains(bernoulli_const_1.samples.rows(), 2); - chains.col(0) - = bernoulli_const_1.samples.col(bernoulli_const_1.samples.cols() - 1); - chains.col(1) - = bernoulli_const_2.samples.col(bernoulli_const_2.samples.cols() - 1); - auto rhat = stan::analyze::split_rank_normalized_rhat(chains); +TEST_F(RankNormalizedRhat, inf_fail) { + chains_theta(0,0) = std::numeric_limits::infinity(); + auto rhat = stan::analyze::split_rank_normalized_rhat(chains_theta); EXPECT_TRUE(std::isnan(rhat.first)); EXPECT_TRUE(std::isnan(rhat.second)); }