From 72eb95060662885a82a218d0c296c3c9e87b6286 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sun, 10 Sep 2017 14:14:36 -0700 Subject: [PATCH 01/12] Changed vector_rng_test_helper tools to work with anywhere between 1-3 parameter rngs (Issue 45) Vectorized all unbounded continuous random number generators and added tests (Issue 45) Changed how skew_normal random numbers are generated (following https://arxiv.org/pdf/0911.2342.pdf) (Issue 527) --- stan/math/prim/scal/prob/cauchy_rng.hpp | 53 +- .../prim/scal/prob/double_exponential_rng.hpp | 65 ++- .../prim/scal/prob/exp_mod_normal_rng.hpp | 63 ++- stan/math/prim/scal/prob/gumbel_rng.hpp | 53 +- stan/math/prim/scal/prob/logistic_rng.hpp | 56 +- stan/math/prim/scal/prob/skew_normal_rng.hpp | 73 ++- stan/math/prim/scal/prob/student_t_rng.hpp | 70 ++- test/unit/math/prim/mat/prob/cauchy_test.cpp | 49 ++ .../prim/mat/prob/double_exponential_test.cpp | 49 ++ .../prim/mat/prob/exp_mod_normal_test.cpp | 18 + test/unit/math/prim/mat/prob/gumbel_test.cpp | 49 ++ .../unit/math/prim/mat/prob/logistic_test.cpp | 49 ++ test/unit/math/prim/mat/prob/normal_test.cpp | 24 +- .../math/prim/mat/prob/skew_normal_test.cpp | 51 ++ .../math/prim/mat/prob/student_t_test.cpp | 48 ++ .../prim/mat/prob/vector_rng_test_helper.hpp | 530 ++++++++++++++---- test/unit/math/prim/scal/prob/util.hpp | 34 +- 17 files changed, 1058 insertions(+), 276 deletions(-) create mode 100644 test/unit/math/prim/mat/prob/cauchy_test.cpp create mode 100644 test/unit/math/prim/mat/prob/double_exponential_test.cpp create mode 100644 test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp create mode 100644 test/unit/math/prim/mat/prob/gumbel_test.cpp create mode 100644 test/unit/math/prim/mat/prob/logistic_test.cpp create mode 100644 test/unit/math/prim/mat/prob/skew_normal_test.cpp create mode 100644 test/unit/math/prim/mat/prob/student_t_test.cpp diff --git a/stan/math/prim/scal/prob/cauchy_rng.hpp b/stan/math/prim/scal/prob/cauchy_rng.hpp index 2e7fd2ca0a9..81bca55b142 100644 --- a/stan/math/prim/scal/prob/cauchy_rng.hpp +++ b/stan/math/prim/scal/prob/cauchy_rng.hpp @@ -1,51 +1,66 @@ #ifndef STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP #define STAN_MATH_PRIM_SCAL_PROB_CAUCHY_RNG_HPP -#include -#include #include #include -#include #include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include namespace stan { namespace math { - /** * Return a pseudorandom Cauchy variate for the given location and scale * using the specified random number generator. * + * mu and sigma can each be either scalars or vectors. All vector inputs + * must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter * @tparam RNG type of random number generator - * @param mu location parameter - * @param sigma positive scale parameter + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) * @param rng random number generator * @return Cauchy random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive + * @throw std::invalid_argument if vector arguments are not the same length */ - template - inline double - cauchy_rng(double mu, - double sigma, + template + inline typename VectorBuilder::type + cauchy_rng(const T_loc &mu, + const T_scale &sigma, RNG& rng) { using boost::variate_generator; using boost::random::cauchy_distribution; static const std::string function = "cauchy_rng"; + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", sigma); - variate_generator > - cauchy_rng(rng, cauchy_distribution<>(mu, sigma)); - return cauchy_rng(); - } + size_t N = max_size(mu, sigma); + + VectorBuilder output(N); + for (size_t n = 0; n < N; n++) { + variate_generator > + cauchy_rng(rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n])); + output[n] = cauchy_rng(); + } + + return output.data(); + } } } #endif diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 3c25118967c..9fbf182d6ba 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -1,58 +1,73 @@ #ifndef STAN_MATH_PRIM_SCAL_PROB_DOUBLE_EXPONENTIAL_RNG_HPP #define STAN_MATH_PRIM_SCAL_PROB_DOUBLE_EXPONENTIAL_RNG_HPP -#include -#include #include #include -#include #include -#include #include - -#include -#include -#include +#include +#include +#include +#include +#include #include namespace stan { namespace math { - /** * Return a pseudorandom double exponential variate with the given location * and scale using the specified random number generator. * + * mu and sigma can each be either scalars or vectors. All vector inputs + * must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter * @tparam RNG class of random number generator - * @param mu location parameter - * @param sigma positive scale parameter + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) * @param rng random number generator * @return double exponential random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive + * @throw std::invalid_argument if vector arguments are not the same length */ - template - inline double - double_exponential_rng(double mu, - double sigma, + template + inline typename VectorBuilder::type + double_exponential_rng(const T_loc &mu, + const T_scale &sigma, RNG& rng) { static const std::string function = "double_exponential_rng"; using boost::variate_generator; using boost::random::uniform_01; - using std::log; using std::abs; + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", sigma); + + size_t N = max_size(mu, sigma); + + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + variate_generator > rng_unit_01(rng, uniform_01<>()); + double a = 0; + double laplaceRN = rng_unit_01(); + if (0.5 - laplaceRN > 0) + a = 1.0; + else if (0.5 - laplaceRN < 0) + a = -1.0; + output[n] = mu_vec[n] - + sigma_vec[n] * a * log1m(2 * abs(0.5 - laplaceRN)); + } - variate_generator > - rng_unit_01(rng, uniform_01<>()); - double a = 0; - double laplaceRN = rng_unit_01(); - if (0.5 - laplaceRN > 0) - a = 1.0; - else if (0.5 - laplaceRN < 0) - a = -1.0; - return mu - sigma * a * log1m(2 * abs(0.5 - laplaceRN)); + return output.data(); } } } diff --git a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp index 7f815eb9eef..4d52d8f1115 100644 --- a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp +++ b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp @@ -3,36 +3,69 @@ #include #include -#include #include -#include -#include +#include #include -#include #include -#include -#include +#include #include #include #include namespace stan { namespace math { - - template - inline double - exp_mod_normal_rng(double mu, - double sigma, - double lambda, + /** + * Return a pseudorandom Exponentially modified normal variate for the + * given location, scale, and inverse scale using the specified random + * number generator. + * + * mu, sigma, and lambda can each be either scalars or vectors. All vector + * inputs must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter + * @tparam T_inv_scale Type of inverse scale parameter + * @tparam RNG type of random number generator + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) + * @param lambda (Sequence of) inverse scale parameter(s) + * @param rng random number generator + * @return Exponentially modified normal random variate + * @throw std::domain_error if mu is infinite, sigma is nonpositive, + * or lambda is nonpositive + * @throw std::invalid_argument if vector arguments are not the same length + */ + template + inline + typename VectorBuilder::type + exp_mod_normal_rng(const T_loc& mu, + const T_scale& sigma, + const T_inv_scale& lambda, RNG& rng) { static const std::string function = "exp_mod_normal_rng"; + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + scalar_seq_view lambda_vec(lambda); + check_finite(function, "Location parameter", mu); - check_positive_finite(function, "Inv_scale parameter", lambda); check_positive_finite(function, "Scale parameter", sigma); + check_positive_finite(function, "Inv_scale parameter", lambda); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", sigma, + "Inv_scale Parameter", lambda); + + size_t N = max_size(mu, sigma, lambda); + + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + output[n] = normal_rng(mu_vec[n], sigma_vec[n], rng) + + exponential_rng(lambda_vec[n], rng); + } - return normal_rng(mu, sigma, rng) - + exponential_rng(lambda, rng); + return output.data(); } } diff --git a/stan/math/prim/scal/prob/gumbel_rng.hpp b/stan/math/prim/scal/prob/gumbel_rng.hpp index 13a6ded98d3..4fac82e770d 100644 --- a/stan/math/prim/scal/prob/gumbel_rng.hpp +++ b/stan/math/prim/scal/prob/gumbel_rng.hpp @@ -3,36 +3,37 @@ #include #include -#include -#include -#include +#include +#include +#include #include -#include -#include -#include -#include #include #include #include namespace stan { namespace math { - /** * Return a pseudorandom Gumbel variate with the given location and scale * using the specified random number generator. * + * mu and sigma can each be either scalars or vectors. All vector inputs + * must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter * @tparam RNG type of random number generator - * @param mu location parameter - * @param beta positive scale parameter + * @param mu (Sequence of) location parameter(s) + * @param beta (Sequence of) scale parameter(s) * @param rng random number generator * @return Gumbel random variate * @throw std::domain_error if mu is infinite or beta is nonpositive. + * @throw std::invalid_argument if vector arguments are not the same length */ - template - inline double - gumbel_rng(double mu, - double beta, + template + inline typename VectorBuilder::type + gumbel_rng(const T_loc &mu, + const T_scale &beta, RNG& rng) { using boost::variate_generator; using boost::uniform_01; @@ -40,13 +41,27 @@ namespace stan { static const std::string function = "gumbel_rng"; check_finite(function, "Location parameter", mu); - check_positive(function, "Scale parameter", beta); + check_positive_finite(function, "Scale parameter", beta); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", beta); - variate_generator > - uniform01_rng(rng, uniform_01<>()); - return mu - beta * std::log(-std::log(uniform01_rng())); - } + scalar_seq_view mu_vec(mu); + scalar_seq_view beta_vec(beta); + + size_t N = max_size(mu, beta); + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + variate_generator > + uniform01_rng(rng, uniform_01<>()); + output[n] = mu_vec[n] - + beta_vec[n] * std::log(-std::log(uniform01_rng())); + } + + return output.data(); + } } } #endif diff --git a/stan/math/prim/scal/prob/logistic_rng.hpp b/stan/math/prim/scal/prob/logistic_rng.hpp index 870087eea54..da1fb64012a 100644 --- a/stan/math/prim/scal/prob/logistic_rng.hpp +++ b/stan/math/prim/scal/prob/logistic_rng.hpp @@ -3,26 +3,37 @@ #include #include -#include #include -#include -#include -#include -#include +#include +#include #include -#include -#include #include #include #include namespace stan { namespace math { - - template - inline double - logistic_rng(double mu, - double sigma, + /** + * Return a pseudorandom Logistic variate for the given location and scale + * using the specified random number generator. + * + * mu and sigma can each be either scalars or vectors. All vector inputs + * must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter + * @tparam RNG type of random number generator + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) + * @param rng random number generator + * @return Logistic random variate + * @throw std::domain_error if mu is infinite or sigma is nonpositive + * @throw std::invalid_argument if vector arguments are not the same length + */ + template + inline typename VectorBuilder::type + logistic_rng(const T_loc &mu, + const T_scale &sigma, RNG& rng) { using boost::variate_generator; using boost::random::exponential_distribution; @@ -31,12 +42,25 @@ namespace stan { check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", sigma); - variate_generator > - exp_rng(rng, exponential_distribution<>(1)); - return mu - sigma * std::log(exp_rng() / exp_rng()); - } + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + + size_t N = max_size(mu, sigma); + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + variate_generator > + exp_rng(rng, exponential_distribution<>(1)); + output[n] = mu_vec[n] - sigma_vec[n] * std::log(exp_rng() / exp_rng()); + } + + return output.data(); + } } } #endif diff --git a/stan/math/prim/scal/prob/skew_normal_rng.hpp b/stan/math/prim/scal/prob/skew_normal_rng.hpp index f119ef678fc..2665a093192 100644 --- a/stan/math/prim/scal/prob/skew_normal_rng.hpp +++ b/stan/math/prim/scal/prob/skew_normal_rng.hpp @@ -1,39 +1,72 @@ #ifndef STAN_MATH_PRIM_SCAL_PROB_SKEW_NORMAL_RNG_HPP #define STAN_MATH_PRIM_SCAL_PROB_SKEW_NORMAL_RNG_HPP -#include -#include +#include #include -#include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include #include namespace stan { namespace math { - - template - inline double - skew_normal_rng(double mu, - double sigma, - double alpha, + /** + * Return a pseudorandom Skew-normal variate for the given location, scale, + * and shape using the specified random number generator. + * + * mu, sigma, and alpha can each be either scalars or vectors. All vector + * inputs must be the same length. + * + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter + * @tparam T_shape Type of shape parameter + * @tparam RNG type of random number generator + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) + * @param alpha (Sequence of) shape parameter(s) + * @param rng random number generator + * @return Skew-normal random variate + * @throw std::domain_error if mu is infinite, sigma is nonpositive, or + * alpha is infinite + * @throw std::invalid_argument if vector arguments are not the same length + */ + template + inline typename VectorBuilder::type + skew_normal_rng(const T_loc& mu, + const T_scale& sigma, + const T_shape& alpha, RNG& rng) { - boost::math::skew_normal_distribution<> dist(mu, sigma, alpha); - static const std::string function = "skew_normal_rng"; check_finite(function, "Location parameter", mu); + check_positive_finite(function, "Scale parameter", sigma); check_finite(function, "Shape parameter", alpha); - check_positive(function, "Scale parameter", sigma); + check_consistent_sizes(function, + "Location parameter", mu, + "Scale Parameter", sigma, + "Shape Parameter", alpha); - return quantile(dist, uniform_rng(0.0, 1.0, rng)); - } + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + scalar_seq_view alpha_vec(alpha); + + size_t N = max_size(mu, sigma, alpha); + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + double r1 = normal_rng(0.0, 1.0, rng), + r2 = normal_rng(0.0, 1.0, rng); + + if (r2 > alpha_vec[n] * r1) + r1 *= -1; + + output[n] = mu_vec[n] + sigma_vec[n] * r1; + } + + return output.data(); + } } } #endif diff --git a/stan/math/prim/scal/prob/student_t_rng.hpp b/stan/math/prim/scal/prob/student_t_rng.hpp index dfb0d3e277b..ed2d4543d5b 100644 --- a/stan/math/prim/scal/prob/student_t_rng.hpp +++ b/stan/math/prim/scal/prob/student_t_rng.hpp @@ -1,46 +1,72 @@ #ifndef STAN_MATH_PRIM_SCAL_PROB_STUDENT_T_RNG_HPP #define STAN_MATH_PRIM_SCAL_PROB_STUDENT_T_RNG_HPP -#include -#include #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include #include +#include +#include #include namespace stan { namespace math { - - template - inline double - student_t_rng(double nu, - double mu, - double sigma, + /** + * Return a pseudorandom student-t variate for the given degrees of freedom, + * location, and scale using the specified random number generator. + * + * nu, mu, and sigma can each be either scalars or vectors. All vector + * inputs must be the same length. + * + * @tparam T_deg Type of degrees of freedom parameter + * @tparam T_loc Type of location parameter + * @tparam T_scale Type of scale parameter + * @tparam RNG type of random number generator + * @param nu (Sequence of) degrees of freedom parameter(s) + * @param mu (Sequence of) location parameter(s) + * @param sigma (Sequence of) scale parameter(s) + * @param rng random number generator + * @return Student-t random variate + * @throw std::domain_error if nu is nonpositive, mu is infinite, or sigma + * is nonpositive + * @throw std::invalid_argument if vector arguments are not the same length + */ + template + inline typename VectorBuilder::type + student_t_rng(const T_deg& nu, + const T_loc& mu, + const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::student_t_distribution; static const std::string function = "student_t_rng"; + scalar_seq_view nu_vec(nu); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + check_positive_finite(function, "Degrees of freedom parameter", nu); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, + "Degrees of freedom parameter", nu, + "Location parameter", mu, + "Scale Parameter", sigma); + + size_t N = max_size(nu, mu, sigma); + + VectorBuilder output(N); + + for (size_t n = 0; n < N; n++) { + variate_generator > + rng_unit_student_t(rng, student_t_distribution<>(nu_vec[n])); + output[n] = mu_vec[n] + sigma_vec[n] * rng_unit_student_t(); + } - variate_generator > - rng_unit_student_t(rng, student_t_distribution<>(nu)); - return mu + sigma * rng_unit_student_t(); + return output.data(); } } diff --git a/test/unit/math/prim/mat/prob/cauchy_test.cpp b/test/unit/math/prim/mat/prob/cauchy_test.cpp new file mode 100644 index 00000000000..bd43b9a23e3 --- /dev/null +++ b/test/unit/math/prim/mat/prob/cauchy_test.cpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct cauchy_rng_wrapper { + template + auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + return stan::math::cauchy_rng(mu, sigma, rng); + } +}; + +std::vector build_quantiles(int N, double mu, double sigma, + double unused) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + boost::math::cauchy_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsCauchy, errorCheck) { + check_dist_throws_all_types(cauchy_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {}, {}); +} + +TEST(ProbDistributionsCauchy, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, cauchy_rng_wrapper{}, build_quantiles, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}, + {}); +} diff --git a/test/unit/math/prim/mat/prob/double_exponential_test.cpp b/test/unit/math/prim/mat/prob/double_exponential_test.cpp new file mode 100644 index 00000000000..130931a60ce --- /dev/null +++ b/test/unit/math/prim/mat/prob/double_exponential_test.cpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct double_exponential_rng_wrapper { + template + auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + return stan::math::double_exponential_rng(mu, sigma, rng); + } +}; + +std::vector build_quantiles(int N, double mu, double sigma, + double unused) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + boost::math::laplace_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsDoubleExponential, errorCheck) { + check_dist_throws_all_types(double_exponential_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {}, {}); +} + +TEST(ProbDistributionsDoubleExponential, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, double_exponential_rng_wrapper{}, build_quantiles, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}, + {}); +} diff --git a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp new file mode 100644 index 00000000000..ec32a1d2422 --- /dev/null +++ b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp @@ -0,0 +1,18 @@ +#include +#include +#include + +struct exp_mod_normal_rng_wrapper { + template + auto operator()(const T1& loc, const T2& scale, const T3& inv_scale, + T_rng& rng) const { + return stan::math::exp_mod_normal_rng(loc, scale, inv_scale, rng); + } +}; + +TEST(ProbDistributionsExpModNormal, errorCheck) { + check_dist_throws_all_types(exp_mod_normal_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {0.2, 1.1, 2.5, 2.0}, {-3.7, -2.5, -1.0, 0.0}); +} diff --git a/test/unit/math/prim/mat/prob/gumbel_test.cpp b/test/unit/math/prim/mat/prob/gumbel_test.cpp new file mode 100644 index 00000000000..0ff6e44e38b --- /dev/null +++ b/test/unit/math/prim/mat/prob/gumbel_test.cpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct gumbel_rng_wrapper { + template + auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + return stan::math::gumbel_rng(mu, sigma, rng); + } +}; + +std::vector build_quantiles(int N, double mu, double sigma, + double unused) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + boost::math::extreme_value_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsGumbel, errorCheck) { + check_dist_throws_all_types(gumbel_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {}, {}); +} + +TEST(ProbDistributionsGumbel, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, gumbel_rng_wrapper{}, build_quantiles, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}, + {}); +} diff --git a/test/unit/math/prim/mat/prob/logistic_test.cpp b/test/unit/math/prim/mat/prob/logistic_test.cpp new file mode 100644 index 00000000000..1727b32c4d4 --- /dev/null +++ b/test/unit/math/prim/mat/prob/logistic_test.cpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct logistic_rng_wrapper { + template + auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + return stan::math::logistic_rng(mu, sigma, rng); + } +}; + +std::vector build_quantiles(int N, double mu, double sigma, + double unused) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + boost::math::logistic_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsLogistic, errorCheck) { + check_dist_throws_all_types(logistic_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {}, {}); +} + +TEST(ProbDistributionsLogistic, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, logistic_rng_wrapper{}, build_quantiles, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}, + {}); +} diff --git a/test/unit/math/prim/mat/prob/normal_test.cpp b/test/unit/math/prim/mat/prob/normal_test.cpp index b804239dd15..00b91a31b11 100644 --- a/test/unit/math/prim/mat/prob/normal_test.cpp +++ b/test/unit/math/prim/mat/prob/normal_test.cpp @@ -15,8 +15,8 @@ using Eigen::Matrix; * pass a templated function as an argument) */ struct normal_rng_wrapper { - template - auto operator()(const T1& mean, const T2& sd, T3& rng) const { + template + auto operator()(const T1& mean, const T2& sd, const T3& unused, T_rng& rng) const { return stan::math::normal_rng(mean, sd, rng); } }; @@ -25,7 +25,8 @@ struct normal_rng_wrapper { * This function builds the quantiles that we will supply to * assert_matches_quantiles to test the normal_rng */ -std::vector build_quantiles(int N, double mu, double sigma) { +std::vector build_quantiles(int N, double mu, double sigma, + double unused) { std::vector quantiles; int K = boost::math::round(2 * std::pow(N, 0.4)); boost::math::normal_distribution<> dist(mu, sigma); @@ -43,14 +44,18 @@ TEST(ProbDistributionsNormal, errorCheck) { /* * This test verifies that normal_rng throws errors in the right places. Test * inputs are picked from the four input initializer lists which are (in - * order): valid values for p1, invalid values for p1, - * valid values for p2, and invalid values for p2. + * order): + * valid values for p1, invalid values for p1, + * valid values for p2, and invalid values for p2, + * valid values for p3 (not used here), and invalid values for p3 (again + * not used). * * It does so for all possible combinations of calling arguments. */ check_dist_throws_all_types(normal_rng_wrapper{}, {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}); + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {}, {}); } TEST(ProbDistributionsNormal, chiSquareGoodnessFitTest) { @@ -60,12 +65,13 @@ TEST(ProbDistributionsNormal, chiSquareGoodnessFitTest) { /* * This test checks that the normal_rng is actually generating numbers from * the correct distributions. Test inputs are picked from the two input - * initializer lists which are (in order): valid values for p1, and valid - * values for p2. + * initializer lists which are (in order): valid values for p1, valid + * values for p2, and valid values for p3 (not used here). * * It does so for all possible combinations of calling arguments. */ check_quantiles_all_types(N, M, normal_rng_wrapper{}, build_quantiles, {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}); + {0.1, 1.0, 2.5, 4.0}, + {}); } diff --git a/test/unit/math/prim/mat/prob/skew_normal_test.cpp b/test/unit/math/prim/mat/prob/skew_normal_test.cpp new file mode 100644 index 00000000000..b57988637a9 --- /dev/null +++ b/test/unit/math/prim/mat/prob/skew_normal_test.cpp @@ -0,0 +1,51 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct skew_normal_rng_wrapper { + template + auto operator()(const T1& mu, const T2& sigma, const T3& alpha, + T_rng& rng) const { + return stan::math::skew_normal_rng(mu, sigma, alpha, rng); + } +}; + +std::vector build_quantiles(int N, double mu, double sigma, + double alpha) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + boost::math::skew_normal_distribution<> dist(mu, sigma, alpha); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsSkewNormal, errorCheck) { + check_dist_throws_all_types(skew_normal_rng_wrapper{}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}, {}); +} + +TEST(ProbDistributionsSkewNormal, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, skew_normal_rng_wrapper{}, build_quantiles, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}, + {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}); +} + diff --git a/test/unit/math/prim/mat/prob/student_t_test.cpp b/test/unit/math/prim/mat/prob/student_t_test.cpp new file mode 100644 index 00000000000..79e75199bfc --- /dev/null +++ b/test/unit/math/prim/mat/prob/student_t_test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include +#include +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +struct student_t_rng_wrapper { + template + auto operator()(const T1& nu, const T2& mu, const T3& sigma, T_rng& rng) const { + return stan::math::student_t_rng(nu, mu, sigma, rng); + } +}; + +std::vector build_quantiles(int N, double nu, double mu, double sigma) { + std::vector quantiles; + int K = boost::math::round(2 * std::pow(N, 0.4)); + + boost::math::students_t_distribution<>dist(nu); + for (int i=1; i(i) / K; + quantiles.push_back(quantile(dist, frac) * sigma + mu); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; +} + +TEST(ProbDistributionsStudentT, errorCheck) { + check_dist_throws_all_types(student_t_rng_wrapper{}, + {1.1, 2.0, 2.5, 2.0}, {-1.7, -0.5, -2.5, 0.0}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}); +} + +TEST(ProbDistributionsStudentT, chiSquareGoodnessFitTest) { + int N = 10000; + int M = 10; + + check_quantiles_all_types(N, M, student_t_rng_wrapper{}, build_quantiles, + {1.1, 2.0, 2.5, 2.0}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {0.1, 1.0, 2.5, 4.0}); +} diff --git a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp index 7f0d3a81478..4678a5d5934 100644 --- a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp +++ b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp @@ -21,6 +21,9 @@ template void assign_parameter_values(T_param& params, const std::vector& values) { + if(values.size() == 0) + return; + // Static cast the size here because sometimes it is unsigned (if T_param // is a std::vector) and sometimes it is signed (when it is an Eigen type) for (size_t i = 0; i < static_cast(params.size()); i++) { @@ -36,6 +39,9 @@ void assign_parameter_values(T_param& params, */ template <> void assign_parameter_values(double& param, const std::vector& values) { + if(values.size() == 0) + return; + param = values[0]; } @@ -61,9 +67,13 @@ void resize_if_vector(double& v, int N) { /* * check_dist_throws feeds the given generate_samples callable various * combinations of valid and invalid parameters (as defined by good_p1_, bad_p1, - * good_p2_, and bad_p2_). For all combinations of valid (good) parameters, - * generate_samples should throw no errors. For all combinations with an invalid - * (bad) parameter, generate_samples should throw domain_errors. + * good_p2_, bad_p2_, good_p3_, and bad_p3_). For all combinations of valid + * (good) parameters, generate_samples should throw no errors. For all + * combinations with an invalid (bad) parameter, generate_samples should throw + * domain_errors. + * + * If good_p2_ or good_p3_ are empty, then it is assumed that those parameters + * are unused and will not be tested. * * The generate_samples callable is also passed various other guaranteed invalid * values like positive infinity, negative infinity, and NaNs (these should also @@ -73,39 +83,55 @@ void resize_if_vector(double& v, int N) { * input vectors. These should cause invalid_argument errors. * * The generate_samples callable must have the signature: - * T_out generate_samples(T_param1 p1, T_param2 p2, T_rng) - * where T_param1 and T_param2 can be any combination of: + * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, T_rng) + * where T_param1, T_param2, and T_param3 can be any combination of: * double, std::vector, VectorXd, or RowVectorXd - * and T_rng is a boost random number generator + * and T_rng is a boost random number generator. p2 and p3 are required for the + * function signature, but they can be ignored if they are not needed. * * The output of generate_samples must be of - * stan::length(out) == stan::max_size(p1, p2) + * stan::length(out) == stan::length(p1) if only p1 is used + * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used + * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used * * Each element of the output, out[i], should correspond to a random variate - * generated according to parameters p1[i] and p2[i] + * generated according to parameters p1[i], p2[i], and p3[i] * * @tparam T_param1 Type of first parameter * @tparam T_param2 Type of second parameter + * @tparam T_param3 Type of third parameter * @tparam T_generate_samples Type of generate_samples functor * @param generate_samples generate_samples functor to test * @param good_p1_ Valid values of first parameter * @param bad_p1_ Invalid values of first parameter * @param good_p2_ Valid values of second parameter + * (leave empty if p2 is unused) * @param bad_p2_ Invalid values of second parameter + * @param good_p3_ Valid values of third parameter + * (leave empty if p3 is unused) + * @param bad_p3_ Invalid values of third parameter */ -template +template void check_dist_throws(T_generate_samples generate_samples, const std::vector& good_p1_, const std::vector& bad_p1_, const std::vector& good_p2_, - const std::vector& bad_p2_) { + const std::vector& bad_p2_, + const std::vector& good_p3_, + const std::vector& bad_p3_) { boost::random::mt19937 rng; T_param1 p1; T_param2 p2; + T_param2 p3; + + bool p2_is_used = good_p2_.size() > 0; + bool p3_is_used = good_p3_.size() > 0; resize_if_vector(p1, 5); // No-op if p1 is a scalar resize_if_vector(p2, 5); // No-op if p2 is a scalar + resize_if_vector(p3, 5); // No-op if p3 is a scalar // Make copies of the input arguments so that we can randomly shuffle them // in the tests @@ -113,32 +139,42 @@ void check_dist_throws(T_generate_samples generate_samples, std::vector bad_p1 = bad_p1_; std::vector good_p2 = good_p2_; std::vector bad_p2 = bad_p2_; + std::vector good_p3 = good_p3_; + std::vector bad_p3 = bad_p3_; // Try a few combinations of parameters that should work for (int i = 0; i < 5; i++) { std::random_shuffle(good_p1.begin(), good_p1.end()); std::random_shuffle(good_p2.begin(), good_p2.end()); + std::random_shuffle(good_p3.begin(), good_p3.end()); assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); - EXPECT_NO_THROW(generate_samples(p1, p2, rng)); + assign_parameter_values(p3, good_p3); + EXPECT_NO_THROW(generate_samples(p1, p2, p3, rng)); } // Now try putting incompatible values in first parameter - if (bad_p1.size() > 0) { - for (auto bad_p1_value : bad_p1) { - assign_parameter_values(p1, { bad_p1_value }); - assign_parameter_values(p2, good_p2); - EXPECT_THROW(generate_samples(p1, p2, rng), std::domain_error); - } + for (auto bad_p1_value : bad_p1) { + assign_parameter_values(p1, { bad_p1_value }); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); } // Now try putting incompatible values in second parameter - if (bad_p2.size() > 0) { - for (auto bad_p2_value : bad_p2) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_p2_value }); - EXPECT_THROW(generate_samples(p1, p2, rng), std::domain_error); - } + for (auto bad_p2_value : bad_p2) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, { bad_p2_value }); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + } + + // Now try putting incompatible values in third parameter + for (auto bad_p3_value : bad_p3) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, { bad_p3_value }); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); } // Make sure sampler throws errors with these values @@ -149,21 +185,61 @@ void check_dist_throws(T_generate_samples generate_samples, for (auto bad_value : bad_values) { assign_parameter_values(p1, { bad_value }); assign_parameter_values(p2, good_p2); - EXPECT_THROW(generate_samples(p1, p2, rng), std::domain_error); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + + if(p2_is_used) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, { bad_value }); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + } + + if(p3_is_used) { + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, { bad_value }); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + } + } + + // Check to make sure _rng rejects vector arguments of different lengths for + // all parameter pairs. + // If p1 is a scalar or the only vector, this test is skipped + resize_if_vector(p1, 3); // No-op if p1 is a scalar + resize_if_vector(p2, 4); // No-op if p2 is a scalar + resize_if_vector(p3, 4); // No-op if p3 is a scalar + if (stan::length(p1) != 1 && + ((p2_is_used && stan::length(p2) != 1) || + (p3_is_used && stan::length(p3) != 1))) { assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_value }); - EXPECT_THROW(generate_samples(p1, p2, rng), std::domain_error); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); } - // Check to make sure _rng rejects vector arguments of different lengths - // If either p1 or p2 are scalars, this test is skipped + // If p2 is a scalar or the only vector, this test is skipped resize_if_vector(p1, 4); // No-op if p1 is a scalar resize_if_vector(p2, 3); // No-op if p2 is a scalar - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - if (stan::length(p1) != 1 && stan::length(p2) != 1) - EXPECT_THROW(generate_samples(p1, p2, rng), std::invalid_argument); + resize_if_vector(p3, 4); // No-op if p3 is a scalar + if (p2_is_used && stan::length(p2) != 1 && + (stan::length(p1) != 1 || (p3_is_used && stan::length(p3) != 1))) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); + } + + // If p3 is a scalar or the only vector, this test is skipped + resize_if_vector(p1, 4); // No-op if p1 is a scalar + resize_if_vector(p2, 4); // No-op if p2 is a scalar + resize_if_vector(p3, 3); // No-op if p3 is a scalar + if (p3_is_used && stan::length(p3) != 1 && + (stan::length(p1) != 1 || (p2_is_used && stan::length(p2)))) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); + } } /* @@ -177,48 +253,151 @@ void check_dist_throws(T_generate_samples generate_samples, * @param bad_p1 Invalid values of first parameter * @param good_p2 Valid values of second parameter * @param bad_p2 Invalid values of second parameter + * @param good_p3 Valid values of third parameter + * @param bad_p3 Invalid values of third parameter */ template void check_dist_throws_all_types(T_generate_samples generate_samples, const std::vector& good_p1, const std::vector& bad_p1, const std::vector& good_p2, - const std::vector& bad_p2) { + const std::vector& bad_p2, + const std::vector& good_p3, + const std::vector& bad_p3) { using Eigen::VectorXd; using Eigen::RowVectorXd; - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws, double> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double, double> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector, double > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd, double> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd, double> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws > + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, double, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, std::vector, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, VectorXd, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws, RowVectorXd> + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + check_dist_throws + (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); } /* @@ -244,71 +423,94 @@ std::vector promote_to_vector(std::vector v) { * generated by the callable generate_samples against the reference * quantiles computed with the given callable generate_quantiles. * + * If good_p2 or good_p3 are empty, then it is assumed that those parameters + * are unused and will not be tested. + * * The generate_samples callable must have the signature: - * T_out generate_samples(T_param1 p1, T_param2 p2, T_rng) - * where T_param1 and T_param2 can be any combination of: + * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, p3, T_rng) + * where T_param1, T_param2, and T_param3 can be any combination of: * double, std::vector, VectorXd, or RowVectorXd - * and T_rng is a boost random number generator + * and T_rng is a boost random number generator. p2 and p3 are required for the + * function signature, but they can be ignored if they are not needed. * * The output of generate_samples must be of - * stan::length(out) == stan::max_size(p1, p2) + * stan::length(out) == stan::length(p1) if only p1 is used + * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used + * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used * * Each element of the output, out[i], should correspond to a random variate - * generated according to parameters p1[i] and p2[i] + * generated according to parameters p1[i], p2[i], and p3[i] * * The generate_quantiles callable must have the signature: - * std::vector generate_samples(int N, double p1, double p2) + * std::vector generate_samples(int N, double p1, double p2, double p3) + * (similarly to above, ignore parameters p2 and p3 if they aren't necessary) * * N is the number of samples that the quantiles will be compared - * against and p1 and p2 the values of the parameters that these quantiles - * correspond to. + * against and p1, p2, and p3 are the values of the parameters that these + * quantiles correspond to. * * The output of generate_quantiles should be the quantiles used for comparing * empirical distributions used by assert_matches_quantiles * - * Parameters for each of these callables are copied from the vectors good_p1 - * and good_p2. + * Parameters for each of these callables are copied from the vectors good_p1, + * good_p2, and good_p3. * * @tparam T_param1 Type of first parameter * @tparam T_param2 Type of second parameter + * @tparam T_param3 Type of third parameter * @tparam T_generate_samples Type of callable for generating samples * @tparam T_generate_quantiles Type of callable for generating quantiles * @param sampler sampler functor to test * @param good_p1 Valid values of first parameter * @param good_p2 Valid values of second parameter + * @param good_p3 Valid values of third parameter */ -template void check_quantiles(int N, int M_vec, T_generate_samples generate_samples, T_generate_quantiles generate_quantiles, const std::vector& good_p1, - const std::vector& good_p2) { + const std::vector& good_p2, + const std::vector& good_p3) { boost::random::mt19937 rng; T_param1 p1; T_param2 p2; + T_param3 p3; resize_if_vector(p1, M_vec); // No-op if p1 is scalar resize_if_vector(p2, M_vec); // No-op if p2 is scalar + resize_if_vector(p3, M_vec); // No-op if p3 is scalar + assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); - int M = stan::max_size(p1, p2); + assign_parameter_values(p3, good_p3); + + bool p2_is_used = good_p2.size() > 0; + bool p3_is_used = good_p3.size() > 0; + + int M = std::max({ stan::length(p1), + (p2_is_used) ? stan::length(p2) : 1, + (p3_is_used) ? stan::length(p3) : 1 }); stan::scalar_seq_view p1_vec(p1); stan::scalar_seq_view p2_vec(p2); + stan::scalar_seq_view p3_vec(p3); std::vector > samples_to_test_transpose; for (int n = 0; n < N; ++n) { - // If p1 and p2 are scalars, the output is a scalar. Need to promote it + // If p1, p2, and p3 are scalars, the output is a scalar. Need to promote it // to a std::vector samples_to_test_transpose. - push_back(promote_to_vector(generate_samples(p1, p2, rng))); + push_back(promote_to_vector(generate_samples(p1, p2, p3, rng))); } + for (int m = 0; m < M; ++m) { std::vector samples_to_test; for (int n = 0; n < N; ++n) { samples_to_test.push_back(samples_to_test_transpose[n][m]); } - std::vector quantiles = generate_quantiles(N, p1_vec[m], p2_vec[m]); + std::vector quantiles = + generate_quantiles(N, p1_vec[m], p2_vec[m], p3_vec[m]); assert_matches_quantiles(samples_to_test, quantiles, 1e-6); } @@ -329,42 +531,142 @@ void check_quantiles_all_types(int N, int M, T_generate_samples generate_samples, T_generate_quantiles generate_quantiles, const std::vector& good_p1, - const std::vector& good_p2) { + const std::vector& good_p2, + const std::vector& good_p3) { using Eigen::VectorXd; using Eigen::RowVectorXd; - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles, VectorXd > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles, RowVectorXd > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles > + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, double, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, std::vector, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, VectorXd, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles, RowVectorXd> + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + check_quantiles + (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); } #endif diff --git a/test/unit/math/prim/scal/prob/util.hpp b/test/unit/math/prim/scal/prob/util.hpp index bcba8137597..b1141686801 100644 --- a/test/unit/math/prim/scal/prob/util.hpp +++ b/test/unit/math/prim/scal/prob/util.hpp @@ -10,18 +10,18 @@ */ void assert_chi_squared(const std::vector& counts, const std::vector& expected, double tolerance) { int bins = counts.size(); - EXPECT_EQ(bins, expected.size()); + EXPECT_EQ(bins, expected.size()); double chi = 0; - for (int i=0; i& counts, const std::vector& samples, const std::vector& quantiles, double tolerance) { int N = samples.size(); - std::vector mysamples = samples; - std::sort(mysamples.begin(), mysamples.end()); + std::vector mysamples = samples; + std::sort(mysamples.begin(), mysamples.end()); - int K = quantiles.size(); - double expected_count = static_cast(N) / K; + int K = quantiles.size(); + double expected_count = static_cast(N) / K; std::vector expected; - for (int i=0; i counts(K); + std::vector counts(K); int current_index = 0; - for (int i=0; i= quantiles[current_index]) { ++current_index; + EXPECT_TRUE(current_index < quantiles.size()); } - ++counts[current_index]; + ++counts[current_index]; } - assert_chi_squared(counts, expected, tolerance); + assert_chi_squared(counts, expected, tolerance); } - #endif From fa7187b1d1db94b130a75434dbc579f2947cd320 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 18 Sep 2017 12:22:26 -0700 Subject: [PATCH 02/12] Changed int to size_t to avoid signed comparison warning (Issue 45) --- test/unit/math/prim/scal/prob/util.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/scal/prob/util.hpp b/test/unit/math/prim/scal/prob/util.hpp index b1141686801..e008b10fca2 100644 --- a/test/unit/math/prim/scal/prob/util.hpp +++ b/test/unit/math/prim/scal/prob/util.hpp @@ -43,7 +43,7 @@ void assert_matches_quantiles(const std::vector& samples, const std::vec } std::vector counts(K); - int current_index = 0; + size_t current_index = 0; for (int i=0; i= quantiles[current_index]) { ++current_index; From ba6c624202dba113b365ed1ad07e2708d80a3245 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 30 Sep 2017 08:55:14 -0700 Subject: [PATCH 03/12] Moved around & symbols and if statements. Avoided cast w/ templating --- stan/math/prim/scal/prob/cauchy_rng.hpp | 4 +-- .../prim/scal/prob/double_exponential_rng.hpp | 4 +-- stan/math/prim/scal/prob/gumbel_rng.hpp | 4 +-- stan/math/prim/scal/prob/logistic_rng.hpp | 4 +-- stan/math/prim/scal/prob/normal_rng.hpp | 4 +-- .../prim/mat/prob/vector_rng_test_helper.hpp | 32 +++++++++++++++---- 6 files changed, 35 insertions(+), 17 deletions(-) diff --git a/stan/math/prim/scal/prob/cauchy_rng.hpp b/stan/math/prim/scal/prob/cauchy_rng.hpp index 81bca55b142..fe31a7fe065 100644 --- a/stan/math/prim/scal/prob/cauchy_rng.hpp +++ b/stan/math/prim/scal/prob/cauchy_rng.hpp @@ -32,8 +32,8 @@ namespace stan { */ template inline typename VectorBuilder::type - cauchy_rng(const T_loc &mu, - const T_scale &sigma, + cauchy_rng(const T_loc& mu, + const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::cauchy_distribution; diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 9fbf182d6ba..fdf864bea2b 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -33,8 +33,8 @@ namespace stan { */ template inline typename VectorBuilder::type - double_exponential_rng(const T_loc &mu, - const T_scale &sigma, + double_exponential_rng(const T_loc& mu, + const T_scale& sigma, RNG& rng) { static const std::string function = "double_exponential_rng"; diff --git a/stan/math/prim/scal/prob/gumbel_rng.hpp b/stan/math/prim/scal/prob/gumbel_rng.hpp index 4fac82e770d..f5dc42985d3 100644 --- a/stan/math/prim/scal/prob/gumbel_rng.hpp +++ b/stan/math/prim/scal/prob/gumbel_rng.hpp @@ -32,8 +32,8 @@ namespace stan { */ template inline typename VectorBuilder::type - gumbel_rng(const T_loc &mu, - const T_scale &beta, + gumbel_rng(const T_loc& mu, + const T_scale& beta, RNG& rng) { using boost::variate_generator; using boost::uniform_01; diff --git a/stan/math/prim/scal/prob/logistic_rng.hpp b/stan/math/prim/scal/prob/logistic_rng.hpp index da1fb64012a..bef1e41e5df 100644 --- a/stan/math/prim/scal/prob/logistic_rng.hpp +++ b/stan/math/prim/scal/prob/logistic_rng.hpp @@ -32,8 +32,8 @@ namespace stan { */ template inline typename VectorBuilder::type - logistic_rng(const T_loc &mu, - const T_scale &sigma, + logistic_rng(const T_loc& mu, + const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::exponential_distribution; diff --git a/stan/math/prim/scal/prob/normal_rng.hpp b/stan/math/prim/scal/prob/normal_rng.hpp index 84d39405c5d..bf054f1459d 100644 --- a/stan/math/prim/scal/prob/normal_rng.hpp +++ b/stan/math/prim/scal/prob/normal_rng.hpp @@ -32,8 +32,8 @@ namespace stan { */ template inline typename VectorBuilder::type - normal_rng(const T_loc &mu, - const T_scale &sigma, + normal_rng(const T_loc& mu, + const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::normal_distribution; diff --git a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp index 4678a5d5934..39e81e5444e 100644 --- a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp +++ b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp @@ -21,12 +21,30 @@ template void assign_parameter_values(T_param& params, const std::vector& values) { - if(values.size() == 0) + if (values.size() == 0) return; - // Static cast the size here because sometimes it is unsigned (if T_param - // is a std::vector) and sometimes it is signed (when it is an Eigen type) - for (size_t i = 0; i < static_cast(params.size()); i++) { + for (int i = 0; i < params.size(); i++) { + params[i] = values[i % values.size()]; + } +} + +/* + * Fill the vector-like variable params with values from the values argument. + * + * values can be shorter than params (in which cause multiple copies of values + * are tiled over the params vector) + * + * @param params Parameter vector to write values to + * @param params Values to copy into params + */ +template <> +void assign_parameter_values(std::vector& params, + const std::vector& values) { + if (values.size() == 0) + return; + + for (size_t i = 0; i < params.size(); i++) { params[i] = values[i % values.size()]; } } @@ -39,7 +57,7 @@ void assign_parameter_values(T_param& params, */ template <> void assign_parameter_values(double& param, const std::vector& values) { - if(values.size() == 0) + if (values.size() == 0) return; param = values[0]; @@ -188,13 +206,13 @@ void check_dist_throws(T_generate_samples generate_samples, assign_parameter_values(p3, good_p3); EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); - if(p2_is_used) { + if (p2_is_used) { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, { bad_value }); EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); } - if(p3_is_used) { + if (p3_is_used) { assign_parameter_values(p2, good_p2); assign_parameter_values(p3, { bad_value }); EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); From 275889c5ba9185ea8cf6ec70b516013e4f63cb6f Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sun, 1 Oct 2017 09:38:26 -0700 Subject: [PATCH 04/12] This commit addresses Bob's comments (from like September 30 or something) 1. Switched to typelists for generating all combinations of tests 2. Random number generators are now tested by building test rigs (that inherit from an abstract base class) and passing those to test functions 3. Wasn't able to easily add tests for ints This is all Issue 45 --- .../prim/scal/prob/double_exponential_rng.hpp | 8 +- .../prim/scal/prob/exp_mod_normal_rng.hpp | 1 - stan/math/prim/scal/prob/gumbel_rng.hpp | 3 +- stan/math/prim/scal/prob/logistic_rng.hpp | 4 +- stan/math/prim/scal/prob/skew_normal_rng.hpp | 13 +- test/unit/math/prim/mat/prob/cauchy_test.cpp | 55 +- .../prim/mat/prob/double_exponential_test.cpp | 55 +- .../prim/mat/prob/exp_mod_normal_test.cpp | 23 +- test/unit/math/prim/mat/prob/gumbel_test.cpp | 55 +- .../unit/math/prim/mat/prob/logistic_test.cpp | 55 +- test/unit/math/prim/mat/prob/normal_test.cpp | 93 +- .../math/prim/mat/prob/skew_normal_test.cpp | 56 +- .../math/prim/mat/prob/student_t_test.cpp | 55 +- .../prim/mat/prob/vector_rng_test_helper.hpp | 954 ++++++++---------- 14 files changed, 655 insertions(+), 775 deletions(-) diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index fdf864bea2b..32b4d72a447 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -55,14 +55,10 @@ namespace stan { VectorBuilder output(N); + variate_generator > rng_unit_01(rng, uniform_01<>()); for (size_t n = 0; n < N; n++) { - variate_generator > rng_unit_01(rng, uniform_01<>()); - double a = 0; double laplaceRN = rng_unit_01(); - if (0.5 - laplaceRN > 0) - a = 1.0; - else if (0.5 - laplaceRN < 0) - a = -1.0; + double a = (0.5 - laplaceRN > 0) ? 1.0 : -1.0; output[n] = mu_vec[n] - sigma_vec[n] * a * log1m(2 * abs(0.5 - laplaceRN)); } diff --git a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp index 4d52d8f1115..99be01bd1d0 100644 --- a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp +++ b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp @@ -67,7 +67,6 @@ namespace stan { return output.data(); } - } } #endif diff --git a/stan/math/prim/scal/prob/gumbel_rng.hpp b/stan/math/prim/scal/prob/gumbel_rng.hpp index f5dc42985d3..1e1596ca7c5 100644 --- a/stan/math/prim/scal/prob/gumbel_rng.hpp +++ b/stan/math/prim/scal/prob/gumbel_rng.hpp @@ -53,9 +53,8 @@ namespace stan { VectorBuilder output(N); + variate_generator > uniform01_rng(rng, uniform_01<>()); for (size_t n = 0; n < N; n++) { - variate_generator > - uniform01_rng(rng, uniform_01<>()); output[n] = mu_vec[n] - beta_vec[n] * std::log(-std::log(uniform01_rng())); } diff --git a/stan/math/prim/scal/prob/logistic_rng.hpp b/stan/math/prim/scal/prob/logistic_rng.hpp index bef1e41e5df..7273a48d49d 100644 --- a/stan/math/prim/scal/prob/logistic_rng.hpp +++ b/stan/math/prim/scal/prob/logistic_rng.hpp @@ -53,9 +53,9 @@ namespace stan { VectorBuilder output(N); + variate_generator > + exp_rng(rng, exponential_distribution<>(1)); for (size_t n = 0; n < N; n++) { - variate_generator > - exp_rng(rng, exponential_distribution<>(1)); output[n] = mu_vec[n] - sigma_vec[n] * std::log(exp_rng() / exp_rng()); } diff --git a/stan/math/prim/scal/prob/skew_normal_rng.hpp b/stan/math/prim/scal/prob/skew_normal_rng.hpp index 2665a093192..a2c24bf5504 100644 --- a/stan/math/prim/scal/prob/skew_normal_rng.hpp +++ b/stan/math/prim/scal/prob/skew_normal_rng.hpp @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include namespace stan { @@ -37,6 +39,9 @@ namespace stan { const T_scale& sigma, const T_shape& alpha, RNG& rng) { + using boost::variate_generator; + using boost::random::normal_distribution; + static const std::string function = "skew_normal_rng"; check_finite(function, "Location parameter", mu); @@ -55,12 +60,14 @@ namespace stan { VectorBuilder output(N); + variate_generator > + norm_rng(rng, normal_distribution<>(0.0, 1.0)); for (size_t n = 0; n < N; n++) { - double r1 = normal_rng(0.0, 1.0, rng), - r2 = normal_rng(0.0, 1.0, rng); + double r1 = norm_rng(), + r2 = norm_rng(); if (r2 > alpha_vec[n] * r1) - r1 *= -1; + r1 = -r1; output[n] = mu_vec[n] + sigma_vec[n] * r1; } diff --git a/test/unit/math/prim/mat/prob/cauchy_test.cpp b/test/unit/math/prim/mat/prob/cauchy_test.cpp index bd43b9a23e3..35f5dabc7d8 100644 --- a/test/unit/math/prim/mat/prob/cauchy_test.cpp +++ b/test/unit/math/prim/mat/prob/cauchy_test.cpp @@ -6,44 +6,39 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; - -struct cauchy_rng_wrapper { +class CauchyTestRig : public VectorRNGTestRig { +public: + CauchyTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + template - auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, + T_rng& rng) const { return stan::math::cauchy_rng(mu, sigma, rng); } -}; -std::vector build_quantiles(int N, double mu, double sigma, - double unused) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::cauchy_distribution<> dist(mu, sigma); - - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + std::vector generate_quantiles(double mu, double sigma, double unused) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::cauchy_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsCauchy, errorCheck) { - check_dist_throws_all_types(cauchy_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {}, {}); + check_dist_throws_all_types(CauchyTestRig()); } TEST(ProbDistributionsCauchy, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, cauchy_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {}); + check_quantiles_all_types(CauchyTestRig()); } diff --git a/test/unit/math/prim/mat/prob/double_exponential_test.cpp b/test/unit/math/prim/mat/prob/double_exponential_test.cpp index 130931a60ce..ba1f76ca576 100644 --- a/test/unit/math/prim/mat/prob/double_exponential_test.cpp +++ b/test/unit/math/prim/mat/prob/double_exponential_test.cpp @@ -6,44 +6,39 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; +class DoubleExponentialTestRig : public VectorRNGTestRig { +public: + DoubleExponentialTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} -struct double_exponential_rng_wrapper { template - auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, + T_rng& rng) const { return stan::math::double_exponential_rng(mu, sigma, rng); } -}; - -std::vector build_quantiles(int N, double mu, double sigma, - double unused) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::laplace_distribution<> dist(mu, sigma); - - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + + std::vector generate_quantiles(double mu, double sigma, double unused) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::laplace_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsDoubleExponential, errorCheck) { - check_dist_throws_all_types(double_exponential_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {}, {}); + check_dist_throws_all_types(DoubleExponentialTestRig()); } TEST(ProbDistributionsDoubleExponential, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, double_exponential_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {}); + check_quantiles_all_types(DoubleExponentialTestRig()); } diff --git a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp index ec32a1d2422..99752ae0e22 100644 --- a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp @@ -2,17 +2,26 @@ #include #include -struct exp_mod_normal_rng_wrapper { +class ExpModNormalTestRig : public VectorRNGTestRig { +public: + ExpModNormalTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {0.2, 1.1, 2.5, 2.0}, {-3.7, -2.5, -1.0, 0.0}) {} + template - auto operator()(const T1& loc, const T2& scale, const T3& inv_scale, - T_rng& rng) const { + auto generate_samples(const T1& loc, const T2& scale, const T3& inv_scale, + T_rng& rng) const { return stan::math::exp_mod_normal_rng(loc, scale, inv_scale, rng); } + + std::vector + generate_quantiles(double loc, double scale, double inv_scale) const { + return std::vector(); + } }; TEST(ProbDistributionsExpModNormal, errorCheck) { - check_dist_throws_all_types(exp_mod_normal_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {0.2, 1.1, 2.5, 2.0}, {-3.7, -2.5, -1.0, 0.0}); + check_dist_throws_all_types(ExpModNormalTestRig()); } diff --git a/test/unit/math/prim/mat/prob/gumbel_test.cpp b/test/unit/math/prim/mat/prob/gumbel_test.cpp index 0ff6e44e38b..223466de8b4 100644 --- a/test/unit/math/prim/mat/prob/gumbel_test.cpp +++ b/test/unit/math/prim/mat/prob/gumbel_test.cpp @@ -6,44 +6,39 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; - -struct gumbel_rng_wrapper { +class GumbelTestRig : public VectorRNGTestRig { +public: + GumbelTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-1.0, -1.5, -2.5, -0.7, 0.0}) {} + template - auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, + T_rng& rng) const { return stan::math::gumbel_rng(mu, sigma, rng); } -}; -std::vector build_quantiles(int N, double mu, double sigma, - double unused) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::extreme_value_distribution<> dist(mu, sigma); - - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + std::vector generate_quantiles(double mu, double sigma, double unused) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::extreme_value_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsGumbel, errorCheck) { - check_dist_throws_all_types(gumbel_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {}, {}); + check_dist_throws_all_types(GumbelTestRig()); } TEST(ProbDistributionsGumbel, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, gumbel_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {}); + check_quantiles_all_types(GumbelTestRig()); } diff --git a/test/unit/math/prim/mat/prob/logistic_test.cpp b/test/unit/math/prim/mat/prob/logistic_test.cpp index 1727b32c4d4..fb211e11183 100644 --- a/test/unit/math/prim/mat/prob/logistic_test.cpp +++ b/test/unit/math/prim/mat/prob/logistic_test.cpp @@ -6,44 +6,39 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; - -struct logistic_rng_wrapper { +class LogisticTestRig : public VectorRNGTestRig { +public: + LogisticTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + template - auto operator()(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, + T_rng& rng) const { return stan::math::logistic_rng(mu, sigma, rng); } -}; -std::vector build_quantiles(int N, double mu, double sigma, - double unused) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::logistic_distribution<> dist(mu, sigma); - - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + std::vector generate_quantiles(double mu, double sigma, double unused) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::logistic_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsLogistic, errorCheck) { - check_dist_throws_all_types(logistic_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {}, {}); + check_dist_throws_all_types(LogisticTestRig()); } TEST(ProbDistributionsLogistic, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, logistic_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {}); + check_quantiles_all_types(LogisticTestRig()); } diff --git a/test/unit/math/prim/mat/prob/normal_test.cpp b/test/unit/math/prim/mat/prob/normal_test.cpp index 00b91a31b11..cc7c93191ac 100644 --- a/test/unit/math/prim/mat/prob/normal_test.cpp +++ b/test/unit/math/prim/mat/prob/normal_test.cpp @@ -9,69 +9,68 @@ using Eigen::Dynamic; using Eigen::Matrix; -/* - * This needs to be a functor instead of a function because we need to pass it - * as an argument to another function (and I'm not sure it's possible to easily - * pass a templated function as an argument) - */ -struct normal_rng_wrapper { +class NormalTestRig : public VectorRNGTestRig { +public: + /* + * The default NormalTestRig constructor initializes the TestRig with + * valid and invalid parameters for a random number generator with two + * arguments. + */ + NormalTestRig() : + VectorRNGTestRig(10000, // Number of samples used for quantiles tests + 10, // Length of vectors for vectorization tests + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + + /* + * This function wraps up the random number generator for testing. + * + * The tested rng can have up to three parameters. Any unused parameters can + * be ignored. + */ template - auto operator()(const T1& mean, const T2& sd, const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& mean, const T2& sd, const T3& unused, + T_rng& rng) const { return stan::math::normal_rng(mean, sd, rng); } -}; - -/* - * This function builds the quantiles that we will supply to - * assert_matches_quantiles to test the normal_rng - */ -std::vector build_quantiles(int N, double mu, double sigma, - double unused) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::normal_distribution<> dist(mu, sigma); - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + /* + * This function builds the quantiles that we will supply to + * assert_matches_quantiles to test the normal_rng + */ + std::vector generate_quantiles(double mu, double sigma, double unused) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::normal_distribution<> dist(mu, sigma); + + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsNormal, errorCheck) { /* - * This test verifies that normal_rng throws errors in the right places. Test - * inputs are picked from the four input initializer lists which are (in - * order): - * valid values for p1, invalid values for p1, - * valid values for p2, and invalid values for p2, - * valid values for p3 (not used here), and invalid values for p3 (again - * not used). + * This test verifies that normal_rng throws errors in the right places. * - * It does so for all possible combinations of calling arguments. + * It does so by calling test_rig::generate_samples for all possible + * combinations of calling arguments. */ - check_dist_throws_all_types(normal_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {}, {}); + check_dist_throws_all_types(NormalTestRig()); } TEST(ProbDistributionsNormal, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - /* * This test checks that the normal_rng is actually generating numbers from - * the correct distributions. Test inputs are picked from the two input - * initializer lists which are (in order): valid values for p1, valid - * values for p2, and valid values for p3 (not used here). + * the correct distributions. Quantiles are computed from + * test_rig::generate_quantiles * * It does so for all possible combinations of calling arguments. */ - check_quantiles_all_types(N, M, normal_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {}); + check_quantiles_all_types(NormalTestRig()); } diff --git a/test/unit/math/prim/mat/prob/skew_normal_test.cpp b/test/unit/math/prim/mat/prob/skew_normal_test.cpp index b57988637a9..dd2c448f622 100644 --- a/test/unit/math/prim/mat/prob/skew_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/skew_normal_test.cpp @@ -6,46 +6,42 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; -struct skew_normal_rng_wrapper { +class SkewNormalTestRig : public VectorRNGTestRig { +public: + SkewNormalTestRig() : + VectorRNGTestRig(10000, 10, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, + {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}, {}) {} + template - auto operator()(const T1& mu, const T2& sigma, const T3& alpha, - T_rng& rng) const { + auto generate_samples(const T1& mu, const T2& sigma, const T3& alpha, + T_rng& rng) const { return stan::math::skew_normal_rng(mu, sigma, alpha, rng); } -}; - -std::vector build_quantiles(int N, double mu, double sigma, - double alpha) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - boost::math::skew_normal_distribution<> dist(mu, sigma, alpha); - for (int i = 1; i < K; ++i) { - double frac = static_cast(i) / K; - quantiles.push_back(quantile(dist, frac)); + std::vector generate_quantiles(double mu, double sigma, double alpha) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + boost::math::skew_normal_distribution<> dist(mu, sigma, alpha); + + for (int i = 1; i < K; ++i) { + double frac = i / K; + quantiles.push_back(quantile(dist, frac)); + } + quantiles.push_back(std::numeric_limits::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsSkewNormal, errorCheck) { - check_dist_throws_all_types(skew_normal_rng_wrapper{}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}, {}); + check_dist_throws_all_types(SkewNormalTestRig()); } TEST(ProbDistributionsSkewNormal, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, skew_normal_rng_wrapper{}, build_quantiles, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}, - {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}); + check_quantiles_all_types(SkewNormalTestRig()); } diff --git a/test/unit/math/prim/mat/prob/student_t_test.cpp b/test/unit/math/prim/mat/prob/student_t_test.cpp index 79e75199bfc..7e10a7a9e17 100644 --- a/test/unit/math/prim/mat/prob/student_t_test.cpp +++ b/test/unit/math/prim/mat/prob/student_t_test.cpp @@ -6,43 +6,40 @@ #include #include -using Eigen::Dynamic; -using Eigen::Matrix; - -struct student_t_rng_wrapper { +class StudentTTestRig : public VectorRNGTestRig { +public: + StudentTTestRig() : + VectorRNGTestRig(10000, 10, + {1.1, 2.0, 2.5, 2.0}, {-1.7, -0.5, -2.5, 0.0}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, + {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + template - auto operator()(const T1& nu, const T2& mu, const T3& sigma, T_rng& rng) const { + auto generate_samples(const T1& nu, const T2& mu, const T3& sigma, + T_rng& rng) const { return stan::math::student_t_rng(nu, mu, sigma, rng); } -}; -std::vector build_quantiles(int N, double nu, double mu, double sigma) { - std::vector quantiles; - int K = boost::math::round(2 * std::pow(N, 0.4)); - - boost::math::students_t_distribution<>dist(nu); - for (int i=1; i(i) / K; - quantiles.push_back(quantile(dist, frac) * sigma + mu); + std::vector generate_quantiles(double nu, double mu, double sigma) + const { + std::vector quantiles; + double K = boost::math::round(2 * std::pow(N_, 0.4)); + + boost::math::students_t_distribution<>dist(nu); + for (int i=1; i::max()); + + return quantiles; } - quantiles.push_back(std::numeric_limits::max()); - - return quantiles; -} +}; TEST(ProbDistributionsStudentT, errorCheck) { - check_dist_throws_all_types(student_t_rng_wrapper{}, - {1.1, 2.0, 2.5, 2.0}, {-1.7, -0.5, -2.5, 0.0}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}); + check_dist_throws_all_types(StudentTTestRig()); } TEST(ProbDistributionsStudentT, chiSquareGoodnessFitTest) { - int N = 10000; - int M = 10; - - check_quantiles_all_types(N, M, student_t_rng_wrapper{}, build_quantiles, - {1.1, 2.0, 2.5, 2.0}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, - {0.1, 1.0, 2.5, 4.0}); + check_quantiles_all_types(StudentTTestRig()); } diff --git a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp index 39e81e5444e..faea885789b 100644 --- a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp +++ b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp @@ -8,6 +8,246 @@ #include #include +/* + * VectorRNGTestRig is an Abstract class for wrapping up random number + * generators in a way so that they can be tested to + * 1) provide interfaces to all the necessary types + * 2) generate the correct distributions + * + * A test rig that inherits from VectorRNGTestRig must implement + * 1) generate_samples + * 2) generate_quantiles + * and must initialize good_p1_, bad_p1_, good_p2_, bad_p2, good_p3_, bad_p3_ + * (depending on how many parameters the wrapped random number generator uses) + * + * generate_samples is a wrapper around the actual random number generator. It + * is expected to have the signature: + * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, T_rng) + * where T_param1, T_param2, and T_param3 can be any combination of: + * double, std::vector, VectorXd, or RowVectorXd + * and T_rng is a boost random number generator. p2 and p3 are required for the + * function signature, but they can be ignored if they are not needed. + * + * Each element of the output, out[i], should correspond to a random variate + * generated according to parameters p1[i], p2[i], and p3[i] + * + * generate_samples is expected to throw domain_errors when passed invalid + * values like positive infinity, negative infinity, and NaNs. + * + * generated_samples is also expected to reject incompatibly sized + * input vectors. These should cause invalid_argument errors. + * + * The output of generate_samples must be of + * stan::length(out) == stan::length(p1) if only p1 is used + * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used + * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used + * It must be defined to take three arguments, but not all must be used. + * + * The generate_quantiles callable must have the signature: + * std::vector generate_samples(double p1, double p2, double p3) + * (similarly to above, ignore parameters p2 and p3 if they aren't necessary) + * + * generate_quantiles should compute the quantiles that will be compared against + * (using assert_matches_quantiles) empirical distributions generated with + * generate_samples (which will be of length N_ for each parameter combination). + * + * p1, p2, and p3 are the values of the parameters that these quantiles + * correspond to. + * + * good_p1_ and bad_p1_ should be initialized to lists of valid and invalid + * parameters for the first parameter of the tested RNG + * + * Likewise, good_p2_, bad_p2_, and good_p3_, bad_p3_ should be initialized. + * + * If testing a distribution that uses only two parameters, leave good_p3_ as + * a length zero vector and it will not be tested. Similarly, if good_p2_ is + * length zero, the second parameter will not be tested either + */ +class VectorRNGTestRig { +public: + int N_; // Number of samples used in the quantiles tests + int M_; // Length of vectors for the vectorization tests + + std::vector good_p1_; + std::vector bad_p1_; + std::vector good_p2_; + std::vector bad_p2_; + std::vector good_p3_; + std::vector bad_p3_; + + /* + * This function wraps up the random number generator for testing. + * + * The tested rng can have up to three parameters. Any unused parameters can + * be ignored. + * + * It *must* be implemented in the child class (it isn't virtual here because + * C++ doesn't allow templated virtual member functions) + */ + template + auto generate_samples(const T1& p1, const T2& p2, const T3& p3, + T_rng& rng) const; + + /* + * This function builds the quantiles that we will supply to + * assert_matches_quantiles to test the random number generator implemented in + * generate_samples + */ + virtual std::vector generate_quantiles(double p1, double p2, + double p3) const = 0; + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector bad_p1, + std::vector good_p2, + std::vector bad_p2, + std::vector good_p3, + std::vector bad_p3) : N_(N), M_(M), + good_p1_(good_p1), + bad_p1_(bad_p1), + good_p2_(good_p2), + bad_p2_(bad_p2), + good_p3_(good_p3), + bad_p3_(bad_p3) { + } + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector bad_p1, + std::vector good_p2, + std::vector bad_p2) : N_(N), M_(M), + good_p1_(good_p1), + bad_p1_(bad_p1), + good_p2_(good_p2), + bad_p2_(bad_p2){ + } + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector bad_p1) : N_(N), M_(M), + good_p1_(good_p1), + bad_p1_(bad_p1) { + } +}; + +/* + * call_permutations_helper is the primary function for the recursive template + * function call_permutations. It calls func(rig), + * and continues the recursion by decrementing K. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_rig Random number generator test rig + * tparam I Loop index for first type + * tparam J Loop index for second type + * tparam K Loop index for third type + */ +template +struct call_permutations_helper { + void operator()(const T_functor& func, const T_rig& rig) const { + func.template operator()::type, + typename std::tuple_element::type, + typename std::tuple_element::type> + (rig); + + call_permutations_helper{}(func, + rig); + } +}; + +/* + * This edge case catches when the right-most argument has completed one + * iteration through all types. It computes func, + * carries from the second argument, and continues the recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_rig Random number generator test rig + * tparam I Loop index for first type + * tparam J Loop index for second type + * tparam K Loop index for third type + */ +template +struct call_permutations_helper { + void operator()(const T_functor& func, const T_rig& rig) const { + func.template operator()::type, + typename std::tuple_element::type, + typename std::tuple_element<0, T_typelist>::type> + (rig); + + call_permutations_helper::value - 1>{}(func, rig); + } +}; + +/* + * This edge case catches when the second argument has completed one iteration + * through all types. It computes func, carries from + * the left-most argument and continues the recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_rig Random number generator test rig + * tparam I Loop index for first type + */ +template +struct call_permutations_helper { + void operator()(const T_functor& func, const T_rig& rig) const { + func.template operator()::type, + typename std::tuple_element<0, T_typelist>::type, + typename std::tuple_element<0, T_typelist>::type> + (rig); + + call_permutations_helper::value - 1, + std::tuple_size::value - 1>{}(func, + rig); + } +}; + +/* + * This edge case catches when all types have been iterated through for each + * template argument. It computes func and ends the + * recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_rig Random number generator test rig + */ +template +struct call_permutations_helper { + void operator()(const T_functor& func, const T_rig& rig) const { + func.template operator()::type, + typename std::tuple_element<0, T_typelist>::type, + typename std::tuple_element<0, T_typelist>::type> + (rig); + } +}; + + +/* + * call_permutations uses template recursion to call the functor func + * with all 3-element permutations of the types in T_typelist + * + * Roughly, this corresponds to: + * for (i in (I - 1):0) // indexes go backwards, like: I - 1, I - 2, ... 0 + * for (j in (J - 1):0) + * for (k in (K - 1):0) + * func(rig) + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_rig Random number generator test rig + */ +template +void call_permutations(const T_functor& func, const T_rig& rig) { + call_permutations_helper + ::value - 1, + std::tuple_size::value - 1, + std::tuple_size::value - 1>{}(func, rig); +} + /* * Fill the vector-like variable params with values from the values argument. * @@ -83,339 +323,170 @@ void resize_if_vector(double& v, int N) { } /* - * check_dist_throws feeds the given generate_samples callable various + * check_dist_throws feeds rig.generate_samples various * combinations of valid and invalid parameters (as defined by good_p1_, bad_p1, * good_p2_, bad_p2_, good_p3_, and bad_p3_). For all combinations of valid * (good) parameters, generate_samples should throw no errors. For all * combinations with an invalid (bad) parameter, generate_samples should throw * domain_errors. * - * If good_p2_ or good_p3_ are empty, then it is assumed that those parameters + * If rig.good_p2_ or rig.good_p3_ are empty, then it is assumed that those parameters * are unused and will not be tested. * - * The generate_samples callable is also passed various other guaranteed invalid + * rig.generate_samples will also be passed various other guaranteed invalid * values like positive infinity, negative infinity, and NaNs (these should also * cause domain_errors). * - * The generated_samples callable is also tested to reject incompatibly sized - * input vectors. These should cause invalid_argument errors. - * - * The generate_samples callable must have the signature: - * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, T_rng) - * where T_param1, T_param2, and T_param3 can be any combination of: - * double, std::vector, VectorXd, or RowVectorXd - * and T_rng is a boost random number generator. p2 and p3 are required for the - * function signature, but they can be ignored if they are not needed. - * - * The output of generate_samples must be of - * stan::length(out) == stan::length(p1) if only p1 is used - * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used - * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used - * - * Each element of the output, out[i], should correspond to a random variate - * generated according to parameters p1[i], p2[i], and p3[i] + * rig.generated_samples is also tested to reject incompatibly sized input + * vectors. These should cause invalid_argument errors. * * @tparam T_param1 Type of first parameter * @tparam T_param2 Type of second parameter * @tparam T_param3 Type of third parameter - * @tparam T_generate_samples Type of generate_samples functor - * @param generate_samples generate_samples functor to test - * @param good_p1_ Valid values of first parameter - * @param bad_p1_ Invalid values of first parameter - * @param good_p2_ Valid values of second parameter - * (leave empty if p2 is unused) - * @param bad_p2_ Invalid values of second parameter - * @param good_p3_ Valid values of third parameter - * (leave empty if p3 is unused) - * @param bad_p3_ Invalid values of third parameter + * @tparam T_rig Type of test rig for random number generator + * @param T_rig Test rig for random number generator */ -template -void check_dist_throws(T_generate_samples generate_samples, - const std::vector& good_p1_, - const std::vector& bad_p1_, - const std::vector& good_p2_, - const std::vector& bad_p2_, - const std::vector& good_p3_, - const std::vector& bad_p3_) { - boost::random::mt19937 rng; - - T_param1 p1; - T_param2 p2; - T_param2 p3; - - bool p2_is_used = good_p2_.size() > 0; - bool p3_is_used = good_p3_.size() > 0; - - resize_if_vector(p1, 5); // No-op if p1 is a scalar - resize_if_vector(p2, 5); // No-op if p2 is a scalar - resize_if_vector(p3, 5); // No-op if p3 is a scalar - - // Make copies of the input arguments so that we can randomly shuffle them - // in the tests - std::vector good_p1 = good_p1_; - std::vector bad_p1 = bad_p1_; - std::vector good_p2 = good_p2_; - std::vector bad_p2 = bad_p2_; - std::vector good_p3 = good_p3_; - std::vector bad_p3 = bad_p3_; - - // Try a few combinations of parameters that should work - for (int i = 0; i < 5; i++) { - std::random_shuffle(good_p1.begin(), good_p1.end()); - std::random_shuffle(good_p2.begin(), good_p2.end()); - std::random_shuffle(good_p3.begin(), good_p3.end()); - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_NO_THROW(generate_samples(p1, p2, p3, rng)); - } +struct check_dist_throws { + template + void operator()(const T_rig& rig) const { + boost::random::mt19937 rng; + + T_param1 p1; + T_param2 p2; + T_param3 p3; + + bool p2_is_used = rig.good_p2_.size() > 0; + bool p3_is_used = rig.good_p3_.size() > 0; + + resize_if_vector(p1, 5); // No-op if p1 is a scalar + resize_if_vector(p2, 5); // No-op if p2 is a scalar + resize_if_vector(p3, 5); // No-op if p3 is a scalar + + // Make copies of the input arguments so that we can randomly shuffle them + // in the tests + std::vector good_p1 = rig.good_p1_; + std::vector bad_p1 = rig.bad_p1_; + std::vector good_p2 = rig.good_p2_; + std::vector bad_p2 = rig.bad_p2_; + std::vector good_p3 = rig.good_p3_; + std::vector bad_p3 = rig.bad_p3_; + + // Try a few combinations of parameters that should work + for (int i = 0; i < 5; i++) { + std::random_shuffle(good_p1.begin(), good_p1.end()); + std::random_shuffle(good_p2.begin(), good_p2.end()); + std::random_shuffle(good_p3.begin(), good_p3.end()); + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_NO_THROW(rig.generate_samples(p1, p2, p3, rng)); + } - // Now try putting incompatible values in first parameter - for (auto bad_p1_value : bad_p1) { - assign_parameter_values(p1, { bad_p1_value }); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); - } + // Now try putting incompatible values in first parameter + for (auto bad_p1_value : bad_p1) { + assign_parameter_values(p1, { bad_p1_value }); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + } - // Now try putting incompatible values in second parameter - for (auto bad_p2_value : bad_p2) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_p2_value }); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); - } + // Now try putting incompatible values in second parameter + for (auto bad_p2_value : bad_p2) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, { bad_p2_value }); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + } - // Now try putting incompatible values in third parameter - for (auto bad_p3_value : bad_p3) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, { bad_p3_value }); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); - } + // Now try putting incompatible values in third parameter + for (auto bad_p3_value : bad_p3) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, { bad_p3_value }); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + } + + // Make sure sampler throws errors with these values + std::vector bad_values = { stan::math::positive_infinity(), + stan::math::negative_infinity(), + stan::math::not_a_number() }; - // Make sure sampler throws errors with these values - std::vector bad_values = { stan::math::positive_infinity(), - stan::math::negative_infinity(), - stan::math::not_a_number() }; + for (auto bad_value : bad_values) { + assign_parameter_values(p1, { bad_value }); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + + if (p2_is_used) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, { bad_value }); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + } + + if (p3_is_used) { + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, { bad_value }); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); + } + } - for (auto bad_value : bad_values) { - assign_parameter_values(p1, { bad_value }); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + // Check to make sure _rng rejects vector arguments of different lengths for + // all parameter pairs. - if (p2_is_used) { + // If p1 is a scalar or the only vector, this test is skipped + resize_if_vector(p1, 3); // No-op if p1 is a scalar + resize_if_vector(p2, 4); // No-op if p2 is a scalar + resize_if_vector(p3, 4); // No-op if p3 is a scalar + if (stan::length(p1) != 1 && + ((p2_is_used && stan::length(p2) != 1) || + (p3_is_used && stan::length(p3) != 1))) { assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_value }); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); } - if (p3_is_used) { + // If p2 is a scalar or the only vector, this test is skipped + resize_if_vector(p1, 4); // No-op if p1 is a scalar + resize_if_vector(p2, 3); // No-op if p2 is a scalar + resize_if_vector(p3, 4); // No-op if p3 is a scalar + if (p2_is_used && stan::length(p2) != 1 && + (stan::length(p1) != 1 || (p3_is_used && stan::length(p3) != 1))) { + assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, { bad_value }); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::domain_error); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); } - } - - // Check to make sure _rng rejects vector arguments of different lengths for - // all parameter pairs. - - // If p1 is a scalar or the only vector, this test is skipped - resize_if_vector(p1, 3); // No-op if p1 is a scalar - resize_if_vector(p2, 4); // No-op if p2 is a scalar - resize_if_vector(p3, 4); // No-op if p3 is a scalar - if (stan::length(p1) != 1 && - ((p2_is_used && stan::length(p2) != 1) || - (p3_is_used && stan::length(p3) != 1))) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); - } - - // If p2 is a scalar or the only vector, this test is skipped - resize_if_vector(p1, 4); // No-op if p1 is a scalar - resize_if_vector(p2, 3); // No-op if p2 is a scalar - resize_if_vector(p3, 4); // No-op if p3 is a scalar - if (p2_is_used && stan::length(p2) != 1 && - (stan::length(p1) != 1 || (p3_is_used && stan::length(p3) != 1))) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); - } - // If p3 is a scalar or the only vector, this test is skipped - resize_if_vector(p1, 4); // No-op if p1 is a scalar - resize_if_vector(p2, 4); // No-op if p2 is a scalar - resize_if_vector(p3, 3); // No-op if p3 is a scalar - if (p3_is_used && stan::length(p3) != 1 && - (stan::length(p1) != 1 || (p2_is_used && stan::length(p2)))) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(generate_samples(p1, p2, p3, rng), std::invalid_argument); + // If p3 is a scalar or the only vector, this test is skipped + resize_if_vector(p1, 4); // No-op if p1 is a scalar + resize_if_vector(p2, 4); // No-op if p2 is a scalar + resize_if_vector(p3, 3); // No-op if p3 is a scalar + if (p3_is_used && stan::length(p3) != 1 && + (stan::length(p1) != 1 || (p2_is_used && stan::length(p2) != 1))) { + assign_parameter_values(p1, good_p1); + assign_parameter_values(p2, good_p2); + assign_parameter_values(p3, good_p3); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); + } } -} +}; /* - * This function calls check_dist_throws with the given generate_samples - * callable and all combinations of double, std::vector, - * Eigen::VectorXd, and Eigen::RowVectorXd as template arguments - * - * @tparam T_generate_samples Type of generate_samples functor - * @param generate_samples generate_samples functor to test - * @param good_p1 Valid values of first parameter - * @param bad_p1 Invalid values of first parameter - * @param good_p2 Valid values of second parameter - * @param bad_p2 Invalid values of second parameter - * @param good_p3 Valid values of third parameter - * @param bad_p3 Invalid values of third parameter + * This function calls check_dist_throws with the given test_rig + * and all combinations of double, std::vector, Eigen::VectorXd, + * and Eigen::RowVectorXd as template arguments + * + * @tparam T_rig Test rig type for random number generator + * @param T_rig Test rig for random number generator */ -template -void check_dist_throws_all_types(T_generate_samples generate_samples, - const std::vector& good_p1, - const std::vector& bad_p1, - const std::vector& good_p2, - const std::vector& bad_p2, - const std::vector& good_p3, - const std::vector& bad_p3) { +template +void check_dist_throws_all_types(const T_rig& rig) { using Eigen::VectorXd; using Eigen::RowVectorXd; - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double, double> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector, double > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd, double> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd, double> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws > - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, double, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, std::vector, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, VectorXd, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws, RowVectorXd> - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); - check_dist_throws - (generate_samples, good_p1, bad_p1, good_p2, bad_p2, good_p3, bad_p3); + call_permutations, VectorXd, + RowVectorXd> >(check_dist_throws{}, rig); } /* @@ -438,253 +509,80 @@ std::vector promote_to_vector(std::vector v) { /* * check_quantiles uses assert_matches_quantiles to check random numbers - * generated by the callable generate_samples against the reference - * quantiles computed with the given callable generate_quantiles. - * - * If good_p2 or good_p3 are empty, then it is assumed that those parameters - * are unused and will not be tested. - * - * The generate_samples callable must have the signature: - * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, p3, T_rng) - * where T_param1, T_param2, and T_param3 can be any combination of: - * double, std::vector, VectorXd, or RowVectorXd - * and T_rng is a boost random number generator. p2 and p3 are required for the - * function signature, but they can be ignored if they are not needed. - * - * The output of generate_samples must be of - * stan::length(out) == stan::length(p1) if only p1 is used - * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used - * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used - * - * Each element of the output, out[i], should correspond to a random variate - * generated according to parameters p1[i], p2[i], and p3[i] - * - * The generate_quantiles callable must have the signature: - * std::vector generate_samples(int N, double p1, double p2, double p3) - * (similarly to above, ignore parameters p2 and p3 if they aren't necessary) - * - * N is the number of samples that the quantiles will be compared - * against and p1, p2, and p3 are the values of the parameters that these - * quantiles correspond to. + * generated by rig.generate_samples against the reference quantiles computed + * with rig.generate_quantiles. * - * The output of generate_quantiles should be the quantiles used for comparing - * empirical distributions used by assert_matches_quantiles - * - * Parameters for each of these callables are copied from the vectors good_p1, - * good_p2, and good_p3. + * If rig.good_p2_ or rig.good_p3_ are empty, then it is assumed that those + * parameters are unused and will not be tested. * * @tparam T_param1 Type of first parameter * @tparam T_param2 Type of second parameter * @tparam T_param3 Type of third parameter - * @tparam T_generate_samples Type of callable for generating samples - * @tparam T_generate_quantiles Type of callable for generating quantiles - * @param sampler sampler functor to test - * @param good_p1 Valid values of first parameter - * @param good_p2 Valid values of second parameter - * @param good_p3 Valid values of third parameter + * @tparam T_rig Type of test rig for random number generator + * @param T_rig Test rig for random number generator */ -template -void check_quantiles(int N, int M_vec, - T_generate_samples generate_samples, - T_generate_quantiles generate_quantiles, - const std::vector& good_p1, - const std::vector& good_p2, - const std::vector& good_p3) { - boost::random::mt19937 rng; - T_param1 p1; - T_param2 p2; - T_param3 p3; - resize_if_vector(p1, M_vec); // No-op if p1 is scalar - resize_if_vector(p2, M_vec); // No-op if p2 is scalar - resize_if_vector(p3, M_vec); // No-op if p3 is scalar - - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - - bool p2_is_used = good_p2.size() > 0; - bool p3_is_used = good_p3.size() > 0; - - int M = std::max({ stan::length(p1), - (p2_is_used) ? stan::length(p2) : 1, - (p3_is_used) ? stan::length(p3) : 1 }); - - stan::scalar_seq_view p1_vec(p1); - stan::scalar_seq_view p2_vec(p2); - stan::scalar_seq_view p3_vec(p3); - - std::vector > samples_to_test_transpose; - for (int n = 0; n < N; ++n) { - // If p1, p2, and p3 are scalars, the output is a scalar. Need to promote it - // to a std::vector - samples_to_test_transpose. - push_back(promote_to_vector(generate_samples(p1, p2, p3, rng))); - } - - for (int m = 0; m < M; ++m) { - std::vector samples_to_test; - for (int n = 0; n < N; ++n) { - samples_to_test.push_back(samples_to_test_transpose[n][m]); +struct check_quantiles { + template + void operator()(const T_rig& rig) const { + boost::random::mt19937 rng; + T_param1 p1; + T_param2 p2; + T_param3 p3; + resize_if_vector(p1, rig.M_); // No-op if p1 is scalar + resize_if_vector(p2, rig.M_); // No-op if p2 is scalar + resize_if_vector(p3, rig.M_); // No-op if p3 is scalar + + assign_parameter_values(p1, rig.good_p1_); + assign_parameter_values(p2, rig.good_p2_); + assign_parameter_values(p3, rig.good_p3_); + + bool p2_is_used = rig.good_p2_.size() > 0; + bool p3_is_used = rig.good_p3_.size() > 0; + + int M = std::max({ stan::length(p1), + (p2_is_used) ? stan::length(p2) : 1, + (p3_is_used) ? stan::length(p3) : 1 }); + + stan::scalar_seq_view p1_vec(p1); + stan::scalar_seq_view p2_vec(p2); + stan::scalar_seq_view p3_vec(p3); + + std::vector > samples_to_test_transpose; + for (int n = 0; n < rig.N_; ++n) { + // If p1, p2, and p3 are scalars, the output is a scalar. Need to promote it + // to a std::vector + samples_to_test_transpose. + push_back(promote_to_vector(rig.generate_samples(p1, p2, p3, rng))); } - std::vector quantiles = - generate_quantiles(N, p1_vec[m], p2_vec[m], p3_vec[m]); - assert_matches_quantiles(samples_to_test, quantiles, 1e-6); + for (int m = 0; m < M; ++m) { + std::vector samples_to_test; + for (int n = 0; n < rig.N_; ++n) { + samples_to_test.push_back(samples_to_test_transpose[n][m]); + } + std::vector quantiles = + rig.generate_quantiles(p1_vec[m], p2_vec[m], p3_vec[m]); + + assert_matches_quantiles(samples_to_test, quantiles, 1e-6); + } } -} +}; /* - * Call check_quantiles with the given arguments and all combinations of + * Call check_quantiles on the given test rig with all combinations of * double, std::vector, Eigen::VectorXd, and Eigen::RowVectorXd * - * @tparam T_generate_samples Type of callable for generating samples - * @tparam T_generate_quantiles Type of callable for generating quantiles - * @param sampler sampler functor to test - * @param good_p1 Valid values of first parameter - * @param good_p2 Valid values of second parameter + * @tparam T_rig Type of test rig for random number generator + * @param T_rig Test rig for random number generator */ -template -void check_quantiles_all_types(int N, int M, - T_generate_samples generate_samples, - T_generate_quantiles generate_quantiles, - const std::vector& good_p1, - const std::vector& good_p2, - const std::vector& good_p3) { +template +void check_quantiles_all_types(const T_rig& rig) { using Eigen::VectorXd; using Eigen::RowVectorXd; - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles > - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, double, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, std::vector, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, VectorXd, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles, RowVectorXd> - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); - check_quantiles - (N, M, generate_samples, generate_quantiles, good_p1, good_p2, good_p3); + call_permutations, VectorXd, + RowVectorXd> >(check_quantiles{}, rig); } #endif From d2281d17094e617148d763f8369354ec1627a852 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Wed, 4 Oct 2017 17:34:57 -0700 Subject: [PATCH 05/12] Added int and std::vector tests (Issue 45) --- .../math/prim/mat/prob/VectorRNGTestRig.hpp | 255 ++++++++++++ test/unit/math/prim/mat/prob/cauchy_test.cpp | 10 +- .../prim/mat/prob/double_exponential_test.cpp | 10 +- .../prim/mat/prob/exp_mod_normal_test.cpp | 15 +- test/unit/math/prim/mat/prob/gumbel_test.cpp | 10 +- .../unit/math/prim/mat/prob/logistic_test.cpp | 10 +- test/unit/math/prim/mat/prob/normal_test.cpp | 14 +- .../math/prim/mat/prob/skew_normal_test.cpp | 15 +- .../math/prim/mat/prob/student_t_test.cpp | 15 +- .../prim/mat/prob/vector_rng_test_helper.hpp | 383 +++++------------- .../scal/meta/apply_template_permutations.hpp | 124 ++++++ 11 files changed, 547 insertions(+), 314 deletions(-) create mode 100644 test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp create mode 100644 test/unit/math/prim/scal/meta/apply_template_permutations.hpp diff --git a/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp new file mode 100644 index 00000000000..4d8e0317b20 --- /dev/null +++ b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp @@ -0,0 +1,255 @@ +#ifndef TEST_UNIT_MATH_PRIM_MAT_PROB_VECTOR_RNG_TEST_RIG_HPP +#define TEST_UNIT_MATH_PRIM_MAT_PROB_VECTOR_RNG_TEST_RIG_HPP + +#include +#include +#include + +/* + * VectorRNGTestRig is an Abstract class for wrapping up random number + * generators in a way so that they can be tested to + * 1) provide interfaces to all the necessary types + * 2) generate the correct distributions + * + * A test rig that inherits from VectorRNGTestRig must implement + * 1) generate_samples + * 2) generate_quantiles + * and must initialize good_p1_, good_p1_int_, bad_p1_, bad_p1_int_, good_p2_, + * good_p2_int_, bad_p2_, bad_p2_int_, good_p3_, good_p3_int_, and bad_p3_ + * bad_p3_int_ (depending on how many parameters the wrapped random number + * generator uses) + * + * generate_samples is a wrapper around the actual random number generator. It + * is expected to have the signature: + * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, T_rng) + * where T_param1, T_param2, and T_param3 can be any combination of: + * int, double, std::vector, std::vector, VectorXd, & RowVectorXd + * and T_rng is a boost random number generator. p2 and p3 are required for the + * function signature, but they can be ignored if they are not needed. + * + * Each element of the output, out[i], should correspond to a random variate + * generated according to parameters p1[i], p2[i], and p3[i] + * + * generate_samples is expected to throw domain_errors when passed invalid + * values like positive infinity, negative infinity, and NaNs. + * + * generated_samples is also expected to reject incompatibly sized + * input vectors. These should cause invalid_argument errors. + * + * The output of generate_samples must be of + * stan::length(out) == stan::length(p1) if only p1 is used + * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used + * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used + * It must be defined to take three arguments, but not all must be used. + * + * The generate_quantiles callable must have the signature: + * std::vector generate_samples(double p1, double p2, double p3) + * (similarly to above, ignore parameters p2 and p3 if they aren't necessary) + * + * generate_quantiles should compute the quantiles that will be compared against + * (using assert_matches_quantiles) empirical distributions generated with + * generate_samples (which will be of length N_ for each parameter combination). + * + * p1, p2, and p3 are the values of the parameters that these quantiles + * correspond to. + * + * good_p1_ and bad_p1_ should be initialized to lists of valid and invalid + * floating point parameters for the first parameter of the tested RNG + * + * good_p1_int_ and bad_p1_int_ should be initialized to lists of valid and + * invalid integer parameters. + * + * good_p2_, good_p2_int_, bad_p2_, bad_p2_int_, good_p3_, good_p3_int_, and + * bad_p3_int_ should be initialized similarly. + * + * If testing a distribution that uses only two parameters, leave good_p3_ and + * good_p3_int as length zero vectors and it will not be tested. Similarly, if + * good_p2_ and good_p2_int_ have length zero, the second parameter will not be + * tested either. + */ +class VectorRNGTestRig { +public: + int N_; // Number of samples used in the quantiles tests + int M_; // Length of vectors for the vectorization tests + + std::vector good_p1_; + std::vector good_p1_int_; + std::vector bad_p1_; + std::vector bad_p1_int_; + std::vector good_p2_; + std::vector good_p2_int_; + std::vector bad_p2_; + std::vector bad_p2_int_; + std::vector good_p3_; + std::vector good_p3_int_; + std::vector bad_p3_; + std::vector bad_p3_int_; + + std::vector always_bad_values_ = { stan::math::positive_infinity(), + stan::math::negative_infinity(), + stan::math::not_a_number() }; + /* + * This function wraps up the random number generator for testing. + * + * The tested rng can have up to three parameters. Any unused parameters can + * be ignored. + * + * It *must* be implemented in the child class (it isn't virtual here because + * C++ doesn't allow templated virtual member functions) + */ + template + auto generate_samples(const T1& p1, const T2& p2, const T3& p3, + T_rng& rng) const; + + /* + * This function builds the quantiles that we will supply to + * assert_matches_quantiles to test the random number generator implemented in + * generate_samples + */ + virtual std::vector generate_quantiles(double p1, double p2, + double p3) const = 0; + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector good_p1_int, + std::vector bad_p1, + std::vector bad_p1_int, + std::vector good_p2, + std::vector good_p2_int, + std::vector bad_p2, + std::vector bad_p2_int, + std::vector good_p3, + std::vector good_p3_int, + std::vector bad_p3, + std::vector bad_p3_int) : N_(N), M_(M), + good_p1_(good_p1), + good_p1_int_(good_p1_int), + bad_p1_(bad_p1), + bad_p1_int_(bad_p1_int), + good_p2_(good_p2), + good_p2_int_(good_p2_int), + bad_p2_(bad_p2), + bad_p2_int_(bad_p2_int), + good_p3_(good_p3), + good_p3_int_(good_p3_int), + bad_p3_(bad_p3), + bad_p3_int_(bad_p3_int) { + } + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector good_p1_int, + std::vector bad_p1, + std::vector bad_p1_int, + std::vector good_p2, + std::vector good_p2_int, + std::vector bad_p2, + std::vector bad_p2_int) : N_(N), M_(M), + good_p1_(good_p1), + good_p1_int_(good_p1_int), + bad_p1_(bad_p1), + bad_p1_int_(bad_p1_int), + good_p2_(good_p2), + good_p2_int_(good_p2_int), + bad_p2_(bad_p2), + bad_p2_int_(bad_p2_int) { + } + + VectorRNGTestRig(int N, int M, + std::vector good_p1, + std::vector good_p1_int, + std::vector bad_p1, + std::vector bad_p1_int) : N_(N), M_(M), + good_p1_(good_p1), + good_p1_int_(good_p1_int), + bad_p1_(bad_p1), + bad_p1_int_(bad_p1_int) { + } + + /* + * If no good values for p2 or p3 are provided, it is assumed that those + * parameters are unused + */ + bool p2_is_used() const { + return good_p2_.size() > 0 || good_p2_int_.size() > 0; + } + + bool p3_is_used() const { + return good_p3_.size() > 0 || good_p3_int_.size() > 0; + } + + /* + * All the accessors work the same. Based on the input template type, they + * return the list of valid/invalid values for whichever parameter + * is accessed. + * + * The bad value accessors always append on the list of always_bad_values + * to whatever was already specified (always_bad_values contains infs and + * nans and stuff that should always cause errors). + * + * @tparam T Template parameter determining which list of parameters to return + * (can be double or int) + * @return List of parameter values + */ + template + std::vector get_good_p1() const { + return good_p1_; + } + + template + std::vector get_bad_p1() const { + return stan::math::append_array(bad_p1_, always_bad_values_); + } + + template + std::vector get_good_p2() const { + return good_p2_; + } + + template + std::vector get_bad_p2() const { + return stan::math::append_array(bad_p2_, always_bad_values_); + } + + template + std::vector get_good_p3() const { + return good_p3_; + } + + template + std::vector get_bad_p3() const { + return stan::math::append_array(bad_p3_, always_bad_values_); + } +}; + +template<> +std::vector VectorRNGTestRig::get_good_p1() const { + return good_p1_int_; +} + +template<> +std::vector VectorRNGTestRig::get_bad_p1() const { + return bad_p1_int_; +} + +template<> +std::vector VectorRNGTestRig::get_good_p2() const { + return good_p2_int_; +} + +template<> +std::vector VectorRNGTestRig::get_bad_p2() const { + return bad_p2_int_; +} + +template<> +std::vector VectorRNGTestRig::get_good_p3() const { + return good_p3_int_; +} + +template<> +std::vector VectorRNGTestRig::get_bad_p3() const { + return bad_p3_int_; +} + +#endif diff --git a/test/unit/math/prim/mat/prob/cauchy_test.cpp b/test/unit/math/prim/mat/prob/cauchy_test.cpp index 35f5dabc7d8..f239e24dd2a 100644 --- a/test/unit/math/prim/mat/prob/cauchy_test.cpp +++ b/test/unit/math/prim/mat/prob/cauchy_test.cpp @@ -10,8 +10,14 @@ class CauchyTestRig : public VectorRNGTestRig { public: CauchyTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, diff --git a/test/unit/math/prim/mat/prob/double_exponential_test.cpp b/test/unit/math/prim/mat/prob/double_exponential_test.cpp index ba1f76ca576..f7f9ecf8514 100644 --- a/test/unit/math/prim/mat/prob/double_exponential_test.cpp +++ b/test/unit/math/prim/mat/prob/double_exponential_test.cpp @@ -10,8 +10,14 @@ class DoubleExponentialTestRig : public VectorRNGTestRig { public: DoubleExponentialTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, diff --git a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp index 99752ae0e22..d6949481690 100644 --- a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp @@ -6,9 +6,18 @@ class ExpModNormalTestRig : public VectorRNGTestRig { public: ExpModNormalTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {0.2, 1.1, 2.5, 2.0}, {-3.7, -2.5, -1.0, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}, + {0.2, 1.1, 2.5, 2.0}, + {1, 2, 3, 4}, + {-3.7, -2.5, -1.0, 0.0}, + {-4, -3, -2, 0}) {} template auto generate_samples(const T1& loc, const T2& scale, const T3& inv_scale, diff --git a/test/unit/math/prim/mat/prob/gumbel_test.cpp b/test/unit/math/prim/mat/prob/gumbel_test.cpp index 223466de8b4..d24ea894fb5 100644 --- a/test/unit/math/prim/mat/prob/gumbel_test.cpp +++ b/test/unit/math/prim/mat/prob/gumbel_test.cpp @@ -10,8 +10,14 @@ class GumbelTestRig : public VectorRNGTestRig { public: GumbelTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-1.0, -1.5, -2.5, -0.7, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-1.0, -1.5, -2.5, -0.7, 0.0}, + {-1, -2, -3, -4, 0}) {} template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, diff --git a/test/unit/math/prim/mat/prob/logistic_test.cpp b/test/unit/math/prim/mat/prob/logistic_test.cpp index fb211e11183..7c97100e7a1 100644 --- a/test/unit/math/prim/mat/prob/logistic_test.cpp +++ b/test/unit/math/prim/mat/prob/logistic_test.cpp @@ -10,8 +10,14 @@ class LogisticTestRig : public VectorRNGTestRig { public: LogisticTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, diff --git a/test/unit/math/prim/mat/prob/normal_test.cpp b/test/unit/math/prim/mat/prob/normal_test.cpp index cc7c93191ac..4a5213946ae 100644 --- a/test/unit/math/prim/mat/prob/normal_test.cpp +++ b/test/unit/math/prim/mat/prob/normal_test.cpp @@ -3,12 +3,10 @@ #include #include #include +#include #include #include -using Eigen::Dynamic; -using Eigen::Matrix; - class NormalTestRig : public VectorRNGTestRig { public: /* @@ -19,8 +17,14 @@ class NormalTestRig : public VectorRNGTestRig { NormalTestRig() : VectorRNGTestRig(10000, // Number of samples used for quantiles tests 10, // Length of vectors for vectorization tests - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, // Valid values for p1 + {-3, -2, -1, 0, 2, 6}, // Valid integer values for p1 + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} /* * This function wraps up the random number generator for testing. diff --git a/test/unit/math/prim/mat/prob/skew_normal_test.cpp b/test/unit/math/prim/mat/prob/skew_normal_test.cpp index dd2c448f622..c3615c15017 100644 --- a/test/unit/math/prim/mat/prob/skew_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/skew_normal_test.cpp @@ -11,9 +11,18 @@ class SkewNormalTestRig : public VectorRNGTestRig { public: SkewNormalTestRig() : VectorRNGTestRig(10000, 10, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}, - {-2.5, -1.7, -1.3, 0.0, 1.0, 4.8}, {}) {} + {-2.5, -1.7, -0.1, 0.0, 2.0}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}, + {-2.0, -1.0, -0.5, 0.0, 0.7, 0.5}, + {-2, -1, 0, 1, 2}, + {}, + {}) {} template auto generate_samples(const T1& mu, const T2& sigma, const T3& alpha, diff --git a/test/unit/math/prim/mat/prob/student_t_test.cpp b/test/unit/math/prim/mat/prob/student_t_test.cpp index 7e10a7a9e17..ba0ee953fd3 100644 --- a/test/unit/math/prim/mat/prob/student_t_test.cpp +++ b/test/unit/math/prim/mat/prob/student_t_test.cpp @@ -10,9 +10,18 @@ class StudentTTestRig : public VectorRNGTestRig { public: StudentTTestRig() : VectorRNGTestRig(10000, 10, - {1.1, 2.0, 2.5, 2.0}, {-1.7, -0.5, -2.5, 0.0}, - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, {}, - {0.1, 1.0, 2.5, 4.0}, {-2.7, -1.5, -0.5, 0.0}) {} + {1.1, 2.0, 2.5, 3.1}, + {1, 2, 3, 4}, + {-1.7, -0.5, -2.5, 0.0}, + {-2, -1, -3, 0}, + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, + {-3, -2, -1, 0, 2, 6}, + {}, + {}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} template auto generate_samples(const T1& nu, const T2& mu, const T3& sigma, diff --git a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp index faea885789b..df01cfbc124 100644 --- a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp +++ b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp @@ -6,271 +6,52 @@ #include #include #include +#include +#include #include /* - * VectorRNGTestRig is an Abstract class for wrapping up random number - * generators in a way so that they can be tested to - * 1) provide interfaces to all the necessary types - * 2) generate the correct distributions - * - * A test rig that inherits from VectorRNGTestRig must implement - * 1) generate_samples - * 2) generate_quantiles - * and must initialize good_p1_, bad_p1_, good_p2_, bad_p2, good_p3_, bad_p3_ - * (depending on how many parameters the wrapped random number generator uses) - * - * generate_samples is a wrapper around the actual random number generator. It - * is expected to have the signature: - * T_out generate_samples(T_param1 p1, T_param2 p2, T_param3, T_rng) - * where T_param1, T_param2, and T_param3 can be any combination of: - * double, std::vector, VectorXd, or RowVectorXd - * and T_rng is a boost random number generator. p2 and p3 are required for the - * function signature, but they can be ignored if they are not needed. - * - * Each element of the output, out[i], should correspond to a random variate - * generated according to parameters p1[i], p2[i], and p3[i] - * - * generate_samples is expected to throw domain_errors when passed invalid - * values like positive infinity, negative infinity, and NaNs. - * - * generated_samples is also expected to reject incompatibly sized - * input vectors. These should cause invalid_argument errors. - * - * The output of generate_samples must be of - * stan::length(out) == stan::length(p1) if only p1 is used - * stan::length(out) == stan::max_size(p1, p2) if p1 and p2 are used - * stan::length(out) == stan::max_size(p1, p2, p3) if all parameters are used - * It must be defined to take three arguments, but not all must be used. - * - * The generate_quantiles callable must have the signature: - * std::vector generate_samples(double p1, double p2, double p3) - * (similarly to above, ignore parameters p2 and p3 if they aren't necessary) - * - * generate_quantiles should compute the quantiles that will be compared against - * (using assert_matches_quantiles) empirical distributions generated with - * generate_samples (which will be of length N_ for each parameter combination). - * - * p1, p2, and p3 are the values of the parameters that these quantiles - * correspond to. - * - * good_p1_ and bad_p1_ should be initialized to lists of valid and invalid - * parameters for the first parameter of the tested RNG - * - * Likewise, good_p2_, bad_p2_, and good_p3_, bad_p3_ should be initialized. - * - * If testing a distribution that uses only two parameters, leave good_p3_ as - * a length zero vector and it will not be tested. Similarly, if good_p2_ is - * length zero, the second parameter will not be tested either - */ -class VectorRNGTestRig { -public: - int N_; // Number of samples used in the quantiles tests - int M_; // Length of vectors for the vectorization tests - - std::vector good_p1_; - std::vector bad_p1_; - std::vector good_p2_; - std::vector bad_p2_; - std::vector good_p3_; - std::vector bad_p3_; - - /* - * This function wraps up the random number generator for testing. - * - * The tested rng can have up to three parameters. Any unused parameters can - * be ignored. - * - * It *must* be implemented in the child class (it isn't virtual here because - * C++ doesn't allow templated virtual member functions) - */ - template - auto generate_samples(const T1& p1, const T2& p2, const T3& p3, - T_rng& rng) const; - - /* - * This function builds the quantiles that we will supply to - * assert_matches_quantiles to test the random number generator implemented in - * generate_samples - */ - virtual std::vector generate_quantiles(double p1, double p2, - double p3) const = 0; - - VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector bad_p1, - std::vector good_p2, - std::vector bad_p2, - std::vector good_p3, - std::vector bad_p3) : N_(N), M_(M), - good_p1_(good_p1), - bad_p1_(bad_p1), - good_p2_(good_p2), - bad_p2_(bad_p2), - good_p3_(good_p3), - bad_p3_(bad_p3) { - } - - VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector bad_p1, - std::vector good_p2, - std::vector bad_p2) : N_(N), M_(M), - good_p1_(good_p1), - bad_p1_(bad_p1), - good_p2_(good_p2), - bad_p2_(bad_p2){ - } - - VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector bad_p1) : N_(N), M_(M), - good_p1_(good_p1), - bad_p1_(bad_p1) { - } -}; - -/* - * call_permutations_helper is the primary function for the recursive template - * function call_permutations. It calls func(rig), - * and continues the recursion by decrementing K. - * - * tparam T_typelist Tuple templated with list of types to iterate through - * tparam T_functor Templated functor - * tparam T_rig Random number generator test rig - * tparam I Loop index for first type - * tparam J Loop index for second type - * tparam K Loop index for third type - */ -template -struct call_permutations_helper { - void operator()(const T_functor& func, const T_rig& rig) const { - func.template operator()::type, - typename std::tuple_element::type, - typename std::tuple_element::type> - (rig); - - call_permutations_helper{}(func, - rig); - } -}; - -/* - * This edge case catches when the right-most argument has completed one - * iteration through all types. It computes func, - * carries from the second argument, and continues the recursion. + * Fill the vector-like variable params with values from the values argument. * - * tparam T_typelist Tuple templated with list of types to iterate through - * tparam T_functor Templated functor - * tparam T_rig Random number generator test rig - * tparam I Loop index for first type - * tparam J Loop index for second type - * tparam K Loop index for third type - */ -template -struct call_permutations_helper { - void operator()(const T_functor& func, const T_rig& rig) const { - func.template operator()::type, - typename std::tuple_element::type, - typename std::tuple_element<0, T_typelist>::type> - (rig); - - call_permutations_helper::value - 1>{}(func, rig); - } -}; - -/* - * This edge case catches when the second argument has completed one iteration - * through all types. It computes func, carries from - * the left-most argument and continues the recursion. + * values can be shorter than params (in which cause multiple copies of values + * are tiled over the params vector) * - * tparam T_typelist Tuple templated with list of types to iterate through - * tparam T_functor Templated functor - * tparam T_rig Random number generator test rig - * tparam I Loop index for first type + * @tparam T_param Type of params vector + * @param params Parameter vector to write values to + * @param params Values to copy into params */ -template -struct call_permutations_helper { - void operator()(const T_functor& func, const T_rig& rig) const { - func.template operator()::type, - typename std::tuple_element<0, T_typelist>::type, - typename std::tuple_element<0, T_typelist>::type> - (rig); - - call_permutations_helper::value - 1, - std::tuple_size::value - 1>{}(func, - rig); - } -}; +template +void assign_parameter_values(T_param& params, + const std::vector& values) { + if (values.size() == 0) + return; -/* - * This edge case catches when all types have been iterated through for each - * template argument. It computes func and ends the - * recursion. - * - * tparam T_typelist Tuple templated with list of types to iterate through - * tparam T_functor Templated functor - * tparam T_rig Random number generator test rig - */ -template -struct call_permutations_helper { - void operator()(const T_functor& func, const T_rig& rig) const { - func.template operator()::type, - typename std::tuple_element<0, T_typelist>::type, - typename std::tuple_element<0, T_typelist>::type> - (rig); + for (int i = 0; i < params.size(); i++) { + params[i] = values[i % values.size()]; } -}; - - -/* - * call_permutations uses template recursion to call the functor func - * with all 3-element permutations of the types in T_typelist - * - * Roughly, this corresponds to: - * for (i in (I - 1):0) // indexes go backwards, like: I - 1, I - 2, ... 0 - * for (j in (J - 1):0) - * for (k in (K - 1):0) - * func(rig) - * - * tparam T_typelist Tuple templated with list of types to iterate through - * tparam T_functor Templated functor - * tparam T_rig Random number generator test rig - */ -template -void call_permutations(const T_functor& func, const T_rig& rig) { - call_permutations_helper - ::value - 1, - std::tuple_size::value - 1, - std::tuple_size::value - 1>{}(func, rig); } /* - * Fill the vector-like variable params with values from the values argument. + * Fill the vector params with values from the values argument. * * values can be shorter than params (in which cause multiple copies of values * are tiled over the params vector) * - * @tparam T_param Type of params vector * @param params Parameter vector to write values to * @param params Values to copy into params */ -template -void assign_parameter_values(T_param& params, +void assign_parameter_values(std::vector& params, const std::vector& values) { if (values.size() == 0) return; - for (int i = 0; i < params.size(); i++) { + for (size_t i = 0; i < params.size(); i++) { params[i] = values[i % values.size()]; } } /* - * Fill the vector-like variable params with values from the values argument. + * Fill the vector params with values from the values argument. * * values can be shorter than params (in which cause multiple copies of values * are tiled over the params vector) @@ -278,9 +59,8 @@ void assign_parameter_values(T_param& params, * @param params Parameter vector to write values to * @param params Values to copy into params */ -template <> -void assign_parameter_values(std::vector& params, - const std::vector& values) { +void assign_parameter_values(std::vector& params, + const std::vector& values) { if (values.size() == 0) return; @@ -295,7 +75,6 @@ void assign_parameter_values(std::vector& params, * @param param Output parameter to write value to * @param params Vector with value to copy into param */ -template <> void assign_parameter_values(double& param, const std::vector& values) { if (values.size() == 0) return; @@ -303,6 +82,19 @@ void assign_parameter_values(double& param, const std::vector& values) { param = values[0]; } +/* + * Assign param the first value of values + * + * @param param Output parameter to write value to + * @param params Vector with value to copy into param + */ +void assign_parameter_values(int& param, const std::vector& values) { + if (values.size() == 0) + return; + + param = values[0]; +} + /* * Resize v to be length N * @@ -322,6 +114,13 @@ template<> void resize_if_vector(double& v, int N) { } +/* + * For ints, resize_if_vector does nothing + */ +template<> +void resize_if_vector(int& v, int N) { +} + /* * check_dist_throws feeds rig.generate_samples various * combinations of valid and invalid parameters (as defined by good_p1_, bad_p1, @@ -355,8 +154,12 @@ struct check_dist_throws { T_param2 p2; T_param3 p3; - bool p2_is_used = rig.good_p2_.size() > 0; - bool p3_is_used = rig.good_p3_.size() > 0; + using T_scalar_param1 = typename stan::scalar_type::type; + using T_scalar_param2 = typename stan::scalar_type::type; + using T_scalar_param3 = typename stan::scalar_type::type; + + bool p2_is_used = rig.p2_is_used(); + bool p3_is_used = rig.p3_is_used(); resize_if_vector(p1, 5); // No-op if p1 is a scalar resize_if_vector(p2, 5); // No-op if p2 is a scalar @@ -364,12 +167,18 @@ struct check_dist_throws { // Make copies of the input arguments so that we can randomly shuffle them // in the tests - std::vector good_p1 = rig.good_p1_; - std::vector bad_p1 = rig.bad_p1_; - std::vector good_p2 = rig.good_p2_; - std::vector bad_p2 = rig.bad_p2_; - std::vector good_p3 = rig.good_p3_; - std::vector bad_p3 = rig.bad_p3_; + std::vector good_p1 = + rig.template get_good_p1(); + std::vector bad_p1 = + rig.template get_bad_p1(); + std::vector good_p2 = + rig.template get_good_p2(); + std::vector bad_p2 = + rig.template get_bad_p2(); + std::vector good_p3 = + rig.template get_good_p3(); + std::vector bad_p3 = + rig.template get_bad_p3(); // Try a few combinations of parameters that should work for (int i = 0; i < 5; i++) { @@ -391,41 +200,21 @@ struct check_dist_throws { } // Now try putting incompatible values in second parameter - for (auto bad_p2_value : bad_p2) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_p2_value }); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); - } - - // Now try putting incompatible values in third parameter - for (auto bad_p3_value : bad_p3) { - assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, { bad_p3_value }); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); - } - - // Make sure sampler throws errors with these values - std::vector bad_values = { stan::math::positive_infinity(), - stan::math::negative_infinity(), - stan::math::not_a_number() }; - - for (auto bad_value : bad_values) { - assign_parameter_values(p1, { bad_value }); - assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, good_p3); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); - - if (p2_is_used) { + if(p2_is_used) { + for (auto bad_p2_value : bad_p2) { assign_parameter_values(p1, good_p1); - assign_parameter_values(p2, { bad_value }); + assign_parameter_values(p2, { bad_p2_value }); + assign_parameter_values(p3, good_p3); EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); } + } - if (p3_is_used) { + // Now try putting incompatible values in third parameter + if(p3_is_used) { + for (auto bad_p3_value : bad_p3) { + assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); - assign_parameter_values(p3, { bad_value }); + assign_parameter_values(p3, { bad_p3_value }); EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::domain_error); } } @@ -485,25 +274,30 @@ void check_dist_throws_all_types(const T_rig& rig) { using Eigen::VectorXd; using Eigen::RowVectorXd; - call_permutations, VectorXd, - RowVectorXd> >(check_dist_throws{}, rig); + apply_template_permutations, double, + std::vector, VectorXd, + RowVectorXd> > + (check_dist_throws{}, rig); } /* * Convert a scalar to a length one vector * + * @tparam T Scalar type * @param v Input scalar * @return vector of length 1 with value v */ -std::vector promote_to_vector(double v) { - return std::vector(1, v); +template +std::vector promote_to_vector(T v) { + return std::vector(1, v); } /* * For arguments that are already vectors, just copy. Probably more efficient * just using std::move but cpplint complained about use of unapproved Rvalues. */ -std::vector promote_to_vector(std::vector v) { +template +std::vector promote_to_vector(std::vector v) { return v; } @@ -533,12 +327,15 @@ struct check_quantiles { resize_if_vector(p2, rig.M_); // No-op if p2 is scalar resize_if_vector(p3, rig.M_); // No-op if p3 is scalar - assign_parameter_values(p1, rig.good_p1_); - assign_parameter_values(p2, rig.good_p2_); - assign_parameter_values(p3, rig.good_p3_); + assign_parameter_values(p1, rig.template get_good_p1 + ::type>()); + assign_parameter_values(p2, rig.template get_good_p2 + ::type>()); + assign_parameter_values(p3, rig.template get_good_p3 + ::type>()); - bool p2_is_used = rig.good_p2_.size() > 0; - bool p3_is_used = rig.good_p3_.size() > 0; + bool p2_is_used = rig.p2_is_used(); + bool p3_is_used = rig.p3_is_used(); int M = std::max({ stan::length(p1), (p2_is_used) ? stan::length(p2) : 1, @@ -581,8 +378,10 @@ void check_quantiles_all_types(const T_rig& rig) { using Eigen::VectorXd; using Eigen::RowVectorXd; - call_permutations, VectorXd, - RowVectorXd> >(check_quantiles{}, rig); + apply_template_permutations, + std::vector, VectorXd, + RowVectorXd > > + (check_quantiles{}, rig); } #endif diff --git a/test/unit/math/prim/scal/meta/apply_template_permutations.hpp b/test/unit/math/prim/scal/meta/apply_template_permutations.hpp new file mode 100644 index 00000000000..19e8d911d73 --- /dev/null +++ b/test/unit/math/prim/scal/meta/apply_template_permutations.hpp @@ -0,0 +1,124 @@ +#ifndef TEST_UNIT_MATH_PRIM_SCAL_META_APPLY_TEMPLATE_PERMUTATIONS_HPP +#define TEST_UNIT_MATH_PRIM_SCAL_META_APPLY_TEMPLATE_PERMUTATIONS_HPP + +#include + +/* + * apply_template_permutations_helper is the entry point for the template recursion + * needed for apply_template_permutations. It calls func(rig), + * and continues the recursion by decrementing K. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_param Parameter to pass to functor + * tparam I Loop index for first type + * tparam J Loop index for second type + * tparam K Loop index for third type + */ +template +struct apply_template_permutations_helper { + void operator()(const T_functor& func, const T_param& param) const { + func.template operator()::type, + typename std::tuple_element::type, + typename std::tuple_element::type> + (param); + + apply_template_permutations_helper{}(func, + param); + } +}; + +/* + * This edge case catches when the paramht-most argument has completed one + * iteration through all types. It computes func, + * carries from the second argument, and continues the recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_param Parameter to pass to functor + * tparam I Loop index for first type + * tparam J Loop index for second type + * tparam K Loop index for third type + */ +template +struct apply_template_permutations_helper { + void operator()(const T_functor& func, const T_param& param) const { + func.template operator()::type, + typename std::tuple_element::type, + typename std::tuple_element<0, T_typelist>::type> + (param); + + apply_template_permutations_helper::value - 1>{}(func, param); + } +}; + +/* + * This edge case catches when the second argument has completed one iteration + * through all types. It computes func, carries from + * the left-most argument and continues the recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_param Parameter to pass to functor + * tparam I Loop index for first type + */ +template +struct apply_template_permutations_helper { + void operator()(const T_functor& func, const T_param& param) const { + func.template operator()::type, + typename std::tuple_element<0, T_typelist>::type, + typename std::tuple_element<0, T_typelist>::type> + (param); + + apply_template_permutations_helper::value - 1, + std::tuple_size::value - 1>{}(func, + param); + } +}; + +/* + * This edge case catches when all types have been iterated through for each + * template argument. It computes func and ends the + * recursion. + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_param Parameter to pass to functor + */ +template +struct apply_template_permutations_helper { + void operator()(const T_functor& func, const T_param& param) const { + func.template operator()::type, + typename std::tuple_element<0, T_typelist>::type, + typename std::tuple_element<0, T_typelist>::type> + (param); + } +}; + + +/* + * apply_template_permutations uses template recursion to call the functor func + * with all 3-element permutations of the types in T_typelist + * + * Roughly, this corresponds to: + * for (i in (I - 1):0) // indexes go backwards, like: I - 1, I - 2, ... 0 + * for (j in (J - 1):0) + * for (k in (K - 1):0) + * func(param) + * + * tparam T_typelist Tuple templated with list of types to iterate through + * tparam T_functor Templated functor + * tparam T_param Parameter to pass to functor + */ +template +void apply_template_permutations(const T_functor& func, const T_param& param) { + apply_template_permutations_helper + ::value - 1, + std::tuple_size::value - 1, + std::tuple_size::value - 1>{}(func, param); +} + +#endif From 8e83e4f2c98859faa85c8d58ee412d8493c6eb35 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Wed, 25 Oct 2017 20:42:51 -0700 Subject: [PATCH 06/12] Cleaned up spacing, changed docs to make interface clearer (Issue 45) --- stan/math/prim/scal/prob/cauchy_rng.hpp | 23 +++--- .../prim/scal/prob/double_exponential_rng.hpp | 38 +++++----- .../prim/scal/prob/exp_mod_normal_rng.hpp | 32 ++++----- stan/math/prim/scal/prob/gumbel_rng.hpp | 32 ++++----- stan/math/prim/scal/prob/logistic_rng.hpp | 28 +++----- stan/math/prim/scal/prob/normal_rng.hpp | 27 +++---- stan/math/prim/scal/prob/skew_normal_rng.hpp | 38 +++++----- stan/math/prim/scal/prob/student_t_rng.hpp | 24 +++---- .../math/prim/mat/prob/VectorRNGTestRig.hpp | 72 +++++++------------ 9 files changed, 125 insertions(+), 189 deletions(-) diff --git a/stan/math/prim/scal/prob/cauchy_rng.hpp b/stan/math/prim/scal/prob/cauchy_rng.hpp index fe31a7fe065..8d68ea42b00 100644 --- a/stan/math/prim/scal/prob/cauchy_rng.hpp +++ b/stan/math/prim/scal/prob/cauchy_rng.hpp @@ -17,8 +17,8 @@ namespace stan { * Return a pseudorandom Cauchy variate for the given location and scale * using the specified random number generator. * - * mu and sigma can each be either scalars or vectors. All vector inputs - * must be the same length. + * mu and sigma can each be a scalar, a std::vector, an Eigen::Vector, or + * an Eigen::RowVector. Any non-scalar inputs must be the same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -28,32 +28,27 @@ namespace stan { * @param rng random number generator * @return Cauchy random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - cauchy_rng(const T_loc& mu, - const T_scale& sigma, - RNG& rng) { + cauchy_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::cauchy_distribution; - static const std::string function = "cauchy_rng"; scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); - check_consistent_sizes(function, - "Location parameter", mu, + check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); - size_t N = max_size(mu, sigma); - - VectorBuilder output(N); - - for (size_t n = 0; n < N; n++) { + for (size_t n = 0; n < N; ++n) { variate_generator > cauchy_rng(rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n])); output[n] = cauchy_rng(); diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 32b4d72a447..7d47766f784 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -8,7 +8,7 @@ #include #include #include -#include +#include #include #include @@ -18,8 +18,8 @@ namespace stan { * Return a pseudorandom double exponential variate with the given location * and scale using the specified random number generator. * - * mu and sigma can each be either scalars or vectors. All vector inputs - * must be the same length. + * mu and sigma can each be a scalar, a std::vector, an Eigen::Vector, or + * an Eigen::RowVector. Any non-scalar inputs must be the same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -29,38 +29,32 @@ namespace stan { * @param rng random number generator * @return double exponential random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - double_exponential_rng(const T_loc& mu, - const T_scale& sigma, - RNG& rng) { - static const std::string function = "double_exponential_rng"; - + double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; - using boost::random::uniform_01; + using boost::random::uniform_real_distribution; using std::abs; + static const std::string function = "double_exponential_rng"; scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); - check_consistent_sizes(function, - "Location parameter", mu, + check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); - size_t N = max_size(mu, sigma); - - VectorBuilder output(N); - - variate_generator > rng_unit_01(rng, uniform_01<>()); - for (size_t n = 0; n < N; n++) { - double laplaceRN = rng_unit_01(); - double a = (0.5 - laplaceRN > 0) ? 1.0 : -1.0; - output[n] = mu_vec[n] - - sigma_vec[n] * a * log1m(2 * abs(0.5 - laplaceRN)); + variate_generator > + z_rng(rng, uniform_real_distribution<>(-1.0, 1.0)); + for (size_t n = 0; n < N; ++n) { + double z = z_rng(); + output[n] = mu_vec[n] - sign(z) * sigma_vec[n] * log(abs(z)); } return output.data(); diff --git a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp index 99be01bd1d0..4653dfc6203 100644 --- a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp +++ b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp @@ -15,12 +15,13 @@ namespace stan { namespace math { /** - * Return a pseudorandom Exponentially modified normal variate for the + * Return a pseudorandom exponentially modified normal variate for the * given location, scale, and inverse scale using the specified random * number generator. * - * mu, sigma, and lambda can each be either scalars or vectors. All vector - * inputs must be the same length. + * mu, sigma, and lambda can each be a scalar, a std::vector, an + * Eigen::Vector, or an Eigen::RowVector. Any non-scalar inputs must be the + * same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -33,37 +34,32 @@ namespace stan { * @return Exponentially modified normal random variate * @throw std::domain_error if mu is infinite, sigma is nonpositive, * or lambda is nonpositive - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - exp_mod_normal_rng(const T_loc& mu, - const T_scale& sigma, - const T_inv_scale& lambda, - RNG& rng) { + exp_mod_normal_rng(const T_loc& mu, const T_scale& sigma, + const T_inv_scale& lambda, RNG& rng) { static const std::string function = "exp_mod_normal_rng"; scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); scalar_seq_view lambda_vec(lambda); + size_t N = max_size(mu, sigma, lambda); + VectorBuilder output(N); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_positive_finite(function, "Inv_scale parameter", lambda); - check_consistent_sizes(function, - "Location parameter", mu, + check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma, "Inv_scale Parameter", lambda); - size_t N = max_size(mu, sigma, lambda); - - VectorBuilder output(N); - - for (size_t n = 0; n < N; n++) { - output[n] = normal_rng(mu_vec[n], sigma_vec[n], rng) + - exponential_rng(lambda_vec[n], rng); - } + for (size_t n = 0; n < N; ++n) + output[n] = normal_rng(mu_vec[n], sigma_vec[n], rng) + + exponential_rng(lambda_vec[n], rng); return output.data(); } diff --git a/stan/math/prim/scal/prob/gumbel_rng.hpp b/stan/math/prim/scal/prob/gumbel_rng.hpp index 1e1596ca7c5..be0bf27f0e5 100644 --- a/stan/math/prim/scal/prob/gumbel_rng.hpp +++ b/stan/math/prim/scal/prob/gumbel_rng.hpp @@ -17,8 +17,8 @@ namespace stan { * Return a pseudorandom Gumbel variate with the given location and scale * using the specified random number generator. * - * mu and sigma can each be either scalars or vectors. All vector inputs - * must be the same length. + * mu and beta can each be a scalar, a std::vector, an Eigen::Vector, or + * an Eigen::RowVector. Any non-scalar inputs must be the same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -28,36 +28,30 @@ namespace stan { * @param rng random number generator * @return Gumbel random variate * @throw std::domain_error if mu is infinite or beta is nonpositive. - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - gumbel_rng(const T_loc& mu, - const T_scale& beta, - RNG& rng) { + gumbel_rng(const T_loc& mu, const T_scale& beta, RNG& rng) { using boost::variate_generator; using boost::uniform_01; - static const std::string function = "gumbel_rng"; - check_finite(function, "Location parameter", mu); - check_positive_finite(function, "Scale parameter", beta); - check_consistent_sizes(function, - "Location parameter", mu, - "Scale Parameter", beta); - scalar_seq_view mu_vec(mu); scalar_seq_view beta_vec(beta); - size_t N = max_size(mu, beta); - VectorBuilder output(N); + check_finite(function, "Location parameter", mu); + check_positive_finite(function, "Scale parameter", beta); + check_consistent_sizes(function, "Location parameter", mu, + "Scale Parameter", beta); + variate_generator > uniform01_rng(rng, uniform_01<>()); - for (size_t n = 0; n < N; n++) { - output[n] = mu_vec[n] - - beta_vec[n] * std::log(-std::log(uniform01_rng())); - } + for (size_t n = 0; n < N; ++n) + output[n] = mu_vec[n] + - beta_vec[n] * std::log(-std::log(uniform01_rng())); return output.data(); } diff --git a/stan/math/prim/scal/prob/logistic_rng.hpp b/stan/math/prim/scal/prob/logistic_rng.hpp index 7273a48d49d..3e2e5aeb556 100644 --- a/stan/math/prim/scal/prob/logistic_rng.hpp +++ b/stan/math/prim/scal/prob/logistic_rng.hpp @@ -17,8 +17,8 @@ namespace stan { * Return a pseudorandom Logistic variate for the given location and scale * using the specified random number generator. * - * mu and sigma can each be either scalars or vectors. All vector inputs - * must be the same length. + * mu and sigma can each be a scalar, a std::vector, an Eigen::Vector, or + * an Eigen::RowVector. Any non-scalar inputs must be the same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -28,36 +28,30 @@ namespace stan { * @param rng random number generator * @return Logistic random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - logistic_rng(const T_loc& mu, - const T_scale& sigma, - RNG& rng) { + logistic_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::exponential_distribution; - static const std::string function = "logistic_rng"; - check_finite(function, "Location parameter", mu); - check_positive_finite(function, "Scale parameter", sigma); - check_consistent_sizes(function, - "Location parameter", mu, - "Scale Parameter", sigma); - scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); + check_finite(function, "Location parameter", mu); + check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, "Location parameter", mu, + "Scale Parameter", sigma); + variate_generator > exp_rng(rng, exponential_distribution<>(1)); - for (size_t n = 0; n < N; n++) { + for (size_t n = 0; n < N; ++n) output[n] = mu_vec[n] - sigma_vec[n] * std::log(exp_rng() / exp_rng()); - } return output.data(); } diff --git a/stan/math/prim/scal/prob/normal_rng.hpp b/stan/math/prim/scal/prob/normal_rng.hpp index bf054f1459d..f7540827226 100644 --- a/stan/math/prim/scal/prob/normal_rng.hpp +++ b/stan/math/prim/scal/prob/normal_rng.hpp @@ -17,8 +17,8 @@ namespace stan { * Return a pseudorandom Normal variate for the given location and scale * using the specified random number generator. * - * mu and sigma can each be either a scalar or vector type (the ones that - * scalar_seq_view knows about). Any vector inputs must be the same length. + * mu and sigma can each be a scalar, a std::vector, an Eigen::Vector, or + * an Eigen::RowVector. Any non-scalar inputs must be the same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -28,32 +28,27 @@ namespace stan { * @param rng random number generator * @return Normal random variate * @throw std::domain_error if mu is infinite or sigma is nonpositive - * @throw std::invalid_argument if mu and sigma are different lengths + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - normal_rng(const T_loc& mu, - const T_scale& sigma, - RNG& rng) { + normal_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::normal_distribution; - static const std::string function = "normal_rng"; - check_finite(function, "Location parameter", mu); - check_positive_finite(function, "Scale parameter", sigma); - check_consistent_sizes(function, - "Location parameter", mu, - "Scale Parameter", sigma); - scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); - for (size_t n = 0; n < N; n++) { + check_finite(function, "Location parameter", mu); + check_positive_finite(function, "Scale parameter", sigma); + check_consistent_sizes(function, "Location parameter", mu, + "Scale Parameter", sigma); + + for (size_t n = 0; n < N; ++n) { variate_generator > norm_rng(rng, normal_distribution<>(mu_vec[n], sigma_vec[n])); output[n] = norm_rng(); diff --git a/stan/math/prim/scal/prob/skew_normal_rng.hpp b/stan/math/prim/scal/prob/skew_normal_rng.hpp index a2c24bf5504..5fff3f3519c 100644 --- a/stan/math/prim/scal/prob/skew_normal_rng.hpp +++ b/stan/math/prim/scal/prob/skew_normal_rng.hpp @@ -17,8 +17,9 @@ namespace stan { * Return a pseudorandom Skew-normal variate for the given location, scale, * and shape using the specified random number generator. * - * mu, sigma, and alpha can each be either scalars or vectors. All vector - * inputs must be the same length. + * mu, sigma, and alpha can each be a scalar, a std::vector, an + * Eigen::Vector, or an Eigen::RowVector. Any non-scalar inputs must be the + * same length. * * @tparam T_loc Type of location parameter * @tparam T_scale Type of scale parameter @@ -31,40 +32,35 @@ namespace stan { * @return Skew-normal random variate * @throw std::domain_error if mu is infinite, sigma is nonpositive, or * alpha is infinite - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - skew_normal_rng(const T_loc& mu, - const T_scale& sigma, - const T_shape& alpha, + skew_normal_rng(const T_loc& mu, const T_scale& sigma, const T_shape& alpha, RNG& rng) { using boost::variate_generator; using boost::random::normal_distribution; - static const std::string function = "skew_normal_rng"; - check_finite(function, "Location parameter", mu); - check_positive_finite(function, "Scale parameter", sigma); - check_finite(function, "Shape parameter", alpha); - check_consistent_sizes(function, - "Location parameter", mu, - "Scale Parameter", sigma, - "Shape Parameter", alpha); - scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); scalar_seq_view alpha_vec(alpha); - size_t N = max_size(mu, sigma, alpha); - VectorBuilder output(N); + check_finite(function, "Location parameter", mu); + check_positive_finite(function, "Scale parameter", sigma); + check_finite(function, "Shape parameter", alpha); + check_consistent_sizes(function, "Location parameter", mu, + "Scale Parameter", sigma, + "Shape Parameter", alpha); + variate_generator > - norm_rng(rng, normal_distribution<>(0.0, 1.0)); - for (size_t n = 0; n < N; n++) { - double r1 = norm_rng(), - r2 = norm_rng(); + norm_rng(rng, normal_distribution<>(0, 1)); + for (size_t n = 0; n < N; ++n) { + double r1 = norm_rng(); + double r2 = norm_rng(); if (r2 > alpha_vec[n] * r1) r1 = -r1; diff --git a/stan/math/prim/scal/prob/student_t_rng.hpp b/stan/math/prim/scal/prob/student_t_rng.hpp index ed2d4543d5b..36f86093d71 100644 --- a/stan/math/prim/scal/prob/student_t_rng.hpp +++ b/stan/math/prim/scal/prob/student_t_rng.hpp @@ -17,8 +17,9 @@ namespace stan { * Return a pseudorandom student-t variate for the given degrees of freedom, * location, and scale using the specified random number generator. * - * nu, mu, and sigma can each be either scalars or vectors. All vector - * inputs must be the same length. + * nu, mu, and sigma can each be a scalar, a std::vector, an + * Eigen::Vector, or an Eigen::RowVector. Any non-scalar inputs must be the + * same length. * * @tparam T_deg Type of degrees of freedom parameter * @tparam T_loc Type of location parameter @@ -31,36 +32,31 @@ namespace stan { * @return Student-t random variate * @throw std::domain_error if nu is nonpositive, mu is infinite, or sigma * is nonpositive - * @throw std::invalid_argument if vector arguments are not the same length + * @throw std::invalid_argument if non-scalars arguments are of different + * lengths */ template inline typename VectorBuilder::type - student_t_rng(const T_deg& nu, - const T_loc& mu, - const T_scale& sigma, + student_t_rng(const T_deg& nu, const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::student_t_distribution; - static const std::string function = "student_t_rng"; scalar_seq_view nu_vec(nu); scalar_seq_view mu_vec(mu); scalar_seq_view sigma_vec(sigma); + size_t N = max_size(nu, mu, sigma); + VectorBuilder output(N); check_positive_finite(function, "Degrees of freedom parameter", nu); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); - check_consistent_sizes(function, - "Degrees of freedom parameter", nu, + check_consistent_sizes(function, "Degrees of freedom parameter", nu, "Location parameter", mu, "Scale Parameter", sigma); - size_t N = max_size(nu, mu, sigma); - - VectorBuilder output(N); - - for (size_t n = 0; n < N; n++) { + for (size_t n = 0; n < N; ++n) { variate_generator > rng_unit_student_t(rng, student_t_distribution<>(nu_vec[n])); output[n] = mu_vec[n] + sigma_vec[n] * rng_unit_student_t(); diff --git a/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp index 4d8e0317b20..bdfcf4ebffe 100644 --- a/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp +++ b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp @@ -110,60 +110,36 @@ class VectorRNGTestRig { double p3) const = 0; VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector good_p1_int, - std::vector bad_p1, - std::vector bad_p1_int, - std::vector good_p2, - std::vector good_p2_int, - std::vector bad_p2, - std::vector bad_p2_int, - std::vector good_p3, - std::vector good_p3_int, - std::vector bad_p3, - std::vector bad_p3_int) : N_(N), M_(M), - good_p1_(good_p1), - good_p1_int_(good_p1_int), - bad_p1_(bad_p1), - bad_p1_int_(bad_p1_int), - good_p2_(good_p2), - good_p2_int_(good_p2_int), - bad_p2_(bad_p2), - bad_p2_int_(bad_p2_int), - good_p3_(good_p3), - good_p3_int_(good_p3_int), - bad_p3_(bad_p3), - bad_p3_int_(bad_p3_int) { + std::vector good_p1, std::vector good_p1_int, + std::vector bad_p1, std::vector bad_p1_int, + std::vector good_p2, std::vector good_p2_int, + std::vector bad_p2, std::vector bad_p2_int, + std::vector good_p3, std::vector good_p3_int, + std::vector bad_p3, std::vector bad_p3_int) + : N_(N), M_(M), good_p1_(good_p1), good_p1_int_(good_p1_int), + bad_p1_(bad_p1), bad_p1_int_(bad_p1_int), + good_p2_(good_p2), good_p2_int_(good_p2_int), + bad_p2_(bad_p2), bad_p2_int_(bad_p2_int), + good_p3_(good_p3), good_p3_int_(good_p3_int), + bad_p3_(bad_p3), bad_p3_int_(bad_p3_int) { } VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector good_p1_int, - std::vector bad_p1, - std::vector bad_p1_int, - std::vector good_p2, - std::vector good_p2_int, - std::vector bad_p2, - std::vector bad_p2_int) : N_(N), M_(M), - good_p1_(good_p1), - good_p1_int_(good_p1_int), - bad_p1_(bad_p1), - bad_p1_int_(bad_p1_int), - good_p2_(good_p2), - good_p2_int_(good_p2_int), - bad_p2_(bad_p2), - bad_p2_int_(bad_p2_int) { + std::vector good_p1, std::vector good_p1_int, + std::vector bad_p1, std::vector bad_p1_int, + std::vector good_p2, std::vector good_p2_int, + std::vector bad_p2, std::vector bad_p2_int) + : N_(N), M_(M), good_p1_(good_p1), good_p1_int_(good_p1_int), + bad_p1_(bad_p1), bad_p1_int_(bad_p1_int), + good_p2_(good_p2), good_p2_int_(good_p2_int), + bad_p2_(bad_p2), bad_p2_int_(bad_p2_int) { } VectorRNGTestRig(int N, int M, - std::vector good_p1, - std::vector good_p1_int, - std::vector bad_p1, - std::vector bad_p1_int) : N_(N), M_(M), - good_p1_(good_p1), - good_p1_int_(good_p1_int), - bad_p1_(bad_p1), - bad_p1_int_(bad_p1_int) { + std::vector good_p1, std::vector good_p1_int, + std::vector bad_p1, std::vector bad_p1_int) + : N_(N), M_(M), good_p1_(good_p1), good_p1_int_(good_p1_int), + bad_p1_(bad_p1), bad_p1_int_(bad_p1_int) { } /* From ee338ffc6ce522ccea96f6dcff003c62c3cd4ae9 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Wed, 25 Oct 2017 21:23:25 -0700 Subject: [PATCH 07/12] Fixed cpplint issues in test (Issue 45) --- .../math/prim/mat/prob/VectorRNGTestRig.hpp | 8 +++--- test/unit/math/prim/mat/prob/cauchy_test.cpp | 6 ++-- .../prim/mat/prob/double_exponential_test.cpp | 6 ++-- .../prim/mat/prob/exp_mod_normal_test.cpp | 1 + test/unit/math/prim/mat/prob/gumbel_test.cpp | 6 ++-- .../unit/math/prim/mat/prob/logistic_test.cpp | 6 ++-- test/unit/math/prim/mat/prob/normal_test.cpp | 12 ++++---- .../math/prim/mat/prob/skew_normal_test.cpp | 6 ++-- .../math/prim/mat/prob/student_t_test.cpp | 8 +++--- .../prim/mat/prob/vector_rng_test_helper.hpp | 24 ++++++++++------ .../scal/meta/apply_template_permutations.hpp | 28 +++++++++++-------- 11 files changed, 62 insertions(+), 49 deletions(-) diff --git a/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp index bdfcf4ebffe..cc0bbdfa59d 100644 --- a/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp +++ b/test/unit/math/prim/mat/prob/VectorRNGTestRig.hpp @@ -69,8 +69,8 @@ */ class VectorRNGTestRig { public: - int N_; // Number of samples used in the quantiles tests - int M_; // Length of vectors for the vectorization tests + int N_; // Number of samples used in the quantiles tests + int M_; // Length of vectors for the vectorization tests std::vector good_p1_; std::vector good_p1_int_; @@ -84,7 +84,7 @@ class VectorRNGTestRig { std::vector good_p3_int_; std::vector bad_p3_; std::vector bad_p3_int_; - + std::vector always_bad_values_ = { stan::math::positive_infinity(), stan::math::negative_infinity(), stan::math::not_a_number() }; @@ -134,7 +134,7 @@ class VectorRNGTestRig { good_p2_(good_p2), good_p2_int_(good_p2_int), bad_p2_(bad_p2), bad_p2_int_(bad_p2_int) { } - + VectorRNGTestRig(int N, int M, std::vector good_p1, std::vector good_p1_int, std::vector bad_p1, std::vector bad_p1_int) diff --git a/test/unit/math/prim/mat/prob/cauchy_test.cpp b/test/unit/math/prim/mat/prob/cauchy_test.cpp index f239e24dd2a..909a8bf93ae 100644 --- a/test/unit/math/prim/mat/prob/cauchy_test.cpp +++ b/test/unit/math/prim/mat/prob/cauchy_test.cpp @@ -18,7 +18,7 @@ class CauchyTestRig : public VectorRNGTestRig { {1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0}, {-3, -2, -1, 0}) {} - + template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { @@ -30,13 +30,13 @@ class CauchyTestRig : public VectorRNGTestRig { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::cauchy_distribution<> dist(mu, sigma); - + for (int i = 1; i < K; ++i) { double frac = i / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/double_exponential_test.cpp b/test/unit/math/prim/mat/prob/double_exponential_test.cpp index f7f9ecf8514..5b4d42bf805 100644 --- a/test/unit/math/prim/mat/prob/double_exponential_test.cpp +++ b/test/unit/math/prim/mat/prob/double_exponential_test.cpp @@ -24,19 +24,19 @@ class DoubleExponentialTestRig : public VectorRNGTestRig { T_rng& rng) const { return stan::math::double_exponential_rng(mu, sigma, rng); } - + std::vector generate_quantiles(double mu, double sigma, double unused) const { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::laplace_distribution<> dist(mu, sigma); - + for (int i = 1; i < K; ++i) { double frac = i / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp index d6949481690..376ff0d893b 100644 --- a/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/exp_mod_normal_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include class ExpModNormalTestRig : public VectorRNGTestRig { public: diff --git a/test/unit/math/prim/mat/prob/gumbel_test.cpp b/test/unit/math/prim/mat/prob/gumbel_test.cpp index d24ea894fb5..5196441d768 100644 --- a/test/unit/math/prim/mat/prob/gumbel_test.cpp +++ b/test/unit/math/prim/mat/prob/gumbel_test.cpp @@ -18,7 +18,7 @@ class GumbelTestRig : public VectorRNGTestRig { {1, 2, 3, 4}, {-1.0, -1.5, -2.5, -0.7, 0.0}, {-1, -2, -3, -4, 0}) {} - + template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { @@ -30,13 +30,13 @@ class GumbelTestRig : public VectorRNGTestRig { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::extreme_value_distribution<> dist(mu, sigma); - + for (int i = 1; i < K; ++i) { double frac = static_cast(i) / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/logistic_test.cpp b/test/unit/math/prim/mat/prob/logistic_test.cpp index 7c97100e7a1..58aa7ca56ee 100644 --- a/test/unit/math/prim/mat/prob/logistic_test.cpp +++ b/test/unit/math/prim/mat/prob/logistic_test.cpp @@ -18,7 +18,7 @@ class LogisticTestRig : public VectorRNGTestRig { {1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0}, {-3, -2, -1, 0}) {} - + template auto generate_samples(const T1& mu, const T2& sigma, const T3& unused, T_rng& rng) const { @@ -30,13 +30,13 @@ class LogisticTestRig : public VectorRNGTestRig { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::logistic_distribution<> dist(mu, sigma); - + for (int i = 1; i < K; ++i) { double frac = i / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/normal_test.cpp b/test/unit/math/prim/mat/prob/normal_test.cpp index 4a5213946ae..881bf21acd6 100644 --- a/test/unit/math/prim/mat/prob/normal_test.cpp +++ b/test/unit/math/prim/mat/prob/normal_test.cpp @@ -15,10 +15,10 @@ class NormalTestRig : public VectorRNGTestRig { * arguments. */ NormalTestRig() : - VectorRNGTestRig(10000, // Number of samples used for quantiles tests - 10, // Length of vectors for vectorization tests - {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, // Valid values for p1 - {-3, -2, -1, 0, 2, 6}, // Valid integer values for p1 + VectorRNGTestRig(10000, // Number of samples used for quantiles tests + 10, // Length of vectors for vectorization tests + {-2.5, -1.7, -0.1, 0.0, 2.0, 5.8}, // Valid values for p1 + {-3, -2, -1, 0, 2, 6}, // Valid integer values for p1 {}, {}, {0.1, 1.0, 2.5, 4.0}, @@ -47,13 +47,13 @@ class NormalTestRig : public VectorRNGTestRig { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::normal_distribution<> dist(mu, sigma); - + for (int i = 1; i < K; ++i) { double frac = i / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/skew_normal_test.cpp b/test/unit/math/prim/mat/prob/skew_normal_test.cpp index c3615c15017..747494ccc4c 100644 --- a/test/unit/math/prim/mat/prob/skew_normal_test.cpp +++ b/test/unit/math/prim/mat/prob/skew_normal_test.cpp @@ -23,7 +23,7 @@ class SkewNormalTestRig : public VectorRNGTestRig { {-2, -1, 0, 1, 2}, {}, {}) {} - + template auto generate_samples(const T1& mu, const T2& sigma, const T3& alpha, T_rng& rng) const { @@ -35,13 +35,13 @@ class SkewNormalTestRig : public VectorRNGTestRig { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); boost::math::skew_normal_distribution<> dist(mu, sigma, alpha); - + for (int i = 1; i < K; ++i) { double frac = i / K; quantiles.push_back(quantile(dist, frac)); } quantiles.push_back(std::numeric_limits::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/student_t_test.cpp b/test/unit/math/prim/mat/prob/student_t_test.cpp index ba0ee953fd3..dd186427a1b 100644 --- a/test/unit/math/prim/mat/prob/student_t_test.cpp +++ b/test/unit/math/prim/mat/prob/student_t_test.cpp @@ -22,7 +22,7 @@ class StudentTTestRig : public VectorRNGTestRig { {1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0}, {-3, -2, -1, 0}) {} - + template auto generate_samples(const T1& nu, const T2& mu, const T3& sigma, T_rng& rng) const { @@ -33,14 +33,14 @@ class StudentTTestRig : public VectorRNGTestRig { const { std::vector quantiles; double K = boost::math::round(2 * std::pow(N_, 0.4)); - + boost::math::students_t_distribution<>dist(nu); - for (int i=1; i::max()); - + return quantiles; } }; diff --git a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp index df01cfbc124..9ef6ed33ef6 100644 --- a/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp +++ b/test/unit/math/prim/mat/prob/vector_rng_test_helper.hpp @@ -8,6 +8,8 @@ #include #include #include +#include +#include #include /* @@ -146,7 +148,8 @@ void resize_if_vector(int& v, int N) { * @param T_rig Test rig for random number generator */ struct check_dist_throws { - template + template void operator()(const T_rig& rig) const { boost::random::mt19937 rng; @@ -200,7 +203,7 @@ struct check_dist_throws { } // Now try putting incompatible values in second parameter - if(p2_is_used) { + if (p2_is_used) { for (auto bad_p2_value : bad_p2) { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, { bad_p2_value }); @@ -210,7 +213,7 @@ struct check_dist_throws { } // Now try putting incompatible values in third parameter - if(p3_is_used) { + if (p3_is_used) { for (auto bad_p3_value : bad_p3) { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); @@ -222,7 +225,7 @@ struct check_dist_throws { // Check to make sure _rng rejects vector arguments of different lengths for // all parameter pairs. - // If p1 is a scalar or the only vector, this test is skipped + // If p1 is a scalar or the only vector, this test is skipped resize_if_vector(p1, 3); // No-op if p1 is a scalar resize_if_vector(p2, 4); // No-op if p2 is a scalar resize_if_vector(p3, 4); // No-op if p3 is a scalar @@ -232,7 +235,8 @@ struct check_dist_throws { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); assign_parameter_values(p3, good_p3); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), + std::invalid_argument); } // If p2 is a scalar or the only vector, this test is skipped @@ -244,7 +248,8 @@ struct check_dist_throws { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); assign_parameter_values(p3, good_p3); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), + std::invalid_argument); } // If p3 is a scalar or the only vector, this test is skipped @@ -256,7 +261,8 @@ struct check_dist_throws { assign_parameter_values(p1, good_p1); assign_parameter_values(p2, good_p2); assign_parameter_values(p3, good_p3); - EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), std::invalid_argument); + EXPECT_THROW(rig.generate_samples(p1, p2, p3, rng), + std::invalid_argument); } } }; @@ -347,8 +353,8 @@ struct check_quantiles { std::vector > samples_to_test_transpose; for (int n = 0; n < rig.N_; ++n) { - // If p1, p2, and p3 are scalars, the output is a scalar. Need to promote it - // to a std::vector + // If p1, p2, and p3 are scalars, the output is a scalar. Need to promote + // it to a std::vector samples_to_test_transpose. push_back(promote_to_vector(rig.generate_samples(p1, p2, p3, rng))); } diff --git a/test/unit/math/prim/scal/meta/apply_template_permutations.hpp b/test/unit/math/prim/scal/meta/apply_template_permutations.hpp index 19e8d911d73..709ecbb97c9 100644 --- a/test/unit/math/prim/scal/meta/apply_template_permutations.hpp +++ b/test/unit/math/prim/scal/meta/apply_template_permutations.hpp @@ -15,7 +15,8 @@ * tparam J Loop index for second type * tparam K Loop index for third type */ -template +template struct apply_template_permutations_helper { void operator()(const T_functor& func, const T_param& param) const { func.template operator()::type, @@ -23,8 +24,8 @@ struct apply_template_permutations_helper { typename std::tuple_element::type> (param); - apply_template_permutations_helper{}(func, - param); + apply_template_permutations_helper + {}(func, param); } }; @@ -40,16 +41,19 @@ struct apply_template_permutations_helper { * tparam J Loop index for second type * tparam K Loop index for third type */ -template -struct apply_template_permutations_helper { +template +struct apply_template_permutations_helper { void operator()(const T_functor& func, const T_param& param) const { func.template operator()::type, typename std::tuple_element::type, typename std::tuple_element<0, T_typelist>::type> (param); - - apply_template_permutations_helper::value - 1>{}(func, param); + + apply_template_permutations_helper + ::value - 1>{}(func, param); } }; @@ -64,13 +68,14 @@ struct apply_template_permutations_helper -struct apply_template_permutations_helper { +struct apply_template_permutations_helper { void operator()(const T_functor& func, const T_param& param) const { func.template operator()::type, typename std::tuple_element<0, T_typelist>::type, typename std::tuple_element<0, T_typelist>::type> (param); - + apply_template_permutations_helper::value - 1, std::tuple_size::value - 1>{}(func, @@ -88,7 +93,8 @@ struct apply_template_permutations_helper -struct apply_template_permutations_helper { +struct apply_template_permutations_helper { void operator()(const T_functor& func, const T_param& param) const { func.template operator()::type, typename std::tuple_element<0, T_typelist>::type, From 9bf8a0ea2d0dd227c1d2c2e66cce728d32cce725 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Fri, 27 Oct 2017 18:24:11 -0700 Subject: [PATCH 08/12] Swapped std::log for log in double_exponential test to get things working --- stan/math/prim/scal/prob/double_exponential_rng.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 7d47766f784..eb215fdc8ca 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -54,7 +54,7 @@ namespace stan { z_rng(rng, uniform_real_distribution<>(-1.0, 1.0)); for (size_t n = 0; n < N; ++n) { double z = z_rng(); - output[n] = mu_vec[n] - sign(z) * sigma_vec[n] * log(abs(z)); + output[n] = mu_vec[n] - sign(z) * sigma_vec[n] * std::log(abs(z)); } return output.data(); From aa207eba210a0458d7b5f8e39be2830b917e05d7 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 28 Oct 2017 10:54:57 -0700 Subject: [PATCH 09/12] Made a couple fixes to get Jenkins tests to pass (errors not reproducing locally) (Issue 45) --- stan/math/prim/scal/prob/double_exponential_rng.hpp | 2 +- test/unit/math/prim/scal/prob/util.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index eb215fdc8ca..7e84a4e0442 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -54,7 +54,7 @@ namespace stan { z_rng(rng, uniform_real_distribution<>(-1.0, 1.0)); for (size_t n = 0; n < N; ++n) { double z = z_rng(); - output[n] = mu_vec[n] - sign(z) * sigma_vec[n] * std::log(abs(z)); + output[n] = mu_vec[n] - stan::math::sign(z) * sigma_vec[n] * std::log(abs(z)); } return output.data(); diff --git a/test/unit/math/prim/scal/prob/util.hpp b/test/unit/math/prim/scal/prob/util.hpp index e25b9475d17..9e6c1b66c02 100644 --- a/test/unit/math/prim/scal/prob/util.hpp +++ b/test/unit/math/prim/scal/prob/util.hpp @@ -47,7 +47,7 @@ void assert_matches_quantiles(const std::vector& samples, } std::vector counts(K); - int current_index = 0; + size_t current_index = 0; for (int i = 0; i < N; ++i) { while (mysamples[i] >= quantiles[current_index]) { ++current_index; From 45176680600788d02b90fdb72dd930b75f8fada7 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 28 Oct 2017 12:51:34 -0700 Subject: [PATCH 10/12] Fix cpplint (Issue 45) --- stan/math/prim/scal/prob/double_exponential_rng.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 7e84a4e0442..31e6dbe0956 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -54,7 +54,8 @@ namespace stan { z_rng(rng, uniform_real_distribution<>(-1.0, 1.0)); for (size_t n = 0; n < N; ++n) { double z = z_rng(); - output[n] = mu_vec[n] - stan::math::sign(z) * sigma_vec[n] * std::log(abs(z)); + output[n] = mu_vec[n] + - stan::math::sign(z) * sigma_vec[n] * std::log(abs(z)); } return output.data(); From eff2de964f6ccacae9383a4dfad01e4caa4d741e Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 28 Oct 2017 17:07:33 -0700 Subject: [PATCH 11/12] Fixed references to std::log and std::abs in double_exponential_rng. Removed reference to stan::math::sign (Issue 45) --- stan/math/prim/scal/prob/double_exponential_rng.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index 31e6dbe0956..43a548ece1f 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -37,7 +37,6 @@ namespace stan { double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) { using boost::variate_generator; using boost::random::uniform_real_distribution; - using std::abs; static const std::string function = "double_exponential_rng"; scalar_seq_view mu_vec(mu); @@ -55,7 +54,7 @@ namespace stan { for (size_t n = 0; n < N; ++n) { double z = z_rng(); output[n] = mu_vec[n] - - stan::math::sign(z) * sigma_vec[n] * std::log(abs(z)); + - ((z > 0) ? 1.0 : -1.0) * sigma_vec[n] * std::log(std::abs(z)); } return output.data(); From 47cdfc9c4f70d20fc100239eb1bbb3b10ab260b1 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 13 Nov 2017 21:16:52 -0800 Subject: [PATCH 12/12] Moved around variable declarations below argument checks (Issue 45) --- stan/math/prim/scal/prob/cauchy_rng.hpp | 10 +++++----- stan/math/prim/scal/prob/double_exponential_rng.hpp | 10 +++++----- stan/math/prim/scal/prob/exp_mod_normal_rng.hpp | 12 ++++++------ stan/math/prim/scal/prob/gumbel_rng.hpp | 10 +++++----- stan/math/prim/scal/prob/logistic_rng.hpp | 10 +++++----- stan/math/prim/scal/prob/normal_rng.hpp | 10 +++++----- stan/math/prim/scal/prob/skew_normal_rng.hpp | 12 ++++++------ stan/math/prim/scal/prob/student_t_rng.hpp | 12 ++++++------ 8 files changed, 43 insertions(+), 43 deletions(-) diff --git a/stan/math/prim/scal/prob/cauchy_rng.hpp b/stan/math/prim/scal/prob/cauchy_rng.hpp index 707ddfff222..453d710b4a3 100644 --- a/stan/math/prim/scal/prob/cauchy_rng.hpp +++ b/stan/math/prim/scal/prob/cauchy_rng.hpp @@ -37,16 +37,16 @@ namespace stan { using boost::random::cauchy_distribution; static const char* function = "cauchy_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); + for (size_t n = 0; n < N; ++n) { variate_generator > cauchy_rng(rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n])); diff --git a/stan/math/prim/scal/prob/double_exponential_rng.hpp b/stan/math/prim/scal/prob/double_exponential_rng.hpp index ae6b5b1b28d..3c7f54942e3 100644 --- a/stan/math/prim/scal/prob/double_exponential_rng.hpp +++ b/stan/math/prim/scal/prob/double_exponential_rng.hpp @@ -38,16 +38,16 @@ namespace stan { using boost::random::uniform_real_distribution; static const char* function = "double_exponential_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); + variate_generator > z_rng(rng, uniform_real_distribution<>(-1.0, 1.0)); for (size_t n = 0; n < N; ++n) { diff --git a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp index 4d81a5d8a6f..95fceaee7fd 100644 --- a/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp +++ b/stan/math/prim/scal/prob/exp_mod_normal_rng.hpp @@ -43,12 +43,6 @@ namespace stan { const T_inv_scale& lambda, RNG& rng) { static const char* function = "exp_mod_normal_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - scalar_seq_view lambda_vec(lambda); - size_t N = max_size(mu, sigma, lambda); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_positive_finite(function, "Inv_scale parameter", lambda); @@ -56,6 +50,12 @@ namespace stan { "Scale Parameter", sigma, "Inv_scale Parameter", lambda); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + scalar_seq_view lambda_vec(lambda); + size_t N = max_size(mu, sigma, lambda); + VectorBuilder output(N); + for (size_t n = 0; n < N; ++n) output[n] = normal_rng(mu_vec[n], sigma_vec[n], rng) + exponential_rng(lambda_vec[n], rng); diff --git a/stan/math/prim/scal/prob/gumbel_rng.hpp b/stan/math/prim/scal/prob/gumbel_rng.hpp index 6180dc8365e..c5e9f809bf8 100644 --- a/stan/math/prim/scal/prob/gumbel_rng.hpp +++ b/stan/math/prim/scal/prob/gumbel_rng.hpp @@ -37,16 +37,16 @@ namespace stan { using boost::uniform_01; static const char* function = "gumbel_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view beta_vec(beta); - size_t N = max_size(mu, beta); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", beta); check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", beta); + scalar_seq_view mu_vec(mu); + scalar_seq_view beta_vec(beta); + size_t N = max_size(mu, beta); + VectorBuilder output(N); + variate_generator > uniform01_rng(rng, uniform_01<>()); for (size_t n = 0; n < N; ++n) output[n] = mu_vec[n] diff --git a/stan/math/prim/scal/prob/logistic_rng.hpp b/stan/math/prim/scal/prob/logistic_rng.hpp index 6b314c8f71f..e84d67fed69 100644 --- a/stan/math/prim/scal/prob/logistic_rng.hpp +++ b/stan/math/prim/scal/prob/logistic_rng.hpp @@ -37,16 +37,16 @@ namespace stan { using boost::random::exponential_distribution; static const char* function = "logistic_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); + variate_generator > exp_rng(rng, exponential_distribution<>(1)); for (size_t n = 0; n < N; ++n) diff --git a/stan/math/prim/scal/prob/normal_rng.hpp b/stan/math/prim/scal/prob/normal_rng.hpp index d11adb21a1e..50ae0b18f82 100644 --- a/stan/math/prim/scal/prob/normal_rng.hpp +++ b/stan/math/prim/scal/prob/normal_rng.hpp @@ -37,16 +37,16 @@ namespace stan { using boost::normal_distribution; static const char* function = "normal_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - size_t N = max_size(mu, sigma); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter", sigma); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + size_t N = max_size(mu, sigma); + VectorBuilder output(N); + for (size_t n = 0; n < N; ++n) { variate_generator > norm_rng(rng, normal_distribution<>(mu_vec[n], sigma_vec[n])); diff --git a/stan/math/prim/scal/prob/skew_normal_rng.hpp b/stan/math/prim/scal/prob/skew_normal_rng.hpp index 38d6607ddc4..deb29a903be 100644 --- a/stan/math/prim/scal/prob/skew_normal_rng.hpp +++ b/stan/math/prim/scal/prob/skew_normal_rng.hpp @@ -42,12 +42,6 @@ namespace stan { using boost::random::normal_distribution; static const char* function = "skew_normal_rng"; - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - scalar_seq_view alpha_vec(alpha); - size_t N = max_size(mu, sigma, alpha); - VectorBuilder output(N); - check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); check_finite(function, "Shape parameter", alpha); @@ -55,6 +49,12 @@ namespace stan { "Scale Parameter", sigma, "Shape Parameter", alpha); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + scalar_seq_view alpha_vec(alpha); + size_t N = max_size(mu, sigma, alpha); + VectorBuilder output(N); + variate_generator > norm_rng(rng, normal_distribution<>(0, 1)); for (size_t n = 0; n < N; ++n) { diff --git a/stan/math/prim/scal/prob/student_t_rng.hpp b/stan/math/prim/scal/prob/student_t_rng.hpp index 6f4e603d002..e778fba939c 100644 --- a/stan/math/prim/scal/prob/student_t_rng.hpp +++ b/stan/math/prim/scal/prob/student_t_rng.hpp @@ -42,12 +42,6 @@ namespace stan { using boost::random::student_t_distribution; static const char* function = "student_t_rng"; - scalar_seq_view nu_vec(nu); - scalar_seq_view mu_vec(mu); - scalar_seq_view sigma_vec(sigma); - size_t N = max_size(nu, mu, sigma); - VectorBuilder output(N); - check_positive_finite(function, "Degrees of freedom parameter", nu); check_finite(function, "Location parameter", mu); check_positive_finite(function, "Scale parameter", sigma); @@ -55,6 +49,12 @@ namespace stan { "Location parameter", mu, "Scale Parameter", sigma); + scalar_seq_view nu_vec(nu); + scalar_seq_view mu_vec(mu); + scalar_seq_view sigma_vec(sigma); + size_t N = max_size(nu, mu, sigma); + VectorBuilder output(N); + for (size_t n = 0; n < N; ++n) { variate_generator > rng_unit_student_t(rng, student_t_distribution<>(nu_vec[n]));