diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 308b084a165..132dd2b6fa1 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -113,6 +113,7 @@ #include #include #include +#include #include #include #include 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..f7ad3739543 --- /dev/null +++ b/stan/math/prim/prob/generalized_normal_lpdf.hpp @@ -0,0 +1,195 @@ +#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 +#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 or not finite + */ +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_finite(function, "Scale parameter", alpha_val); + + // With β = +∞ this could be defined to be uniform, but we don't support that. + check_positive_finite(function, "Shape parameter", beta_val); + + if (size_zero(y, mu, alpha, beta)) { + return 0; + } + if constexpr (!include_summand::value) { + return 0; + } + + const auto& one_plus_inv_beta + = to_ref_if::value>(inv(beta_val) + 1); + const auto& residual + = to_ref_if::value>(y_val - mu_val); + const auto& inverse_scale + = to_ref_if::value>(inv(alpha_val)); + const auto& scaled_residual = to_ref(residual / alpha_val); + const auto& abs_scaled_residual = to_ref(abs(scaled_residual)); + const auto& scaled_residual_is_zero + = eval(value_of_rec(abs_scaled_residual) == 0); + auto abs_scaled_residual_or_one = eval(abs_scaled_residual); + if constexpr (is_eigen_v>) { + abs_scaled_residual_or_one + = scaled_residual_is_zero.select(1.0, abs_scaled_residual_or_one); + } else { + if (scaled_residual_is_zero) { + abs_scaled_residual_or_one += 1; + } + } + auto abs_scaled_residual_pow_beta + = eval(pow(abs_scaled_residual_or_one, beta_val)); + // At beta == 2 the kernel is exactly square(scaled_residual). Retaining that + // expression at the tie preserves the normal Hessian; differentiating the + // generic abs/pow expression would instead inherit abs'(0) == 0. + auto beta_is_two = eval(value_of_rec(beta_val) == 2.0); + auto apply_zero_residual_branch = [&](auto& generic_value, + const auto& beta_two_value) { + if constexpr (is_eigen_v>) { + if constexpr (is_eigen_v>) { + generic_value = scaled_residual_is_zero.select( + beta_is_two.select(beta_two_value, 0.0), generic_value); + } else if (beta_is_two) { + generic_value + = scaled_residual_is_zero.select(beta_two_value, generic_value); + } else { + generic_value = scaled_residual_is_zero.select(0.0, generic_value); + } + } else if (scaled_residual_is_zero) { + if constexpr (is_eigen_v>) { + generic_value = beta_is_two.select(beta_two_value, 0.0); + } else { + generic_value = beta_is_two ? beta_two_value : 0; + } + } + }; + auto beta_two_kernel = eval(square(scaled_residual) + 0.0 * beta_val); + apply_zero_residual_branch(abs_scaled_residual_pow_beta, beta_two_kernel); + const size_t num_terms = max_size(y, mu, alpha, beta); + + T_partials_return logp = -sum(abs_scaled_residual_pow_beta); + + if constexpr (include_summand::value) { + logp -= LOG_TWO * num_terms; + } + if constexpr (include_summand::value) { + logp -= sum(log(alpha_val)) * (num_terms / math::size(alpha)); + } + if constexpr (include_summand::value) { + logp -= sum(lgamma(one_plus_inv_beta)) * (num_terms / math::size(beta)); + } + + auto ops_partials + = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); + + if constexpr (!is_constant_all::value) { + // At y == mu, the derivative is 0 for beta > 1 and undefined otherwise. + // Use abs_scaled_residual_or_one (0 replaced with 1) to avoid NaN from + // pow(0, beta-1), and preserve the exact normal curvature when beta == 2. + auto kernel_derivative_wrt_residual + = eval(sign(residual) * beta_val + * pow(abs_scaled_residual_or_one, beta_val - 1) * inverse_scale); + const auto& beta_two_kernel_derivative + = eval(2.0 * scaled_residual / alpha_val + 0.0 * beta_val); + apply_zero_residual_branch(kernel_derivative_wrt_residual, + beta_two_kernel_derivative); + if constexpr (!is_constant::value) { + partials<0>(ops_partials) = -kernel_derivative_wrt_residual; + } + if constexpr (!is_constant::value) { + partials<1>(ops_partials) = std::move(kernel_derivative_wrt_residual); + } + } + if constexpr (!is_constant::value) { + partials<2>(ops_partials) + = (beta_val * abs_scaled_residual_pow_beta - 1) * inverse_scale; + } + if constexpr (!is_constant::value) { + // multiply_log(0, 0) = 0 by convention, but fvar autodiff of + // multiply_log at (0, 0) produces NaN. abs_scaled_residual_or_one avoids + // this. + partials<3>(ops_partials) + = digamma(one_plus_inv_beta) * inv_square(beta_val) + - multiply_log(abs_scaled_residual_pow_beta, + abs_scaled_residual_or_one); + } + + 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/test/prob/generalized_normal/generalized_normal_test.hpp b/test/prob/generalized_normal/generalized_normal_test.hpp new file mode 100644 index 00000000000..fd17690ce0e --- /dev/null +++ b/test/prob/generalized_normal/generalized_normal_test.hpp @@ -0,0 +1,228 @@ +// 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 + + param[0] = 0.5; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1; // beta + parameters.push_back(param); + log_prob.push_back( + -1.1931471805599453094172321214581765681); // expected log_prob + + param[0] = 0.5; // y + param[1] = 0; // mu + param[2] = 1; // alpha + param[3] = 1.5; // beta + parameters.push_back(param); + log_prob.push_back( + -0.94438573819257864983001127257011830807); // 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; + + auto base = abs(y - mu) / alpha; + bool at_zero = stan::math::value_of_rec(base) == 0; + if (at_zero) + base += 1; + auto pow_term = pow(base, beta); + if (at_zero) + pow_term = 0; + return -LOG_TWO - log(alpha) - lgamma(1.0 + inv(beta)) - pow_term; + } +}; + +TEST(ProbDistributionsGeneralizedNormal, VectorWithYEqualsMu) { + using Eigen::VectorXd; + using stan::math::generalized_normal_lpdf; + using stan::math::var; + + // y[1] == mu[1], other elements differ + VectorXd y(3), mu(3); + y << -1.0, 0.5, 2.0; + mu << 0.0, 0.5, 0.0; + double alpha = 1.0; + double beta = 1.5; + + double lp = generalized_normal_lpdf(y, mu, alpha, beta); + EXPECT_TRUE(std::isfinite(lp)); + + // Same with var types to check autodiff + std::vector y_v(y.data(), y.data() + y.size()); + std::vector mu_v(mu.data(), mu.data() + mu.size()); + var alpha_v = alpha; + var beta_v = beta; + var lp_v = generalized_normal_lpdf(y_v, mu_v, alpha_v, beta_v); + lp_v.grad(); + for (size_t i = 0; i < y_v.size(); ++i) { + EXPECT_TRUE(std::isfinite(y_v[i].adj())); + EXPECT_TRUE(std::isfinite(mu_v[i].adj())); + } + EXPECT_TRUE(std::isfinite(alpha_v.adj())); + EXPECT_TRUE(std::isfinite(beta_v.adj())); +} diff --git a/test/unit/math/mix/prob/generalized_normal_lpdf_test.cpp b/test/unit/math/mix/prob/generalized_normal_lpdf_test.cpp new file mode 100644 index 00000000000..75b1b6829b0 --- /dev/null +++ b/test/unit/math/mix/prob/generalized_normal_lpdf_test.cpp @@ -0,0 +1,340 @@ +#include +#include +#include +#include + +namespace generalized_normal_lpdf_test { + +// expect_ad takes at most three autodiff arguments, so each factory holds +// one of (y, mu, alpha, beta) primitive and rotates the other three. +auto f_y_mu_alpha(double beta) { + return [beta](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + }; +} + +auto f_y_mu_beta(double alpha) { + return [alpha](const auto& y, const auto& mu, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + }; +} + +auto f_y_alpha_beta(double mu) { + return [mu](const auto& y, const auto& alpha, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + }; +} + +auto f_mu_alpha_beta(double y) { + return [y](const auto& mu, const auto& alpha, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + }; +} + +// Every argument is under autodiff in at least one of the four rotations. +void expect_all_rotations(double y, double mu, double alpha, double beta) { + stan::test::expect_ad(f_y_mu_alpha(beta), y, mu, alpha); + stan::test::expect_ad(f_y_mu_beta(alpha), y, mu, beta); + stan::test::expect_ad(f_y_alpha_beta(mu), y, alpha, beta); + stan::test::expect_ad(f_mu_alpha_beta(y), mu, alpha, beta); +} + +} // namespace generalized_normal_lpdf_test + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_interior) { + using generalized_normal_lpdf_test::expect_all_rotations; + + expect_all_rotations(1.3, 0.4, 1.7, 2.0); + expect_all_rotations(-2.5, 0.8, 0.6, 3.0); + expect_all_rotations(0.0, -1.1, 2.3, 1.4); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_beta_boundaries) { + using generalized_normal_lpdf_test::expect_all_rotations; + + // beta == 1 is the double exponential, beta == 2 the normal. + expect_all_rotations(1.7, 0.3, 1.2, 1.0); + expect_all_rotations(1.7, 0.3, 1.2, 2.0); + + // beta < 1 gives a cusp at y == mu; beta large approaches the uniform. + expect_all_rotations(1.7, 0.3, 1.2, 0.5); + expect_all_rotations(0.4, 0.3, 1.2, 8.0); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_alpha_boundaries) { + using generalized_normal_lpdf_test::expect_all_rotations; + + expect_all_rotations(0.9, 0.5, 0.05, 2.0); + expect_all_rotations(0.9, 0.5, 40.0, 2.0); + expect_all_rotations(0.9, 0.5, 40.0, 0.7); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_tails) { + using generalized_normal_lpdf_test::expect_all_rotations; + + // Large (|y - mu| / alpha) ^ beta, where pow() dominates the value. + expect_all_rotations(12.0, 0.0, 1.0, 2.0); + expect_all_rotations(-9.5, 1.5, 0.8, 1.0); + expect_all_rotations(30.0, 0.0, 1.0, 0.5); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_y_equals_mu_smooth) { + using generalized_normal_lpdf_test::expect_all_rotations; + + // expect_ad checks up to third derivatives, and d3/dy3 log p only exists + // at y == mu for beta > 3, so the tie is only tested there. + expect_all_rotations(0.5, 0.5, 1.3, 4.0); + expect_all_rotations(-2.0, -2.0, 0.4, 4.5); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_y_equals_mu_normal) { + using generalized_normal_lpdf_test::f_y_mu_alpha; + + // With beta held fixed at 2 the distribution is normal(mu, + // alpha / sqrt(2)), so all y, mu, and alpha derivatives exist. + stan::test::expect_ad(f_y_mu_alpha(2.0), 0.5, 0.5, 1.3); + stan::test::expect_ad(f_y_mu_alpha(2.0), 0.0, 0.0, 1.0); + + // Same requirement inside a vectorized call, where only y[0] == mu[0]. + auto f_vec = [](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, 2.0); + }; + Eigen::VectorXd y(2), mu(2), alpha(2); + y << 0.5, 0.3; + mu << 0.5, -0.2; + alpha << 1.0, 1.4; + stan::test::expect_ad(f_vec, y, mu, alpha); + + // The same boundary handling applies elementwise when beta is a vector. + Eigen::VectorXd beta(2); + beta << 2.0, 3.5; + auto f_vec_beta = [&beta](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, beta); + }; + stan::test::expect_ad(f_vec_beta, y, mu, alpha); +} + +TEST_F(AgradRev, + mathMixScalFun_generalized_normal_lpdf_y_equals_mu_variable_beta) { + using stan::math::generalized_normal_lpdf; + + // Pure beta derivatives are smooth at the tie. In particular, the kernel + // contribution r^beta * log(r) is zero by continuity. + auto f_beta = [](const auto& beta) { + return generalized_normal_lpdf(0.5, 0.5, 1.3, beta); + }; + stan::test::expect_ad(f_beta, 2.0); + + // The joint value, gradient, and Hessian also exist at beta == 2. Do not use + // expect_ad on all four arguments here: it additionally checks the mixed + // third derivative d3/(d beta d y2), which diverges at y == mu. + auto f_joint = [](const auto& x) { + return generalized_normal_lpdf(x(0), x(1), x(2), x(3)); + }; + Eigen::VectorXd x(4); + x << 0.5, 0.5, 1.3, 2.0; + double fx; + Eigen::VectorXd grad; + Eigen::MatrixXd hessian; + stan::math::hessian(f_joint, x, fx, grad, hessian); + + const double alpha = x(2); + EXPECT_NEAR(-stan::math::LOG_TWO - std::log(alpha) - stan::math::lgamma(1.5), + fx, 1e-12); + EXPECT_NEAR(0.0, grad(0), 1e-12); + EXPECT_NEAR(0.0, grad(1), 1e-12); + EXPECT_NEAR(-1.0 / alpha, grad(2), 1e-12); + EXPECT_NEAR(stan::math::digamma(1.5) / 4.0, grad(3), 1e-12); + + const double normal_curvature = -2.0 / (alpha * alpha); + EXPECT_NEAR(normal_curvature, hessian(0, 0), 1e-12); + EXPECT_NEAR(-normal_curvature, hessian(0, 1), 1e-12); + EXPECT_NEAR(-normal_curvature, hessian(1, 0), 1e-12); + EXPECT_NEAR(normal_curvature, hessian(1, 1), 1e-12); + EXPECT_NEAR(0.0, hessian(0, 3), 1e-12); + EXPECT_NEAR(0.0, hessian(1, 3), 1e-12); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_y_equals_mu_cusp) { + using stan::math::generalized_normal_lpdf; + using stan::math::var; + + // At the tie, the first derivative is undefined for beta <= 1 and the + // second derivative is not finite for beta < 2. Check Stan's zero-gradient + // convention and the remaining parameter derivatives directly. + for (double beta : {0.5, 1.0, 1.5}) { + var y = 0.75; + var mu = 0.75; + var alpha = 1.4; + var beta_v = beta; + var lp = generalized_normal_lpdf(y, mu, alpha, beta_v); + + double expect_lp = -stan::math::LOG_TWO - std::log(1.4) + - stan::math::lgamma(1.0 + 1.0 / beta); + EXPECT_FLOAT_EQ(expect_lp, lp.val()); + + lp.grad(); + EXPECT_FLOAT_EQ(0.0, y.adj()); + EXPECT_FLOAT_EQ(0.0, mu.adj()); + EXPECT_FLOAT_EQ(-1.0 / 1.4, alpha.adj()); + EXPECT_FLOAT_EQ(stan::math::digamma(1.0 + 1.0 / beta) / (beta * beta), + beta_v.adj()); + stan::math::recover_memory(); + } +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_vectorized) { + auto f_vec = [](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, 3.0); + }; + auto f_vec_beta = [](const auto& y, const auto& mu, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, 1.1, beta); + }; + + Eigen::VectorXd y(3), mu(3), alpha(3), beta(3); + y << -1.0, 0.6, 2.0; + mu << 0.0, 0.5, -0.4; + alpha << 1.0, 2.5, 0.7; + beta << 3.5, 4.0, 5.5; + + stan::test::expect_ad(f_vec, y, mu, alpha); + stan::test::expect_ad(f_vec, y, mu[0], alpha); + stan::test::expect_ad(f_vec, y, mu, alpha[0]); + stan::test::expect_ad(f_vec, y[0], mu, alpha); + stan::test::expect_ad(f_vec, y[1], mu[0], alpha); + + stan::test::expect_ad(f_vec_beta, y, mu, beta); + stan::test::expect_ad(f_vec_beta, y, mu, beta[0]); + stan::test::expect_ad(f_vec_beta, y[1], mu[0], beta); + + std::vector y_std{-1.0, 0.6, 2.0}; + std::vector mu_std{0.0, 0.5, -0.4}; + std::vector alpha_std{1.0, 2.5, 0.7}; + stan::test::expect_ad(f_vec, y_std, mu_std, alpha_std); + stan::test::expect_ad(f_vec, y_std, mu_std[0], alpha_std); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_vectorized_tie) { + // y[1] == mu[1] exercises the y == mu branch inside a vectorized call; + // beta > 3 keeps the third derivative at the tie well defined. + auto f_vec = [](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, 4.0); + }; + auto f_vec_beta = [](const auto& y, const auto& mu, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, 1.1, beta); + }; + + Eigen::VectorXd y(3), mu(3), alpha(3), beta(3); + y << -1.0, 0.5, 2.0; + mu << 0.0, 0.5, -0.4; + alpha << 1.0, 2.5, 0.7; + beta << 3.5, 4.0, 5.5; + + stan::test::expect_ad(f_vec, y, mu, alpha); + stan::test::expect_ad(f_vec, y, mu, alpha[0]); + stan::test::expect_ad(f_vec, y[1], mu[1], alpha); + stan::test::expect_ad(f_vec_beta, y, mu, beta); + stan::test::expect_ad(f_vec_beta, y[1], mu[1], beta); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_propto) { + using stan::math::generalized_normal_lpdf; + using stan::math::var; + + // expect_ad cannot be used here: with all-double arguments propto=true + // drops every term and returns 0, so no value comparison is possible. + for (double y_val : {1.3, 0.4}) { + Eigen::Matrix x(4); + x << y_val, 0.4, 1.7, 2.5; + var lp_propto = generalized_normal_lpdf(x(0), x(1), x(2), x(3)); + lp_propto.grad(); + Eigen::VectorXd grad_propto(4); + for (int i = 0; i < 4; ++i) + grad_propto(i) = x(i).adj(); + double val_propto = lp_propto.val(); + stan::math::recover_memory(); + + Eigen::Matrix x2(4); + x2 << y_val, 0.4, 1.7, 2.5; + var lp_full = generalized_normal_lpdf(x2(0), x2(1), x2(2), x2(3)); + lp_full.grad(); + + // Only the -log(2) normalizing constant is dropped when propto=true. + EXPECT_FLOAT_EQ(lp_full.val() + stan::math::LOG_TWO, val_propto); + for (int i = 0; i < 4; ++i) + EXPECT_FLOAT_EQ(x2(i).adj(), grad_propto(i)); + stan::math::recover_memory(); + } +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_errors) { + using generalized_normal_lpdf_test::f_mu_alpha_beta; + using generalized_normal_lpdf_test::f_y_alpha_beta; + using generalized_normal_lpdf_test::f_y_mu_alpha; + using generalized_normal_lpdf_test::f_y_mu_beta; + + double inf = std::numeric_limits::infinity(); + double nan = std::numeric_limits::quiet_NaN(); + + // alpha must be positive, beta must be positive, mu must be finite, + // y must not be NaN -- every instantiation must throw where double does. + stan::test::expect_ad(f_y_mu_alpha(2.0), 1.0, 0.0, 0.0); + stan::test::expect_ad(f_y_mu_alpha(2.0), 1.0, 0.0, -1.5); + stan::test::expect_ad(f_y_mu_beta(1.0), 1.0, 0.0, 0.0); + stan::test::expect_ad(f_y_mu_beta(1.0), 1.0, 0.0, -2.0); + stan::test::expect_ad(f_y_alpha_beta(0.0), nan, 1.0, 2.0); + stan::test::expect_ad(f_mu_alpha_beta(1.0), inf, 1.0, 2.0); + stan::test::expect_ad(f_mu_alpha_beta(1.0), -inf, 1.0, 2.0); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_infinite_params) { + using stan::math::generalized_normal_lpdf; + double inf = std::numeric_limits::infinity(); + + // Infinite scale and shape parameters are outside the supported domain. + EXPECT_THROW(generalized_normal_lpdf(1.0, 0.0, inf, 2.0), std::domain_error); + EXPECT_THROW(generalized_normal_lpdf(1.0, 0.0, 1.0, inf), std::domain_error); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_subnormal_alpha) { + using stan::math::generalized_normal_lpdf; + + // inv(alpha) overflows for subnormal alpha, so abs(diff) * inv_alpha is + // 0 * inf = NaN at y == mu and the |y - mu| == 0 test never fires. + for (double a : {1e-309, 1e-310, std::numeric_limits::denorm_min()}) { + double expected + = -stan::math::LOG_TWO - std::log(a) - stan::math::lgamma(1.5); + EXPECT_FLOAT_EQ(expected, generalized_normal_lpdf(0.0, 0.0, a, 2.0)); + } +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_matvar) { + auto f = [](const auto& y, const auto& mu, const auto& alpha) { + return stan::math::generalized_normal_lpdf(y, mu, alpha, 3.0); + }; + auto f_beta = [](const auto& y, const auto& mu, const auto& beta) { + return stan::math::generalized_normal_lpdf(y, mu, 1.1, beta); + }; + + Eigen::VectorXd y(3), mu(3), alpha(3), beta(3); + y << -1.0, 0.6, 2.0; + mu << 0.0, 0.5, -0.4; + alpha << 1.0, 2.5, 0.7; + beta << 3.5, 4.0, 5.5; + + stan::test::expect_ad_matvar(f, y, mu, alpha); + stan::test::expect_ad_matvar(f, y, mu, alpha[0]); + stan::test::expect_ad_matvar(f_beta, y, mu, beta); +} + +TEST_F(AgradRev, mathMixScalFun_generalized_normal_lpdf_size_errors) { + using stan::math::generalized_normal_lpdf; + + Eigen::VectorXd y(3), mu(2); + y << 0.1, 0.2, 0.3; + mu << 0.0, 0.0; + EXPECT_THROW(generalized_normal_lpdf(y, mu, 1.0, 2.0), std::invalid_argument); + + Eigen::VectorXd empty(0); + EXPECT_FLOAT_EQ(0.0, generalized_normal_lpdf(empty, 0.0, 1.0, 2.0)); +}