From 420012011f698f375f7a86794ccbaf134865efdf Mon Sep 17 00:00:00 2001 From: jachym Date: Mon, 3 Mar 2025 21:42:40 +0100 Subject: [PATCH 1/9] Generalized normal lpdf --- .../prim/prob/double_exponential_lpdf.hpp | 20 ++- .../prim/prob/generalized_normal_lpdf.hpp | 138 ++++++++++++++++++ stan/math/prim/prob/normal_lpdf.hpp | 14 +- .../double_exponential_test.hpp | 1 - test/prob/normal/normal_test.hpp | 1 - 5 files changed, 153 insertions(+), 21 deletions(-) create mode 100644 stan/math/prim/prob/generalized_normal_lpdf.hpp diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp index 059eebdeb4e..1592e786d38 100644 --- a/stan/math/prim/prob/double_exponential_lpdf.hpp +++ b/stan/math/prim/prob/double_exponential_lpdf.hpp @@ -3,11 +3,9 @@ #include #include -#include -#include #include #include -#include +#include #include #include #include @@ -73,11 +71,11 @@ return_type_t double_exponential_lpdf( const auto& inv_sigma = to_ref(inv(sigma_val)); const auto& y_m_mu = to_ref_if::value>(y_val - mu_val); - const auto& abs_diff_y_mu = fabs(y_m_mu); + const auto& abs_diff_y_mu = abs(y_m_mu); const auto& scaled_diff - = to_ref_if::value>(abs_diff_y_mu * inv_sigma); + = to_ref_if::value>(abs_diff_y_mu * inv_sigma); - size_t N = max_size(y, mu, sigma); + const size_t N = max_size(y, mu, sigma); if (include_summand::value) { logp -= N * LOG_TWO; } @@ -89,16 +87,16 @@ return_type_t double_exponential_lpdf( if (!is_constant_all::value) { const auto& diff_sign = sign(y_m_mu); const auto& rep_deriv - = to_ref_if<(!is_constant_all::value - && !is_constant_all::value)>(diff_sign * inv_sigma); - if (!is_constant_all::value) { + = to_ref_if<(!is_constant::value + && !is_constant::value)>(diff_sign * inv_sigma); + if (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; } - if (!is_constant_all::value) { + if (!is_constant::value) { partials<1>(ops_partials) = rep_deriv; } } - if (!is_constant_all::value) { + if (!is_constant::value) { partials<2>(ops_partials) = inv_sigma * (scaled_diff - 1); } diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp new file mode 100644 index 00000000000..a17e8452dc2 --- /dev/null +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -0,0 +1,138 @@ +#ifndef STAN_MATH_PRIM_PROB_GENERALIZED_NORMAL_LPDF_HPP +#define STAN_MATH_PRIM_PROB_GENERALIZED_NORMAL_LPDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The log of the generalized normal density for the specified scalar(s) given + * the specified location, scale and shape parameters. y, mu, alpha, or beta can + * each be either a scalar or a vector. Any vector inputs must be the same length. + * + *

The result log probability is defined to be the sum of the + * log probabilities for each observation/mean/scale/shape tuple. + * + * @tparam T_y type of scalar + * @tparam T_loc type of location parameter + * @tparam T_scale type of scale parameter + * @tparam T_shape type of shape parameter + * @param y (Sequence of) scalar(s) + * @param mu (Sequence of) location parameter(s) + * @param alpha (Sequence of) scale parameter(s) + * @param beta (Sequence of) shape parameter(s) + * @return The log of the product of the densities + * @throw std::domain_error if alpha or beta is not positive + */ +template * = nullptr> +inline return_type_t generalized_normal_lpdf( + T_y&& y, T_loc&& mu, T_scale&& alpha, T_shape&& beta) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_if_not_constant_t; + using T_mu_ref = ref_type_if_not_constant_t; + using T_alpha_ref = ref_type_if_not_constant_t; + using T_beta_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "generalized_normal_lpdf"; + check_consistent_sizes(function, "Random variable", y, "Location parameter", mu, + "Scale parameter", alpha, "Shape parameter", beta); + + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_alpha_ref alpha_ref = std::forward(alpha); + T_beta_ref beta_ref = std::forward(beta); + + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref)); + decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref)); + decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref)); + + check_not_nan(function, "Random variable", y_val); + check_finite(function, "Location parameter", mu_val); + check_positive(function, "Scale parameter", alpha_val); + check_positive(function, "Shape parameter", beta_val); + + if (size_zero(y, mu, alpha, beta)) { + return 0; + } + if (!include_summand::value) { + return 0; + } + + auto ops_partials = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); + const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); + const auto& diff = to_ref_if::value>(y_val - mu_val); + const auto& scaled_abs_diff = to_ref_if::value>(abs(diff) / alpha_val); + const auto& scaled_abs_diff_pow = to_ref(pow(scaled_abs_diff, beta_val)); + const size_t N = max_size(y, mu, alpha, beta); + + T_partials_return logp = -sum(scaled_abs_diff_pow); + + if (include_summand::value) { + logp -= LOG_TWO * N; + } + if (include_summand::value) { + logp -= sum(log(alpha_val)) * (N / math::size(alpha)); + } + if (include_summand::value) { + logp -= sum(lgamma(inv_beta1p)) * (N / math::size(beta)); + } + + if (!is_constant_all::value) { + const auto& beta_scaled_abs_diff_pow = to_ref_if<(!is_constant_all::value + + !is_constant_all::value + + !is_constant_all::value + >= 2)>(beta_val * scaled_abs_diff_pow); + + if (!is_constant_all::value) { + const auto& deriv = to_ref_if<(!is_constant::value + && !is_constant::value)>(sign(diff) * beta_scaled_abs_diff_pow / alpha_val); + if (!is_constant::value) { + partials<0>(ops_partials) = -deriv; + } + if (!is_constant::value) { + partials<1>(ops_partials) = std::move(deriv); + } + } + if (!is_constant::value) { + partials<2>(ops_partials) = (beta_scaled_abs_diff_pow - 1) / alpha_val; + } + } + if (!is_constant::value) { + partials<3>(ops_partials) = digamma(inv_beta1p) / square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); + } + + return ops_partials.build(logp); +} + +template +inline return_type_t generalized_normal_lpdf( + T_y&& y, T_loc&& mu, T_scale&& alpha, T_shape&& beta) { + return generalized_normal_lpdf(std::forward(y), + std::forward(mu), + std::forward(alpha), + std::forward(beta)); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index f53c521a74e..83242104591 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -3,8 +3,6 @@ #include #include -#include -#include #include #include #include @@ -87,17 +85,17 @@ inline return_type_t normal_lpdf(T_y&& y, T_loc&& mu, } if (!is_constant_all::value) { - auto scaled_diff = to_ref_if::value - + !is_constant_all::value - + !is_constant_all::value + auto scaled_diff = to_ref_if::value + + !is_constant::value + + !is_constant::value >= 2>(inv_sigma * y_scaled); - if (!is_constant_all::value) { + if (!is_constant::value) { partials<0>(ops_partials) = -scaled_diff; } - if (!is_constant_all::value) { + if (!is_constant::value) { partials<2>(ops_partials) = inv_sigma * y_scaled_sq - inv_sigma; } - if (!is_constant_all::value) { + if (!is_constant::value) { partials<1>(ops_partials) = std::move(scaled_diff); } } diff --git a/test/prob/double_exponential/double_exponential_test.hpp b/test/prob/double_exponential/double_exponential_test.hpp index db5aa235a59..be206376758 100644 --- a/test/prob/double_exponential/double_exponential_test.hpp +++ b/test/prob/double_exponential/double_exponential_test.hpp @@ -2,7 +2,6 @@ #include #include -using stan::math::var; using std::numeric_limits; using std::vector; diff --git a/test/prob/normal/normal_test.hpp b/test/prob/normal/normal_test.hpp index 81e2ca75927..df818f10e02 100644 --- a/test/prob/normal/normal_test.hpp +++ b/test/prob/normal/normal_test.hpp @@ -3,7 +3,6 @@ #include #include -using stan::math::var; using std::numeric_limits; using std::vector; From 84ea5eddf1d7fa3fc0511fe18bec43048b61ff39 Mon Sep 17 00:00:00 2001 From: jachym Date: Tue, 4 Mar 2025 16:32:19 +0100 Subject: [PATCH 2/9] Generalized normal lpdf - fix corner cases --- .../prim/prob/double_exponential_lpdf.hpp | 14 +++---- .../prim/prob/generalized_normal_lpdf.hpp | 42 ++++++++++--------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp index 1592e786d38..2169a27afda 100644 --- a/stan/math/prim/prob/double_exponential_lpdf.hpp +++ b/stan/math/prim/prob/double_exponential_lpdf.hpp @@ -57,7 +57,6 @@ return_type_t double_exponential_lpdf( return 0.0; } - T_partials_return logp(0.0); auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); @@ -71,10 +70,10 @@ return_type_t double_exponential_lpdf( const auto& inv_sigma = to_ref(inv(sigma_val)); const auto& y_m_mu = to_ref_if::value>(y_val - mu_val); - const auto& abs_diff_y_mu = abs(y_m_mu); - const auto& scaled_diff - = to_ref_if::value>(abs_diff_y_mu * inv_sigma); + const auto& scaled_abs_diff + = to_ref_if::value>(abs(y_m_mu) * inv_sigma); + T_partials_return logp = -sum(scaled_abs_diff); const size_t N = max_size(y, mu, sigma); if (include_summand::value) { logp -= N * LOG_TWO; @@ -82,13 +81,12 @@ return_type_t double_exponential_lpdf( if (include_summand::value) { logp -= sum(log(sigma_val)) * N / math::size(sigma); } - logp -= sum(scaled_diff); if (!is_constant_all::value) { - const auto& diff_sign = sign(y_m_mu); const auto& rep_deriv = to_ref_if<(!is_constant::value - && !is_constant::value)>(diff_sign * inv_sigma); + && !is_constant::value)>(sign(y_m_mu) * inv_sigma); + if (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; } @@ -97,7 +95,7 @@ return_type_t double_exponential_lpdf( } } if (!is_constant::value) { - partials<2>(ops_partials) = inv_sigma * (scaled_diff - 1); + partials<2>(ops_partials) = inv_sigma * (scaled_abs_diff - 1); } return ops_partials.build(logp); diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index a17e8452dc2..3102d12ba98 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -8,9 +8,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -78,11 +80,11 @@ inline return_type_t generalized_normal_lpdf( return 0; } - auto ops_partials = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); - const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); + const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); // sus T_shape const auto& diff = to_ref_if::value>(y_val - mu_val); - const auto& scaled_abs_diff = to_ref_if::value>(abs(diff) / alpha_val); - const auto& scaled_abs_diff_pow = to_ref(pow(scaled_abs_diff, beta_val)); + const auto& inv_alpha = to_ref(inv(alpha_val)); + const auto& scaled_abs_diff = to_ref_if::value>(abs(diff) * inv_alpha); + const auto& scaled_abs_diff_pow = to_ref_if::value>(pow(scaled_abs_diff, beta_val)); const size_t N = max_size(y, mu, alpha, beta); T_partials_return logp = -sum(scaled_abs_diff_pow); @@ -97,26 +99,26 @@ inline return_type_t generalized_normal_lpdf( logp -= sum(lgamma(inv_beta1p)) * (N / math::size(beta)); } - if (!is_constant_all::value) { - const auto& beta_scaled_abs_diff_pow = to_ref_if<(!is_constant_all::value - + !is_constant_all::value - + !is_constant_all::value - >= 2)>(beta_val * scaled_abs_diff_pow); + auto ops_partials = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); - if (!is_constant_all::value) { - const auto& deriv = to_ref_if<(!is_constant::value - && !is_constant::value)>(sign(diff) * beta_scaled_abs_diff_pow / alpha_val); - if (!is_constant::value) { - partials<0>(ops_partials) = -deriv; - } - if (!is_constant::value) { - partials<1>(ops_partials) = std::move(deriv); - } + if (!is_constant_all::value) { + // note: The partial derivatives for y, mu are undefined when y == mu && beta < 1. + // The derivative limit as mu -> y goes: + // to 0 from both sides if beta > 1 (defined as 0) + // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as 0, consistent with double_exponential_lpdf) + // to +∞ from right but -∞ from left as y -> mu if beta < 1 (undefined) + const auto& rep_deriv = to_ref_if::value + && !is_constant::value>(sign(diff) * beta_val * pow(scaled_abs_diff, beta_val-1) * inv_alpha); + if (!is_constant::value) { + partials<0>(ops_partials) = -rep_deriv; } - if (!is_constant::value) { - partials<2>(ops_partials) = (beta_scaled_abs_diff_pow - 1) / alpha_val; + if (!is_constant::value) { + partials<1>(ops_partials) = std::move(rep_deriv); } } + if (!is_constant::value) { + partials<2>(ops_partials) = (beta_val * scaled_abs_diff_pow - 1) * inv_alpha; + } if (!is_constant::value) { partials<3>(ops_partials) = digamma(inv_beta1p) / square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); } From 9d9b1cf353c3c3a42f151d086cd2980561d0815c Mon Sep 17 00:00:00 2001 From: jachym Date: Tue, 4 Mar 2025 16:39:16 +0100 Subject: [PATCH 3/9] Generalized normal lpdf - fix corner cases --- stan/math/prim/prob/generalized_normal_lpdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index 3102d12ba98..5286d287c4e 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -120,7 +120,7 @@ inline return_type_t generalized_normal_lpdf( partials<2>(ops_partials) = (beta_val * scaled_abs_diff_pow - 1) * inv_alpha; } if (!is_constant::value) { - partials<3>(ops_partials) = digamma(inv_beta1p) / square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); + partials<3>(ops_partials) = digamma(inv_beta1p) * inv_square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); } return ops_partials.build(logp); From 383b51d97c1919408332e764d5d6ecbe6f63e5d3 Mon Sep 17 00:00:00 2001 From: jachym Date: Tue, 4 Mar 2025 16:50:10 +0100 Subject: [PATCH 4/9] comments --- stan/math/prim/prob/generalized_normal_lpdf.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index 5286d287c4e..7ffbadcc6fb 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -71,7 +71,7 @@ inline return_type_t generalized_normal_lpdf( check_not_nan(function, "Random variable", y_val); check_finite(function, "Location parameter", mu_val); check_positive(function, "Scale parameter", alpha_val); - check_positive(function, "Shape parameter", beta_val); + check_positive(function, "Shape parameter", beta_val); // With β = +∞ this could be defined to be uniform, but we don't support that. if (size_zero(y, mu, alpha, beta)) { return 0; @@ -80,7 +80,7 @@ inline return_type_t generalized_normal_lpdf( return 0; } - const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); // sus T_shape + const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); const auto& diff = to_ref_if::value>(y_val - mu_val); const auto& inv_alpha = to_ref(inv(alpha_val)); const auto& scaled_abs_diff = to_ref_if::value>(abs(diff) * inv_alpha); From 180b1538870c99e45d7dc4b4cf627440819a0071 Mon Sep 17 00:00:00 2001 From: jachym Date: Tue, 4 Mar 2025 17:02:45 +0100 Subject: [PATCH 5/9] test file --- .../generalized_normal_test.hpp | 168 ++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 test/prob/generalized_normal/generalized_normal_test.hpp diff --git a/test/prob/generalized_normal/generalized_normal_test.hpp b/test/prob/generalized_normal/generalized_normal_test.hpp new file mode 100644 index 00000000000..bfa98b97d14 --- /dev/null +++ b/test/prob/generalized_normal/generalized_normal_test.hpp @@ -0,0 +1,168 @@ +// Arguments: Doubles, Doubles, Doubles, Doubles +#include +#include +#include +#include +#include +#include + + +using std::numeric_limits; +using std::vector; + +class AgradDistributionGeneralizedNormal : public AgradDistributionTest { + public: + void valid_values(vector >& parameters, + vector& log_prob) { + vector param(4); + + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta + parameters.push_back(param); + log_prob.push_back(-0.57236494292470008707171367567652935582); // expected log_prob + + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta + parameters.push_back(param); + log_prob.push_back(-1.5723649429247000870717136756765293558); // expected log_prob + + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta + parameters.push_back(param); + log_prob.push_back(-4.5723649429247000870717136756765293558); // expected log_prob + + param[0] = -3.5; // y + param[1] = 1.9; // mu + param[2] = 7.2; // alpha + param[3] = 2; // beta + parameters.push_back(param); + log_prob.push_back(-3.1089459689467097140959090592117462617); // expected log_prob + + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta + parameters.push_back(param); + log_prob.push_back(-0.69314718055994530941723212145817656808); // expected log_prob + + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta + parameters.push_back(param); + log_prob.push_back(-1.6931471805599453094172321214581765681); // expected log_prob + + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta + parameters.push_back(param); + log_prob.push_back(-2.6931471805599453094172321214581765681); // expected log_prob + + param[0] = -3.5; // y + param[1] = 1.9; // mu + param[2] = 7.2; // alpha + param[3] = 1; // beta + parameters.push_back(param); + log_prob.push_back(-3.4172282065819549364414275049933934740); // expected log_prob + + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta + parameters.push_back(param); + log_prob.push_back(-0.59083234759930449611508182336583846717); // expected log_prob + + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta + parameters.push_back(param); + log_prob.push_back(-1.5908323475993044961150818233658384672); // expected log_prob + + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta + parameters.push_back(param); + log_prob.push_back(-3.4192594723454945937184592717852346243); // expected log_prob + + param[0] = -3.5; // y + param[1] = 1.9; // mu + param[2] = 7.2; // alpha + param[3] = 1.5; // beta + parameters.push_back(param); + log_prob.push_back(-3.2144324264596431082120695849657575107); // expected log_prob + + } + + void invalid_values(vector& index, vector& value) { + // y + + // mu + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + // alpha + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(-numeric_limits::infinity()); + + // beta + index.push_back(3U); + value.push_back(0.0); + + index.push_back(3U); + value.push_back(-1.0); + + index.push_back(3U); + value.push_back(-numeric_limits::infinity()); + } + + template + stan::return_type_t + log_prob(const T_y& y, const T_loc& mu, const T_scale& alpha, + const T_shape& beta, const T5&, const T6&) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + } + + template + stan::return_type_t + log_prob(const T_y& y, const T_loc& mu, const T_scale& alpha, + const T_shape& beta, const T5&, const T6&) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + } + + template + stan::return_type_t log_prob_function( + const T_y& y, const T_loc& mu, const T_scale& alpha, const T_shape& beta, + const T5&, const T6&) { + using stan::math::abs; + using stan::math::inv; + using stan::math::lgamma; + using stan::math::log; + using stan::math::LOG_TWO; + + return - LOG_TWO + - log(alpha) + - lgamma(1.0 + inv(beta)) + - pow(abs(y - mu) / alpha, beta); + } +}; From 80c4d61942e8d0680fdcae3d1b4e0669978b8095 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 4 Mar 2025 11:06:27 -0500 Subject: [PATCH 6/9] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../prim/prob/double_exponential_lpdf.hpp | 4 +- .../prim/prob/generalized_normal_lpdf.hpp | 58 +++++--- stan/math/prim/prob/normal_lpdf.hpp | 8 +- .../generalized_normal_test.hpp | 126 ++++++++++-------- 4 files changed, 110 insertions(+), 86 deletions(-) diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp index 2169a27afda..aaad44850aa 100644 --- a/stan/math/prim/prob/double_exponential_lpdf.hpp +++ b/stan/math/prim/prob/double_exponential_lpdf.hpp @@ -84,8 +84,8 @@ return_type_t double_exponential_lpdf( if (!is_constant_all::value) { const auto& rep_deriv - = to_ref_if<(!is_constant::value - && !is_constant::value)>(sign(y_m_mu) * inv_sigma); + = to_ref_if<(!is_constant::value && !is_constant::value)>( + sign(y_m_mu) * inv_sigma); if (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index 7ffbadcc6fb..2c478c92aca 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -27,7 +27,8 @@ namespace math { /** \ingroup prob_dists * The log of the generalized normal density for the specified scalar(s) given * the specified location, scale and shape parameters. y, mu, alpha, or beta can - * each be either a scalar or a vector. Any vector inputs must be the same length. + * each be either a scalar or a vector. Any vector inputs must be the same + * length. * *

The result log probability is defined to be the sum of the * log probabilities for each observation/mean/scale/shape tuple. @@ -55,8 +56,8 @@ inline return_type_t generalized_normal_lpdf( using T_alpha_ref = ref_type_if_not_constant_t; using T_beta_ref = ref_type_if_not_constant_t; static constexpr const char* function = "generalized_normal_lpdf"; - check_consistent_sizes(function, "Random variable", y, "Location parameter", mu, - "Scale parameter", alpha, "Shape parameter", beta); + check_consistent_sizes(function, "Random variable", y, "Location parameter", + mu, "Scale parameter", alpha, "Shape parameter", beta); T_y_ref y_ref = std::forward(y); T_mu_ref mu_ref = std::forward(mu); @@ -71,7 +72,9 @@ inline return_type_t generalized_normal_lpdf( check_not_nan(function, "Random variable", y_val); check_finite(function, "Location parameter", mu_val); check_positive(function, "Scale parameter", alpha_val); - check_positive(function, "Shape parameter", beta_val); // With β = +∞ this could be defined to be uniform, but we don't support that. + check_positive(function, "Shape parameter", + beta_val); // With β = +∞ this could be defined to be uniform, + // but we don't support that. if (size_zero(y, mu, alpha, beta)) { return 0; @@ -80,11 +83,17 @@ inline return_type_t generalized_normal_lpdf( return 0; } - const auto& inv_beta1p = to_ref_if::value>(inv(beta_val) + 1); - const auto& diff = to_ref_if::value>(y_val - mu_val); + const auto& inv_beta1p + = to_ref_if::value>(inv(beta_val) + 1); + const auto& diff + = to_ref_if::value>(y_val - mu_val); const auto& inv_alpha = to_ref(inv(alpha_val)); - const auto& scaled_abs_diff = to_ref_if::value>(abs(diff) * inv_alpha); - const auto& scaled_abs_diff_pow = to_ref_if::value>(pow(scaled_abs_diff, beta_val)); + const auto& scaled_abs_diff + = to_ref_if::value>(abs(diff) + * inv_alpha); + const auto& scaled_abs_diff_pow + = to_ref_if::value>( + pow(scaled_abs_diff, beta_val)); const size_t N = max_size(y, mu, alpha, beta); T_partials_return logp = -sum(scaled_abs_diff_pow); @@ -99,16 +108,21 @@ inline return_type_t generalized_normal_lpdf( logp -= sum(lgamma(inv_beta1p)) * (N / math::size(beta)); } - auto ops_partials = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); + auto ops_partials + = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); if (!is_constant_all::value) { - // note: The partial derivatives for y, mu are undefined when y == mu && beta < 1. - // The derivative limit as mu -> y goes: + // note: The partial derivatives for y, mu are undefined when y == mu && + // beta < 1. The derivative limit as mu -> y goes: // to 0 from both sides if beta > 1 (defined as 0) - // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as 0, consistent with double_exponential_lpdf) - // to +∞ from right but -∞ from left as y -> mu if beta < 1 (undefined) - const auto& rep_deriv = to_ref_if::value - && !is_constant::value>(sign(diff) * beta_val * pow(scaled_abs_diff, beta_val-1) * inv_alpha); + // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as + // 0, consistent with double_exponential_lpdf) to +∞ from right but -∞ + // from left as y -> mu if beta < 1 (undefined) + const auto& rep_deriv + = to_ref_if < !is_constant::value + && !is_constant::value + > (sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) + * inv_alpha); if (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; } @@ -117,10 +131,13 @@ inline return_type_t generalized_normal_lpdf( } } if (!is_constant::value) { - partials<2>(ops_partials) = (beta_val * scaled_abs_diff_pow - 1) * inv_alpha; + partials<2>(ops_partials) + = (beta_val * scaled_abs_diff_pow - 1) * inv_alpha; } if (!is_constant::value) { - partials<3>(ops_partials) = digamma(inv_beta1p) * inv_square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); + partials<3>(ops_partials) + = digamma(inv_beta1p) * inv_square(beta_val) + - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); } return ops_partials.build(logp); @@ -129,10 +146,9 @@ inline return_type_t generalized_normal_lpdf( template inline return_type_t generalized_normal_lpdf( T_y&& y, T_loc&& mu, T_scale&& alpha, T_shape&& beta) { - return generalized_normal_lpdf(std::forward(y), - std::forward(mu), - std::forward(alpha), - std::forward(beta)); + return generalized_normal_lpdf( + std::forward(y), std::forward(mu), + std::forward(alpha), std::forward(beta)); } } // namespace math diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index 83242104591..bb4b7eb9429 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -85,10 +85,10 @@ inline return_type_t normal_lpdf(T_y&& y, T_loc&& mu, } if (!is_constant_all::value) { - auto scaled_diff = to_ref_if::value - + !is_constant::value - + !is_constant::value - >= 2>(inv_sigma * y_scaled); + auto scaled_diff + = to_ref_if::value + !is_constant::value + + !is_constant::value + >= 2>(inv_sigma * y_scaled); if (!is_constant::value) { partials<0>(ops_partials) = -scaled_diff; } diff --git a/test/prob/generalized_normal/generalized_normal_test.hpp b/test/prob/generalized_normal/generalized_normal_test.hpp index bfa98b97d14..ee26e3eb486 100644 --- a/test/prob/generalized_normal/generalized_normal_test.hpp +++ b/test/prob/generalized_normal/generalized_normal_test.hpp @@ -6,7 +6,6 @@ #include #include - using std::numeric_limits; using std::vector; @@ -16,90 +15,101 @@ class AgradDistributionGeneralizedNormal : public AgradDistributionTest { vector& log_prob) { vector param(4); - param[0] = 0; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 2; // beta + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta parameters.push_back(param); - log_prob.push_back(-0.57236494292470008707171367567652935582); // expected log_prob + log_prob.push_back( + -0.57236494292470008707171367567652935582); // expected log_prob - param[0] = 1; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 2; // beta + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta parameters.push_back(param); - log_prob.push_back(-1.5723649429247000870717136756765293558); // expected log_prob + log_prob.push_back( + -1.5723649429247000870717136756765293558); // expected log_prob - param[0] = -2; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 2; // beta + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 2; // beta parameters.push_back(param); - log_prob.push_back(-4.5723649429247000870717136756765293558); // expected log_prob + log_prob.push_back( + -4.5723649429247000870717136756765293558); // expected log_prob param[0] = -3.5; // y param[1] = 1.9; // mu param[2] = 7.2; // alpha param[3] = 2; // beta parameters.push_back(param); - log_prob.push_back(-3.1089459689467097140959090592117462617); // expected log_prob + log_prob.push_back( + -3.1089459689467097140959090592117462617); // expected log_prob - param[0] = 0; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1; // beta + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta parameters.push_back(param); - log_prob.push_back(-0.69314718055994530941723212145817656808); // expected log_prob + log_prob.push_back( + -0.69314718055994530941723212145817656808); // expected log_prob - param[0] = 1; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1; // beta + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta parameters.push_back(param); - log_prob.push_back(-1.6931471805599453094172321214581765681); // expected log_prob + log_prob.push_back( + -1.6931471805599453094172321214581765681); // expected log_prob - param[0] = -2; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1; // beta + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta parameters.push_back(param); - log_prob.push_back(-2.6931471805599453094172321214581765681); // expected log_prob + log_prob.push_back( + -2.6931471805599453094172321214581765681); // expected log_prob param[0] = -3.5; // y param[1] = 1.9; // mu param[2] = 7.2; // alpha param[3] = 1; // beta parameters.push_back(param); - log_prob.push_back(-3.4172282065819549364414275049933934740); // expected log_prob + log_prob.push_back( + -3.4172282065819549364414275049933934740); // expected log_prob - param[0] = 0; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1.5; // beta + param[0] = 0; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta parameters.push_back(param); - log_prob.push_back(-0.59083234759930449611508182336583846717); // expected log_prob + log_prob.push_back( + -0.59083234759930449611508182336583846717); // expected log_prob - param[0] = 1; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1.5; // beta + param[0] = 1; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta parameters.push_back(param); - log_prob.push_back(-1.5908323475993044961150818233658384672); // expected log_prob + log_prob.push_back( + -1.5908323475993044961150818233658384672); // expected log_prob - param[0] = -2; // y - param[1] = 0; // mu - param[2] = 1; // alpha - param[3] = 1.5; // beta + param[0] = -2; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta parameters.push_back(param); - log_prob.push_back(-3.4192594723454945937184592717852346243); // expected log_prob + log_prob.push_back( + -3.4192594723454945937184592717852346243); // expected log_prob param[0] = -3.5; // y param[1] = 1.9; // mu param[2] = 7.2; // alpha param[3] = 1.5; // beta parameters.push_back(param); - log_prob.push_back(-3.2144324264596431082120695849657575107); // expected log_prob - + log_prob.push_back( + -3.2144324264596431082120695849657575107); // expected log_prob } void invalid_values(vector& index, vector& value) { @@ -135,17 +145,17 @@ class AgradDistributionGeneralizedNormal : public AgradDistributionTest { template - stan::return_type_t - log_prob(const T_y& y, const T_loc& mu, const T_scale& alpha, - const T_shape& beta, const T5&, const T6&) { + stan::return_type_t log_prob( + const T_y& y, const T_loc& mu, const T_scale& alpha, const T_shape& beta, + const T5&, const T6&) { return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); } template - stan::return_type_t - log_prob(const T_y& y, const T_loc& mu, const T_scale& alpha, - const T_shape& beta, const T5&, const T6&) { + stan::return_type_t log_prob( + const T_y& y, const T_loc& mu, const T_scale& alpha, const T_shape& beta, + const T5&, const T6&) { return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); } @@ -160,9 +170,7 @@ class AgradDistributionGeneralizedNormal : public AgradDistributionTest { using stan::math::log; using stan::math::LOG_TWO; - return - LOG_TWO - - log(alpha) - - lgamma(1.0 + inv(beta)) + return -LOG_TWO - log(alpha) - lgamma(1.0 + inv(beta)) - pow(abs(y - mu) / alpha, beta); } }; From 3911299b448bc110d82dabad92da57886185b577 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Wed, 5 Mar 2025 12:19:07 +0100 Subject: [PATCH 7/9] restore develop state for changes in normal and double_exponential --- .../prim/prob/double_exponential_lpdf.hpp | 28 +++++++++++-------- stan/math/prim/prob/normal_lpdf.hpp | 16 ++++++----- .../double_exponential_test.hpp | 1 + test/prob/normal/normal_test.hpp | 1 + 4 files changed, 27 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp index aaad44850aa..059eebdeb4e 100644 --- a/stan/math/prim/prob/double_exponential_lpdf.hpp +++ b/stan/math/prim/prob/double_exponential_lpdf.hpp @@ -3,9 +3,11 @@ #include #include +#include +#include #include #include -#include +#include #include #include #include @@ -57,6 +59,7 @@ return_type_t double_exponential_lpdf( return 0.0; } + T_partials_return logp(0.0); auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); @@ -70,32 +73,33 @@ return_type_t double_exponential_lpdf( const auto& inv_sigma = to_ref(inv(sigma_val)); const auto& y_m_mu = to_ref_if::value>(y_val - mu_val); - const auto& scaled_abs_diff - = to_ref_if::value>(abs(y_m_mu) * inv_sigma); + const auto& abs_diff_y_mu = fabs(y_m_mu); + const auto& scaled_diff + = to_ref_if::value>(abs_diff_y_mu * inv_sigma); - T_partials_return logp = -sum(scaled_abs_diff); - const size_t N = max_size(y, mu, sigma); + size_t N = max_size(y, mu, sigma); if (include_summand::value) { logp -= N * LOG_TWO; } if (include_summand::value) { logp -= sum(log(sigma_val)) * N / math::size(sigma); } + logp -= sum(scaled_diff); if (!is_constant_all::value) { + const auto& diff_sign = sign(y_m_mu); const auto& rep_deriv - = to_ref_if<(!is_constant::value && !is_constant::value)>( - sign(y_m_mu) * inv_sigma); - - if (!is_constant::value) { + = to_ref_if<(!is_constant_all::value + && !is_constant_all::value)>(diff_sign * inv_sigma); + if (!is_constant_all::value) { partials<0>(ops_partials) = -rep_deriv; } - if (!is_constant::value) { + if (!is_constant_all::value) { partials<1>(ops_partials) = rep_deriv; } } - if (!is_constant::value) { - partials<2>(ops_partials) = inv_sigma * (scaled_abs_diff - 1); + if (!is_constant_all::value) { + partials<2>(ops_partials) = inv_sigma * (scaled_diff - 1); } return ops_partials.build(logp); diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index bb4b7eb9429..f53c521a74e 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -3,6 +3,8 @@ #include #include +#include +#include #include #include #include @@ -85,17 +87,17 @@ inline return_type_t normal_lpdf(T_y&& y, T_loc&& mu, } if (!is_constant_all::value) { - auto scaled_diff - = to_ref_if::value + !is_constant::value - + !is_constant::value - >= 2>(inv_sigma * y_scaled); - if (!is_constant::value) { + auto scaled_diff = to_ref_if::value + + !is_constant_all::value + + !is_constant_all::value + >= 2>(inv_sigma * y_scaled); + if (!is_constant_all::value) { partials<0>(ops_partials) = -scaled_diff; } - if (!is_constant::value) { + if (!is_constant_all::value) { partials<2>(ops_partials) = inv_sigma * y_scaled_sq - inv_sigma; } - if (!is_constant::value) { + if (!is_constant_all::value) { partials<1>(ops_partials) = std::move(scaled_diff); } } diff --git a/test/prob/double_exponential/double_exponential_test.hpp b/test/prob/double_exponential/double_exponential_test.hpp index be206376758..db5aa235a59 100644 --- a/test/prob/double_exponential/double_exponential_test.hpp +++ b/test/prob/double_exponential/double_exponential_test.hpp @@ -2,6 +2,7 @@ #include #include +using stan::math::var; using std::numeric_limits; using std::vector; diff --git a/test/prob/normal/normal_test.hpp b/test/prob/normal/normal_test.hpp index df818f10e02..81e2ca75927 100644 --- a/test/prob/normal/normal_test.hpp +++ b/test/prob/normal/normal_test.hpp @@ -3,6 +3,7 @@ #include #include +using stan::math::var; using std::numeric_limits; using std::vector; From 800eb41cc33e04093bee307aadee46e71d8bd00f Mon Sep 17 00:00:00 2001 From: jachym Date: Sat, 8 Mar 2025 00:06:24 +0100 Subject: [PATCH 8/9] Incorporate changes requested by @SteveBronder --- .../prim/prob/generalized_normal_lpdf.hpp | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index 7ffbadcc6fb..d993deb5d01 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -76,7 +76,7 @@ inline return_type_t generalized_normal_lpdf( if (size_zero(y, mu, alpha, beta)) { return 0; } - if (!include_summand::value) { + if constexpr (!include_summand::value) { return 0; } @@ -89,37 +89,38 @@ inline return_type_t generalized_normal_lpdf( T_partials_return logp = -sum(scaled_abs_diff_pow); - if (include_summand::value) { + if constexpr (include_summand::value) { logp -= LOG_TWO * N; } - if (include_summand::value) { + if constexpr (include_summand::value) { logp -= sum(log(alpha_val)) * (N / math::size(alpha)); } - if (include_summand::value) { + if constexpr (include_summand::value) { logp -= sum(lgamma(inv_beta1p)) * (N / math::size(beta)); } auto ops_partials = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); - if (!is_constant_all::value) { - // note: The partial derivatives for y, mu are undefined when y == mu && beta < 1. - // The derivative limit as mu -> y goes: - // to 0 from both sides if beta > 1 (defined as 0) - // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as 0, consistent with double_exponential_lpdf) - // to +∞ from right but -∞ from left as y -> mu if beta < 1 (undefined) - const auto& rep_deriv = to_ref_if::value - && !is_constant::value>(sign(diff) * beta_val * pow(scaled_abs_diff, beta_val-1) * inv_alpha); - if (!is_constant::value) { + if constexpr (!is_constant_all::value) { + // note: The partial derivatives for y, μ are undefined when + // y == μ && beta < 1. + // The derivative limit as y → μ (i.e. diff → 0) has the following cases: + // β > 1: 0 from both sides (defined as 0) + // β == 1: +1/α from right, but -1/α from left (defined as + // 0, consistent with double_exponential_lpdf) + // β < 1: -∞ from left as y → μ, but +∞ from right (undefined) + auto rep_deriv = eval(sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) * inv_alpha); + if constexpr (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; } - if (!is_constant::value) { + if constexpr (!is_constant::value) { partials<1>(ops_partials) = std::move(rep_deriv); } } - if (!is_constant::value) { + if constexpr (!is_constant::value) { partials<2>(ops_partials) = (beta_val * scaled_abs_diff_pow - 1) * inv_alpha; } - if (!is_constant::value) { + if constexpr (!is_constant::value) { partials<3>(ops_partials) = digamma(inv_beta1p) * inv_square(beta_val) - multiply_log(scaled_abs_diff_pow, scaled_abs_diff); } From 4dab04576995561e0c8f3aacf0b01c8549dc56b6 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 7 Mar 2025 18:24:10 -0500 Subject: [PATCH 9/9] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/generalized_normal_lpdf.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/prob/generalized_normal_lpdf.hpp b/stan/math/prim/prob/generalized_normal_lpdf.hpp index 5cd73cf4369..919597d9418 100644 --- a/stan/math/prim/prob/generalized_normal_lpdf.hpp +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -74,7 +74,7 @@ inline return_type_t generalized_normal_lpdf( check_positive(function, "Scale parameter", alpha_val); // With β = +∞ this could be defined to be uniform, but we don't support that. - check_positive(function, "Shape parameter", beta_val); + check_positive(function, "Shape parameter", beta_val); if (size_zero(y, mu, alpha, beta)) { return 0; @@ -112,16 +112,15 @@ inline return_type_t generalized_normal_lpdf( = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); if constexpr (!is_constant_all::value) { - // note: The partial derivatives for y, μ are undefined when - // y == μ && beta < 1. + // note: The partial derivatives for y, μ are undefined when + // y == μ && beta < 1. // The derivative limit as y → μ (i.e. diff → 0) has the following cases: // β > 1: 0 from both sides (defined as 0) // β == 1: +1/α from right, but -1/α from left (defined as - // 0, consistent with double_exponential_lpdf) + // 0, consistent with double_exponential_lpdf) // β < 1: -∞ from left as y → μ, but +∞ from right (undefined) - auto rep_deriv = eval( - sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) * inv_alpha - ); + auto rep_deriv = eval(sign(diff) * beta_val + * pow(scaled_abs_diff, beta_val - 1) * inv_alpha); if constexpr (!is_constant::value) { partials<0>(ops_partials) = -rep_deriv; }