From 0a68fdb1c243a45b779904a00729aa3bb239bd7d Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 14 Jun 2025 11:35:01 +0800 Subject: [PATCH 1/7] Initial version of gaussian copula lpdf --- stan/math/prim/fun/size.hpp | 6 + stan/math/prim/fun/size_mvt.hpp | 10 ++ stan/math/prim/fun/vector_seq_view.hpp | 7 +- stan/math/prim/prob.hpp | 1 + .../prob/gaussian_copula_cholesky_lpdf.hpp | 140 ++++++++++++++++++ 5 files changed, 161 insertions(+), 3 deletions(-) create mode 100644 stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp diff --git a/stan/math/prim/fun/size.hpp b/stan/math/prim/fun/size.hpp index a9005d97b39..8a5418ecc24 100644 --- a/stan/math/prim/fun/size.hpp +++ b/stan/math/prim/fun/size.hpp @@ -34,6 +34,12 @@ inline int64_t size(const T& m) { return m.size(); } + +template * = nullptr> +inline int64_t size(const T& /*t*/) { + return std::tuple_size>{}; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/fun/size_mvt.hpp b/stan/math/prim/fun/size_mvt.hpp index 26b55e145d8..d766b500b9d 100644 --- a/stan/math/prim/fun/size_mvt.hpp +++ b/stan/math/prim/fun/size_mvt.hpp @@ -36,11 +36,21 @@ int64_t size_mvt(const MatrixT& /* unused */) { return 1; } +template * = nullptr> +int64_t size_mvt(const TupleT& /* unused */) { + return 1; +} + template * = nullptr> int64_t size_mvt(const std::vector& x) { return x.size(); } +template * = nullptr> +int64_t size_mvt(const std::vector& x) { + return x.size(); +} + template * = nullptr> int64_t size_mvt(const std::vector& x) { return x.size(); diff --git a/stan/math/prim/fun/vector_seq_view.hpp b/stan/math/prim/fun/vector_seq_view.hpp index 8c76979927a..b11d8e291f3 100644 --- a/stan/math/prim/fun/vector_seq_view.hpp +++ b/stan/math/prim/fun/vector_seq_view.hpp @@ -8,12 +8,13 @@ namespace stan { namespace internal { template -using is_matrix_or_std_vector - = math::disjunction, is_std_vector>; +using is_matrix_or_std_vector_or_tuple + = math::disjunction, is_std_vector, math::is_tuple>; template using is_scalar_container = math::disjunction< is_matrix, + math::is_tuple, math::conjunction, is_stan_scalar>>>; } // namespace internal @@ -76,7 +77,7 @@ class vector_seq_view>> { */ template class vector_seq_view< - T, require_std_vector_vt> { + T, require_std_vector_vt> { public: explicit vector_seq_view(const T& v) noexcept : v_(v) {} inline auto size() const noexcept { return v_.size(); } diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index e04c7e87172..0bbc6b30391 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -111,6 +111,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp new file mode 100644 index 00000000000..78a1e7c9a37 --- /dev/null +++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp @@ -0,0 +1,140 @@ +#ifndef STAN_MATH_PRIM_PROB_GAUSSIAN_COPULA_CHOLESKY_LPDF_HPP +#define STAN_MATH_PRIM_PROB_GAUSSIAN_COPULA_CHOLESKY_LPDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +template +auto apply_cdfs_elem_impl(Ty_elem&& y_elem, Ttuple&& cdf_tuple, std::index_sequence){ + auto&& cdf_func = std::get<0>(cdf_tuple); + return cdf_func( + std::forward(y_elem), + std::get(std::forward(cdf_tuple))... + ); +} + + +template +auto apply_cdfs_impl(Ty&& y, Ttuple&& cdf_tuple, std::index_sequence){ + plain_type_t res(y.size()); + + // Use initializer-list expansion to assign the result of each CDF + // to the respective element in the result vector, as we need a constant + // expression for indexing the tuple of CDF functors + static_cast(std::initializer_list{( + res[I] = apply_cdfs_elem_impl( + std::forward(y[I]), + std::get(cdf_tuple), + std::make_index_sequence< + // Using size - 1 as the first element is the functor to apply + std::tuple_size< + std::remove_reference_t< + decltype(std::get(cdf_tuple))>>{} - 1>{}), + 0)...}); + + return res; +} + +template +auto apply_cdfs(Ty&& y, Ttuple&& cdf_tuple){ + return apply_cdfs_impl( + std::forward(y), + std::forward(cdf_tuple), + std::make_index_sequence< + std::tuple_size>{}>{} + ); +} +} + +/* + vectors ~ gaussian_copula_cholesky(cdf_funs_tuple, chol) + (cdf_fun, ...) +*/ + +template > +T_return gaussian_copula_cholesky_lpdf( + const T_y& y, const T_cdf_fun_tuple& cdf_fun_tuple, const T_chol chol) { + static constexpr const char* function = "gaussian_copula_cholesky_lpdf"; + + using T_y_ref = ref_type_t; + using T_chochol_ref = ref_type_t; + using T_cdf_ref = ref_type_t; + + check_consistent_sizes_mvt(function, "y", y, "cdf_fun_tuple", cdf_fun_tuple); + const size_t size_mvt_y = size_mvt(y); + const size_t size_mvt_cdf_tuple = size_mvt(cdf_fun_tuple); + if (size_mvt_y == 0) { + return 0; + } + T_y_ref y_ref = y; + T_chochol_ref chol_ref = chol; + T_cdf_ref cdf_tuple_ref = cdf_fun_tuple; + + vector_seq_view y_vec(y_ref); + vector_seq_view cdf_tuple_vec(cdf_tuple_ref); + const size_t size_vec = max_size_mvt(y, cdf_fun_tuple); + + const int size_y = math::size(y_vec[0]); + const int size_cdf_tuple = math::size(cdf_tuple_vec[0]); + + // check size consistency of all random variables y + for (size_t i = 1; i < size_mvt_y; i++) { + check_size_match(function, + "Size of one of the vectors of " + "the random variable", + math::size(y_vec[i]), + "Size of the first vector of the " + "random variable", + size_y); + } + // check size consistency of all CDF tuples + for (size_t i = 1; i < size_mvt_cdf_tuple; i++) { + check_size_match(function, + "Size of one of the vectors of " + "the CDF functor tuple", + math::size(cdf_tuple_vec[i]), + "Size of the first vector of the " + "location variable", + size_cdf_tuple); + } + + check_size_match(function, "Size of random variable", size_y, + "size of CDF functor tuple", size_cdf_tuple); + check_size_match(function, "Size of random variable", size_y, + "rows of covariance parameter", chol.rows()); + check_size_match(function, "Size of random variable", size_y, + "columns of covariance parameter", chol.cols()); + + for (size_t i = 0; i < size_vec; i++) { + check_not_nan(function, "Random variable", y_vec[i]); + } + check_cholesky_factor(function, "Cholesky decomposition of a variance matrix", + chol_ref); + + promote_scalar_t> q(size_vec); + T_return lp(0); + for (size_t i = 0; i < size_vec; i++) { + q[i] = to_vector(inv_Phi(internal::apply_cdfs(y_vec[i], cdf_tuple_vec[i]))); + lp -= std_normal_lpdf(q[i]); + } + const std::vector zero_vec(size_vec, rep_vector(0, size_y)); + lp += multi_normal_cholesky_lpdf(q, zero_vec, chol_ref); + return lp; +} + +} // namespace math +} // namespace stan +#endif From 95a938822bfa2add75ad670b89baeba48d82c4b6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 15 Jun 2025 12:12:33 +0800 Subject: [PATCH 2/7] Update doc, operate on LCDF scale --- .../prob/gaussian_copula_cholesky_lpdf.hpp | 129 +++++++++++------- 1 file changed, 79 insertions(+), 50 deletions(-) diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp index 78a1e7c9a37..a515479faab 100644 --- a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp @@ -3,7 +3,6 @@ #include #include -#include #include #include #include @@ -11,84 +10,106 @@ #include #include #include +#include namespace stan { namespace math { namespace internal { template -auto apply_cdfs_elem_impl(Ty_elem&& y_elem, Ttuple&& cdf_tuple, std::index_sequence){ - auto&& cdf_func = std::get<0>(cdf_tuple); - return cdf_func( +auto apply_lcdfs_elem_impl( + Ty_elem&& y_elem, Ttuple&& lcdf_tuple, std::index_sequence) { + auto&& lcdf_func = std::get<0>(lcdf_tuple); + return lcdf_func( std::forward(y_elem), - std::get(std::forward(cdf_tuple))... + std::get(std::forward(lcdf_tuple))... ); } template -auto apply_cdfs_impl(Ty&& y, Ttuple&& cdf_tuple, std::index_sequence){ - plain_type_t res(y.size()); - - // Use initializer-list expansion to assign the result of each CDF - // to the respective element in the result vector, as we need a constant - // expression for indexing the tuple of CDF functors - static_cast(std::initializer_list{( - res[I] = apply_cdfs_elem_impl( +std::vector> apply_lcdfs_impl( + Ty&& y, Ttuple&& lcdf_tuple, std::index_sequence) { + // Use initializer-list expansion to return a vector with the result of + // applying each LCDF functor to the respective input. + // We cannot use apply_scalar_unary here as we need a constant + // expression for indexing the tuple of LCDF functors + return { + apply_lcdfs_elem_impl( std::forward(y[I]), - std::get(cdf_tuple), - std::make_index_sequence< - // Using size - 1 as the first element is the functor to apply - std::tuple_size< - std::remove_reference_t< - decltype(std::get(cdf_tuple))>>{} - 1>{}), - 0)...}); - - return res; + std::get(lcdf_tuple), + std::make_index_sequence< + // Using size - 1 as the first element is the functor to apply + std::tuple_size< + std::remove_reference_t< + decltype(std::get(lcdf_tuple))>>{} - 1>{} + )... + }; } template -auto apply_cdfs(Ty&& y, Ttuple&& cdf_tuple){ - return apply_cdfs_impl( +auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple){ + return apply_lcdfs_impl( std::forward(y), - std::forward(cdf_tuple), - std::make_index_sequence< - std::tuple_size>{}>{} + std::forward(lcdf_tuple), + std::make_index_sequence< + std::tuple_size>{}>{} ); } } -/* - vectors ~ gaussian_copula_cholesky(cdf_funs_tuple, chol) - (cdf_fun, ...) -*/ -template > +/** \ingroup multivar_dists + * The log of the Gaussian copula density for the given y and + * a Cholesky factor L of the variance matrix. + * + * As the Gaussian copula requires inputs to be in the unit interval, + * users are required to provide a tuple (of the same size as y) where each + * element is a tuple containing a functor for calculating the LCDF and any + * additional parameters required by the LCDF functor. + * + * This version of the function is vectorized on y and the tuple of LCDF + * functor tuples. + * + * @param y A scalar vector + * @param lcdf_fun_tuple A tuple of tuples, where each inner tuple follows the + * structure {functor, param1, param2, ...}, where the functor has signature + * (y[i], param1, param2, ...) and returns the log CDF for the given y[i]. + * @param chol The Cholesky decomposition of a variance matrix + * of the multivariate normal distribution + * @return The log of the gaussian copula density. + * @throw std::domain_error if LL' is not square, not symmetric, + * or not semi-positive definite. + * @tparam T_y Type of scalar. + * @tparam T_lcdf_fun_tuple Type of tuple of LCDF functor tuples. + * @tparam T_chol Type of cholesky factor. + */ +template > T_return gaussian_copula_cholesky_lpdf( - const T_y& y, const T_cdf_fun_tuple& cdf_fun_tuple, const T_chol chol) { + const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { static constexpr const char* function = "gaussian_copula_cholesky_lpdf"; using T_y_ref = ref_type_t; - using T_chochol_ref = ref_type_t; - using T_cdf_ref = ref_type_t; + using T_chol_ref = ref_type_t; + using T_lcdf_ref = ref_type_t; - check_consistent_sizes_mvt(function, "y", y, "cdf_fun_tuple", cdf_fun_tuple); + check_consistent_sizes_mvt(function, "y", y, "lcdf_fun_tuple", lcdf_fun_tuple); const size_t size_mvt_y = size_mvt(y); - const size_t size_mvt_cdf_tuple = size_mvt(cdf_fun_tuple); + const size_t size_mvt_lcdf_tuple = size_mvt(lcdf_fun_tuple); if (size_mvt_y == 0) { return 0; } T_y_ref y_ref = y; - T_chochol_ref chol_ref = chol; - T_cdf_ref cdf_tuple_ref = cdf_fun_tuple; + T_chol_ref chol_ref = chol; + T_lcdf_ref lcdf_tuple_ref = lcdf_fun_tuple; vector_seq_view y_vec(y_ref); - vector_seq_view cdf_tuple_vec(cdf_tuple_ref); - const size_t size_vec = max_size_mvt(y, cdf_fun_tuple); + vector_seq_view lcdf_tuple_vec(lcdf_tuple_ref); + const size_t size_vec = max_size_mvt(y, lcdf_fun_tuple); const int size_y = math::size(y_vec[0]); - const int size_cdf_tuple = math::size(cdf_tuple_vec[0]); + const int size_lcdf_tuple = math::size(lcdf_tuple_vec[0]); // check size consistency of all random variables y for (size_t i = 1; i < size_mvt_y; i++) { @@ -101,18 +122,18 @@ T_return gaussian_copula_cholesky_lpdf( size_y); } // check size consistency of all CDF tuples - for (size_t i = 1; i < size_mvt_cdf_tuple; i++) { + for (size_t i = 1; i < size_mvt_lcdf_tuple; i++) { check_size_match(function, "Size of one of the vectors of " "the CDF functor tuple", - math::size(cdf_tuple_vec[i]), + math::size(lcdf_tuple_vec[i]), "Size of the first vector of the " "location variable", - size_cdf_tuple); + size_lcdf_tuple); } check_size_match(function, "Size of random variable", size_y, - "size of CDF functor tuple", size_cdf_tuple); + "size of CDF functor tuple", size_lcdf_tuple); check_size_match(function, "Size of random variable", size_y, "rows of covariance parameter", chol.rows()); check_size_match(function, "Size of random variable", size_y, @@ -127,14 +148,22 @@ T_return gaussian_copula_cholesky_lpdf( promote_scalar_t> q(size_vec); T_return lp(0); for (size_t i = 0; i < size_vec; i++) { - q[i] = to_vector(inv_Phi(internal::apply_cdfs(y_vec[i], cdf_tuple_vec[i]))); - lp -= std_normal_lpdf(q[i]); + const auto& u = internal::apply_lcdfs(y_vec[i], lcdf_tuple_vec[i]); + check_bounded(function, "CDF-transformed inputs", u, NEGATIVE_INFTY, 0); + q[i] = to_vector(std_normal_log_qf(u)); + lp -= std_normal_lpdf(q[i]); } const std::vector zero_vec(size_vec, rep_vector(0, size_y)); - lp += multi_normal_cholesky_lpdf(q, zero_vec, chol_ref); + lp += multi_normal_cholesky_lpdf(q, zero_vec, chol_ref); return lp; } +template > +T_return gaussian_copula_cholesky_lpdf( + const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { + return gaussian_copula_cholesky_lpdf(y, lcdf_fun_tuple, chol); +} } // namespace math } // namespace stan #endif From 6307acbcc0e75606c663013b61d109bbf13848ff Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 15 Jun 2025 12:51:26 +0800 Subject: [PATCH 3/7] Prim test --- .../prob/gaussian_copula_cholesky_lpdf.hpp | 6 +- .../gaussian_copula_cholesky_lpdf_test.cpp | 123 ++++++++++++++++++ 2 files changed, 126 insertions(+), 3 deletions(-) create mode 100644 test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp index a515479faab..936e8aebc5e 100644 --- a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp @@ -125,7 +125,7 @@ T_return gaussian_copula_cholesky_lpdf( for (size_t i = 1; i < size_mvt_lcdf_tuple; i++) { check_size_match(function, "Size of one of the vectors of " - "the CDF functor tuple", + "the LCDF functor tuple", math::size(lcdf_tuple_vec[i]), "Size of the first vector of the " "location variable", @@ -133,7 +133,7 @@ T_return gaussian_copula_cholesky_lpdf( } check_size_match(function, "Size of random variable", size_y, - "size of CDF functor tuple", size_lcdf_tuple); + "size of LCDF functor tuple", size_lcdf_tuple); check_size_match(function, "Size of random variable", size_y, "rows of covariance parameter", chol.rows()); check_size_match(function, "Size of random variable", size_y, @@ -149,7 +149,7 @@ T_return gaussian_copula_cholesky_lpdf( T_return lp(0); for (size_t i = 0; i < size_vec; i++) { const auto& u = internal::apply_lcdfs(y_vec[i], lcdf_tuple_vec[i]); - check_bounded(function, "CDF-transformed inputs", u, NEGATIVE_INFTY, 0); + check_bounded(function, "LCDF-transformed inputs", u, NEGATIVE_INFTY, 0); q[i] = to_vector(std_normal_log_qf(u)); lp -= std_normal_lpdf(q[i]); } diff --git a/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp new file mode 100644 index 00000000000..e884cab58e5 --- /dev/null +++ b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp @@ -0,0 +1,123 @@ +#include +#include + +TEST(ProbDistributionsGaussianCopulaCholesky, EqualsMultiNormalCholesky) { + using stan::math::gaussian_copula_cholesky_lpdf; + using stan::math::multi_normal_cholesky_lpdf; + using stan::math::diag_pre_multiply; + + Eigen::VectorXd y(2); + y << 1, 3; + + Eigen::VectorXd mu(2); + mu << 0.1, 0; + + Eigen::VectorXd sd(2); + sd << 2, 1; + + Eigen::MatrixXd chol(2, 2); + chol << 1, 0, 0.5, 0.8660254037844386; + + auto norm_lcdf = [](auto&& y, auto&& mu, auto&& sigma) { + return stan::math::normal_lcdf(y, mu, sigma); + }; + + auto std_norm_lcdf = [](auto&& y) { + return stan::math::std_normal_lcdf(y); + }; + + // y[0] ~ normal(0.1, 2) + // y[1] ~ normal(0, 1) + auto lcdf_functors = std::make_tuple( + std::make_tuple(norm_lcdf, mu[0], sd[0]), + std::make_tuple(std_norm_lcdf) + ); + + double log_prob = + stan::math::normal_lpdf(y[0], 0.1, 2) + + stan::math::std_normal_lpdf(y[1]) + + gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); + + + double expected_log_prob = multi_normal_cholesky_lpdf( + y, mu, diag_pre_multiply(sd, chol)); + + EXPECT_FLOAT_EQ(log_prob, expected_log_prob); +} + +TEST(ProbDistributionsGaussianCopulaCholesky, NonNormalMarginals) { + Eigen::VectorXd y(2); + y << 10, 6; + + Eigen::MatrixXd chol(2, 2); + chol << 1, 0, 0.2, 0.9797958971; + + auto gamma_lcdf = [](auto&& y, auto&& shape, auto&& scale) { + return stan::math::gamma_lcdf(y, shape, scale); + }; + + auto exp_lcdf = [](auto&& y, auto&& scale) { + return stan::math::exponential_lcdf(y, scale); + }; + + // y[0] ~ gamma(2, 1) + // y[1] ~ exponential(2) + auto lcdf_functors = std::make_tuple( + std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple(exp_lcdf, 2.0) + ); + + double log_prob = + stan::math::gamma_lpdf(y[0], 2.0, 1.0) + + stan::math::exponential_lpdf(y[1], 2.0) + + stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); + + EXPECT_FLOAT_EQ(log_prob, -16.61005941); +} + +TEST(ProbDistributionsGaussianCopulaCholesky, Errors) { + Eigen::VectorXd y(2); + y << 10, 6; + + Eigen::VectorXd small_y(1); + small_y << 10; + + Eigen::MatrixXd chol(2, 2); + chol << 1, 0, 0.2, 0.9797958971; + + auto gamma_lcdf = [](auto&& y, auto&& shape, auto&& scale) { + return stan::math::gamma_lcdf(y, shape, scale); + }; + auto exp_lcdf = [](auto&& y, auto&& scale) { + return stan::math::exponential_lcdf(y, scale); + }; + + auto lcdf_functors = std::make_tuple( + std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple(exp_lcdf, 2.0) + ); + + auto small_lcdf_functors = std::make_tuple( + std::make_tuple(gamma_lcdf, 2.0, 1.0) + ); + + auto invalid_lcdf_functors = std::make_tuple( + std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple([](auto&& y){ return y * 10; }) + ); + + EXPECT_THROW( + stan::math::gaussian_copula_cholesky_lpdf(y, small_lcdf_functors, chol), + std::invalid_argument + ); + + EXPECT_THROW( + stan::math::gaussian_copula_cholesky_lpdf(small_y, lcdf_functors, chol), + std::invalid_argument + ); + + EXPECT_THROW( + stan::math::gaussian_copula_cholesky_lpdf(y, invalid_lcdf_functors, chol), + std::domain_error + ); +} From 71dbf0029e4fa2bf3fe632853864f66a67954d40 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 15 Jun 2025 20:21:47 +0800 Subject: [PATCH 4/7] Initial mix, fwd failing --- .../prob/gaussian_copula_cholesky_lpdf.hpp | 71 +++++++++---------- .../prob/gaussian_copula_cholesky_test.cpp | 35 +++++++++ 2 files changed, 67 insertions(+), 39 deletions(-) create mode 100644 test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp index 936e8aebc5e..9562f006253 100644 --- a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp @@ -26,35 +26,17 @@ auto apply_lcdfs_elem_impl( ); } - -template -std::vector> apply_lcdfs_impl( - Ty&& y, Ttuple&& lcdf_tuple, std::index_sequence) { - // Use initializer-list expansion to return a vector with the result of - // applying each LCDF functor to the respective input. - // We cannot use apply_scalar_unary here as we need a constant - // expression for indexing the tuple of LCDF functors - return { - apply_lcdfs_elem_impl( - std::forward(y[I]), - std::get(lcdf_tuple), - std::make_index_sequence< - // Using size - 1 as the first element is the functor to apply - std::tuple_size< - std::remove_reference_t< - decltype(std::get(lcdf_tuple))>>{} - 1>{} - )... - }; -} - template -auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple){ - return apply_lcdfs_impl( - std::forward(y), - std::forward(lcdf_tuple), - std::make_index_sequence< - std::tuple_size>{}>{} - ); +auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple) { + return index_apply>{}>( + [&y, &lcdf_tuple](auto... Is) { + return std::make_tuple( + apply_lcdfs_elem_impl( + y[Is], std::get(lcdf_tuple), std::make_index_sequence< + std::tuple_size(lcdf_tuple))>>{} - 1>{} + )... + ); + }); } } @@ -84,9 +66,8 @@ auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple){ * @tparam T_lcdf_fun_tuple Type of tuple of LCDF functor tuples. * @tparam T_chol Type of cholesky factor. */ -template > -T_return gaussian_copula_cholesky_lpdf( +template +auto gaussian_copula_cholesky_lpdf( const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { static constexpr const char* function = "gaussian_copula_cholesky_lpdf"; @@ -97,15 +78,17 @@ T_return gaussian_copula_cholesky_lpdf( check_consistent_sizes_mvt(function, "y", y, "lcdf_fun_tuple", lcdf_fun_tuple); const size_t size_mvt_y = size_mvt(y); const size_t size_mvt_lcdf_tuple = size_mvt(lcdf_fun_tuple); - if (size_mvt_y == 0) { - return 0; - } T_y_ref y_ref = y; T_chol_ref chol_ref = chol; T_lcdf_ref lcdf_tuple_ref = lcdf_fun_tuple; vector_seq_view y_vec(y_ref); vector_seq_view lcdf_tuple_vec(lcdf_tuple_ref); + using T_return = return_type_t; + + if (size_mvt_y == 0) { + return T_return(0); + } const size_t size_vec = max_size_mvt(y, lcdf_fun_tuple); const int size_y = math::size(y_vec[0]); @@ -148,9 +131,20 @@ T_return gaussian_copula_cholesky_lpdf( promote_scalar_t> q(size_vec); T_return lp(0); for (size_t i = 0; i < size_vec; i++) { - const auto& u = internal::apply_lcdfs(y_vec[i], lcdf_tuple_vec[i]); + const auto& y_i = y_vec[i]; + const auto& func_i = lcdf_tuple_vec[i]; + + const auto& res = internal::apply_lcdfs(y_i, func_i); + const auto& u = index_apply>{}>( + [&res, size_y](auto... Is) { + Eigen::Matrix u_inner(size_y); + static_cast(std::initializer_list{ + (u_inner[Is] = std::get(res), 0)... + }); + return u_inner; + }); check_bounded(function, "LCDF-transformed inputs", u, NEGATIVE_INFTY, 0); - q[i] = to_vector(std_normal_log_qf(u)); + q[i] = std_normal_log_qf(u); lp -= std_normal_lpdf(q[i]); } const std::vector zero_vec(size_vec, rep_vector(0, size_y)); @@ -158,9 +152,8 @@ T_return gaussian_copula_cholesky_lpdf( return lp; } -template > -T_return gaussian_copula_cholesky_lpdf( +template +auto gaussian_copula_cholesky_lpdf( const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { return gaussian_copula_cholesky_lpdf(y, lcdf_fun_tuple, chol); } diff --git a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp new file mode 100644 index 00000000000..0abcf4b931c --- /dev/null +++ b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp @@ -0,0 +1,35 @@ +#include + +TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { + // Bind functors for compatibility with AD framework + auto f = [](const auto func1, const auto func2) { + return [=](const auto& y, const auto& args, const auto& sigma) { + auto lcdf_functors = std::make_tuple( + std::make_tuple(func1, args[0]), + std::make_tuple(func2, args[1]) + ); + auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose()); + auto L = stan::math::cholesky_decompose(sigma_sym); + return stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, L); + }; + }; + + // y[0] ~ gamma(2, 1) + // y[1] ~ std_normal() + Eigen::VectorXd y1(2); + y1 << 1.0, 0.1; + Eigen::VectorXd args(2); + args << 0, 2.0; + + auto func1 = [](const auto& y, const auto& mu) { + return stan::math::normal_lcdf(y, mu, 1.0); + }; + auto func2 = [](const auto& y) { + return stan::math::std_normal_lcdf(y); + }; + + Eigen::MatrixXd Sigma22(2, 2); + Sigma22 << 2.0, 0.5, 0.5, 1.1; + + stan::test::expect_ad(f(func2, func2), y1, args, Sigma22); +} From d362e3c331fb12d7722f3a775139f76474ce5aa0 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 15 Jun 2025 08:43:32 -0400 Subject: [PATCH 5/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/size.hpp | 1 - stan/math/prim/fun/vector_seq_view.hpp | 3 +- .../prob/gaussian_copula_cholesky_lpdf.hpp | 59 +++++++-------- .../prob/gaussian_copula_cholesky_test.cpp | 10 +-- .../gaussian_copula_cholesky_lpdf_test.cpp | 71 ++++++++----------- 5 files changed, 62 insertions(+), 82 deletions(-) diff --git a/stan/math/prim/fun/size.hpp b/stan/math/prim/fun/size.hpp index 8a5418ecc24..34bb9f42358 100644 --- a/stan/math/prim/fun/size.hpp +++ b/stan/math/prim/fun/size.hpp @@ -34,7 +34,6 @@ inline int64_t size(const T& m) { return m.size(); } - template * = nullptr> inline int64_t size(const T& /*t*/) { return std::tuple_size>{}; diff --git a/stan/math/prim/fun/vector_seq_view.hpp b/stan/math/prim/fun/vector_seq_view.hpp index b11d8e291f3..a28b916c7e3 100644 --- a/stan/math/prim/fun/vector_seq_view.hpp +++ b/stan/math/prim/fun/vector_seq_view.hpp @@ -13,8 +13,7 @@ using is_matrix_or_std_vector_or_tuple template using is_scalar_container = math::disjunction< - is_matrix, - math::is_tuple, + is_matrix, math::is_tuple, math::conjunction, is_stan_scalar>>>; } // namespace internal diff --git a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp index 9562f006253..2a068ced450 100644 --- a/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp +++ b/stan/math/prim/prob/gaussian_copula_cholesky_lpdf.hpp @@ -17,29 +17,25 @@ namespace math { namespace internal { template -auto apply_lcdfs_elem_impl( - Ty_elem&& y_elem, Ttuple&& lcdf_tuple, std::index_sequence) { +auto apply_lcdfs_elem_impl(Ty_elem&& y_elem, Ttuple&& lcdf_tuple, + std::index_sequence) { auto&& lcdf_func = std::get<0>(lcdf_tuple); - return lcdf_func( - std::forward(y_elem), - std::get(std::forward(lcdf_tuple))... - ); + return lcdf_func(std::forward(y_elem), + std::get(std::forward(lcdf_tuple))...); } template auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple) { return index_apply>{}>( - [&y, &lcdf_tuple](auto... Is) { - return std::make_tuple( - apply_lcdfs_elem_impl( - y[Is], std::get(lcdf_tuple), std::make_index_sequence< - std::tuple_size(lcdf_tuple))>>{} - 1>{} - )... - ); - }); -} + [&y, &lcdf_tuple](auto... Is) { + return std::make_tuple(apply_lcdfs_elem_impl( + y[Is], std::get(lcdf_tuple), + std::make_index_sequence(lcdf_tuple))>>{} + - 1>{})...); + }); } - +} // namespace internal /** \ingroup multivar_dists * The log of the Gaussian copula density for the given y and @@ -67,15 +63,17 @@ auto apply_lcdfs(Ty&& y, Ttuple&& lcdf_tuple) { * @tparam T_chol Type of cholesky factor. */ template -auto gaussian_copula_cholesky_lpdf( - const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { +auto gaussian_copula_cholesky_lpdf(const T_y& y, + const T_lcdf_fun_tuple& lcdf_fun_tuple, + const T_chol chol) { static constexpr const char* function = "gaussian_copula_cholesky_lpdf"; using T_y_ref = ref_type_t; using T_chol_ref = ref_type_t; using T_lcdf_ref = ref_type_t; - check_consistent_sizes_mvt(function, "y", y, "lcdf_fun_tuple", lcdf_fun_tuple); + check_consistent_sizes_mvt(function, "y", y, "lcdf_fun_tuple", + lcdf_fun_tuple); const size_t size_mvt_y = size_mvt(y); const size_t size_mvt_lcdf_tuple = size_mvt(lcdf_fun_tuple); T_y_ref y_ref = y; @@ -84,7 +82,9 @@ auto gaussian_copula_cholesky_lpdf( vector_seq_view y_vec(y_ref); vector_seq_view lcdf_tuple_vec(lcdf_tuple_ref); - using T_return = return_type_t; + using T_return = return_type_t< + T_y, decltype(internal::apply_lcdfs(y_vec[0], lcdf_tuple_vec[0])), + T_chol>; if (size_mvt_y == 0) { return T_return(0); @@ -135,14 +135,14 @@ auto gaussian_copula_cholesky_lpdf( const auto& func_i = lcdf_tuple_vec[i]; const auto& res = internal::apply_lcdfs(y_i, func_i); - const auto& u = index_apply>{}>( - [&res, size_y](auto... Is) { - Eigen::Matrix u_inner(size_y); - static_cast(std::initializer_list{ - (u_inner[Is] = std::get(res), 0)... + const auto& u = index_apply< + std::tuple_size>{}>( + [&res, size_y](auto... Is) { + Eigen::Matrix u_inner(size_y); + static_cast(std::initializer_list{ + (u_inner[Is] = std::get(res), 0)...}); + return u_inner; }); - return u_inner; - }); check_bounded(function, "LCDF-transformed inputs", u, NEGATIVE_INFTY, 0); q[i] = std_normal_log_qf(u); lp -= std_normal_lpdf(q[i]); @@ -153,8 +153,9 @@ auto gaussian_copula_cholesky_lpdf( } template -auto gaussian_copula_cholesky_lpdf( - const T_y& y, const T_lcdf_fun_tuple& lcdf_fun_tuple, const T_chol chol) { +auto gaussian_copula_cholesky_lpdf(const T_y& y, + const T_lcdf_fun_tuple& lcdf_fun_tuple, + const T_chol chol) { return gaussian_copula_cholesky_lpdf(y, lcdf_fun_tuple, chol); } } // namespace math diff --git a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp index 0abcf4b931c..8ff40ccd063 100644 --- a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp +++ b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp @@ -4,10 +4,8 @@ TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { // Bind functors for compatibility with AD framework auto f = [](const auto func1, const auto func2) { return [=](const auto& y, const auto& args, const auto& sigma) { - auto lcdf_functors = std::make_tuple( - std::make_tuple(func1, args[0]), - std::make_tuple(func2, args[1]) - ); + auto lcdf_functors = std::make_tuple(std::make_tuple(func1, args[0]), + std::make_tuple(func2, args[1])); auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose()); auto L = stan::math::cholesky_decompose(sigma_sym); return stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, L); @@ -24,9 +22,7 @@ TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { auto func1 = [](const auto& y, const auto& mu) { return stan::math::normal_lcdf(y, mu, 1.0); }; - auto func2 = [](const auto& y) { - return stan::math::std_normal_lcdf(y); - }; + auto func2 = [](const auto& y) { return stan::math::std_normal_lcdf(y); }; Eigen::MatrixXd Sigma22(2, 2); Sigma22 << 2.0, 0.5, 0.5, 1.1; diff --git a/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp index e884cab58e5..b5a71abc6f2 100644 --- a/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp +++ b/test/unit/math/prim/prob/gaussian_copula_cholesky_lpdf_test.cpp @@ -2,9 +2,9 @@ #include TEST(ProbDistributionsGaussianCopulaCholesky, EqualsMultiNormalCholesky) { + using stan::math::diag_pre_multiply; using stan::math::gaussian_copula_cholesky_lpdf; using stan::math::multi_normal_cholesky_lpdf; - using stan::math::diag_pre_multiply; Eigen::VectorXd y(2); y << 1, 3; @@ -22,25 +22,19 @@ TEST(ProbDistributionsGaussianCopulaCholesky, EqualsMultiNormalCholesky) { return stan::math::normal_lcdf(y, mu, sigma); }; - auto std_norm_lcdf = [](auto&& y) { - return stan::math::std_normal_lcdf(y); - }; + auto std_norm_lcdf = [](auto&& y) { return stan::math::std_normal_lcdf(y); }; // y[0] ~ normal(0.1, 2) // y[1] ~ normal(0, 1) - auto lcdf_functors = std::make_tuple( - std::make_tuple(norm_lcdf, mu[0], sd[0]), - std::make_tuple(std_norm_lcdf) - ); - - double log_prob = - stan::math::normal_lpdf(y[0], 0.1, 2) + - stan::math::std_normal_lpdf(y[1]) + - gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); + auto lcdf_functors = std::make_tuple(std::make_tuple(norm_lcdf, mu[0], sd[0]), + std::make_tuple(std_norm_lcdf)); + double log_prob = stan::math::normal_lpdf(y[0], 0.1, 2) + + stan::math::std_normal_lpdf(y[1]) + + gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); - double expected_log_prob = multi_normal_cholesky_lpdf( - y, mu, diag_pre_multiply(sd, chol)); + double expected_log_prob + = multi_normal_cholesky_lpdf(y, mu, diag_pre_multiply(sd, chol)); EXPECT_FLOAT_EQ(log_prob, expected_log_prob); } @@ -62,15 +56,13 @@ TEST(ProbDistributionsGaussianCopulaCholesky, NonNormalMarginals) { // y[0] ~ gamma(2, 1) // y[1] ~ exponential(2) - auto lcdf_functors = std::make_tuple( - std::make_tuple(gamma_lcdf, 2.0, 1.0), - std::make_tuple(exp_lcdf, 2.0) - ); + auto lcdf_functors = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple(exp_lcdf, 2.0)); - double log_prob = - stan::math::gamma_lpdf(y[0], 2.0, 1.0) + - stan::math::exponential_lpdf(y[1], 2.0) + - stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); + double log_prob + = stan::math::gamma_lpdf(y[0], 2.0, 1.0) + + stan::math::exponential_lpdf(y[1], 2.0) + + stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, chol); EXPECT_FLOAT_EQ(log_prob, -16.61005941); } @@ -92,32 +84,25 @@ TEST(ProbDistributionsGaussianCopulaCholesky, Errors) { return stan::math::exponential_lcdf(y, scale); }; - auto lcdf_functors = std::make_tuple( - std::make_tuple(gamma_lcdf, 2.0, 1.0), - std::make_tuple(exp_lcdf, 2.0) - ); + auto lcdf_functors = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple(exp_lcdf, 2.0)); - auto small_lcdf_functors = std::make_tuple( - std::make_tuple(gamma_lcdf, 2.0, 1.0) - ); + auto small_lcdf_functors + = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0)); - auto invalid_lcdf_functors = std::make_tuple( - std::make_tuple(gamma_lcdf, 2.0, 1.0), - std::make_tuple([](auto&& y){ return y * 10; }) - ); + auto invalid_lcdf_functors + = std::make_tuple(std::make_tuple(gamma_lcdf, 2.0, 1.0), + std::make_tuple([](auto&& y) { return y * 10; })); EXPECT_THROW( - stan::math::gaussian_copula_cholesky_lpdf(y, small_lcdf_functors, chol), - std::invalid_argument - ); + stan::math::gaussian_copula_cholesky_lpdf(y, small_lcdf_functors, chol), + std::invalid_argument); EXPECT_THROW( - stan::math::gaussian_copula_cholesky_lpdf(small_y, lcdf_functors, chol), - std::invalid_argument - ); + stan::math::gaussian_copula_cholesky_lpdf(small_y, lcdf_functors, chol), + std::invalid_argument); EXPECT_THROW( - stan::math::gaussian_copula_cholesky_lpdf(y, invalid_lcdf_functors, chol), - std::domain_error - ); + stan::math::gaussian_copula_cholesky_lpdf(y, invalid_lcdf_functors, chol), + std::domain_error); } From 8b6806676b2a2c0692aed703b3e885d9ee66ad49 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 16 Jun 2025 23:02:05 +0800 Subject: [PATCH 6/7] Update test construction --- .../math/mix/prob/gaussian_copula_cholesky_test.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp index 8ff40ccd063..c373ac7cd7e 100644 --- a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp +++ b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp @@ -4,8 +4,8 @@ TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { // Bind functors for compatibility with AD framework auto f = [](const auto func1, const auto func2) { return [=](const auto& y, const auto& args, const auto& sigma) { - auto lcdf_functors = std::make_tuple(std::make_tuple(func1, args[0]), - std::make_tuple(func2, args[1])); + auto lcdf_functors = std::make_tuple(std::make_tuple(func1, args[0], args[1]), + std::make_tuple(func2)); auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose()); auto L = stan::math::cholesky_decompose(sigma_sym); return stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, L); @@ -17,15 +17,15 @@ TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { Eigen::VectorXd y1(2); y1 << 1.0, 0.1; Eigen::VectorXd args(2); - args << 0, 2.0; + args << 2, 1; - auto func1 = [](const auto& y, const auto& mu) { - return stan::math::normal_lcdf(y, mu, 1.0); + auto func1 = [](const auto& y, const auto& shape, const auto& scale) { + return stan::math::gamma_lcdf(y, shape, scale); }; auto func2 = [](const auto& y) { return stan::math::std_normal_lcdf(y); }; Eigen::MatrixXd Sigma22(2, 2); Sigma22 << 2.0, 0.5, 0.5, 1.1; - stan::test::expect_ad(f(func2, func2), y1, args, Sigma22); + stan::test::expect_ad(f(func1, func2), y1, args, Sigma22); } From c85abc54822e17a32fc309d78b82bb7aedc331b1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 16 Jun 2025 11:03:07 -0400 Subject: [PATCH 7/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp index c373ac7cd7e..60eb8f8bef8 100644 --- a/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp +++ b/test/unit/math/mix/prob/gaussian_copula_cholesky_test.cpp @@ -4,8 +4,8 @@ TEST_F(AgradRev, ProbDistributionsGaussCopulaCholesky) { // Bind functors for compatibility with AD framework auto f = [](const auto func1, const auto func2) { return [=](const auto& y, const auto& args, const auto& sigma) { - auto lcdf_functors = std::make_tuple(std::make_tuple(func1, args[0], args[1]), - std::make_tuple(func2)); + auto lcdf_functors = std::make_tuple( + std::make_tuple(func1, args[0], args[1]), std::make_tuple(func2)); auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose()); auto L = stan::math::cholesky_decompose(sigma_sym); return stan::math::gaussian_copula_cholesky_lpdf(y, lcdf_functors, L);