diff --git a/stan/math/prim/prob/beta_cdf.hpp b/stan/math/prim/prob/beta_cdf.hpp index b099c8df6ab..2e51fbdf842 100644 --- a/stan/math/prim/prob/beta_cdf.hpp +++ b/stan/math/prim/prob/beta_cdf.hpp @@ -4,9 +4,8 @@ #include #include #include +#include #include -#include -#include #include #include #include @@ -59,9 +58,6 @@ inline return_type_t beta_cdf( scalar_seq_view y_vec(y_ref); scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); - size_t size_alpha = stan::math::size(alpha); - size_t size_beta = stan::math::size(beta); - size_t size_alpha_beta = max_size(alpha, beta); size_t N = max_size(y, alpha, beta); // Explicit return for extreme values @@ -72,33 +68,10 @@ inline return_type_t beta_cdf( } } - VectorBuilder, T_partials_return, T_scale_succ> - digamma_alpha(size_alpha); - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < size_alpha; n++) { - digamma_alpha[n] = digamma(alpha_vec.val(n)); - } - } - - VectorBuilder, T_partials_return, T_scale_fail> - digamma_beta(size_beta); - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < size_beta; n++) { - digamma_beta[n] = digamma(beta_vec.val(n)); - } - } - - VectorBuilder, - T_partials_return, T_scale_succ, T_scale_fail> - digamma_sum(size_alpha_beta); - if constexpr (is_any_autodiff_v) { - for (size_t n = 0; n < size_alpha_beta; n++) { - digamma_sum[n] = digamma(alpha_vec.val(n) + beta_vec.val(n)); - } - } - + constexpr bool calc_alpha = is_autodiff_v; + constexpr bool calc_beta = is_autodiff_v; for (size_t n = 0; n < N; n++) { - const T_partials_return y_dbl = y_vec.val(n); + const auto y_dbl = y_vec.val(n); // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -106,30 +79,37 @@ inline return_type_t beta_cdf( continue; } - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - const T_partials_return Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl); - const T_partials_return inv_Pn - = is_constant_all::value ? 0 : inv(Pn); + const auto alpha_dbl = alpha_vec.val(n); + const auto beta_dbl = beta_vec.val(n); + const auto Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl); P *= Pn; - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += inc_beta_ddz(alpha_dbl, beta_dbl, y_dbl) * inv_Pn; - } - - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += inc_beta_dda(alpha_dbl, beta_dbl, y_dbl, digamma_alpha[n], - digamma_sum[n]) - * inv_Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += inc_beta_ddb(alpha_dbl, beta_dbl, y_dbl, digamma_beta[n], - digamma_sum[n]) - * inv_Pn; + if constexpr (is_any_autodiff_v) { + const auto inv_Pn = inv(Pn); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] + += inc_beta_ddz(alpha_dbl, beta_dbl, y_dbl) * inv_Pn; + } + + if constexpr (calc_alpha && calc_beta) { + const auto [g1, g2] = grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0); + partials<1>(ops_partials)[n] += g1 * inv_Pn; + partials<2>(ops_partials)[n] += g2 * inv_Pn; + } else if constexpr (calc_alpha) { + partials<1>(ops_partials)[n] + += grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), 0.0, + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; + } else if constexpr (calc_beta) { + partials<2>(ops_partials)[n] + += grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, 0.0, digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; + } } } diff --git a/stan/math/prim/prob/beta_lccdf.hpp b/stan/math/prim/prob/beta_lccdf.hpp index 66232243893..ed4234a6f5a 100644 --- a/stan/math/prim/prob/beta_lccdf.hpp +++ b/stan/math/prim/prob/beta_lccdf.hpp @@ -67,40 +67,17 @@ inline return_type_t beta_lccdf( scalar_seq_view y_vec(y_ref); scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); - size_t size_alpha = stan::math::size(alpha); - size_t size_beta = stan::math::size(beta_param); - size_t size_alpha_beta = max_size(alpha, beta_param); size_t N = max_size(y, alpha, beta_param); - VectorBuilder, - T_partials_return, T_scale_succ> - digamma_alpha(size_alpha); - VectorBuilder, - T_partials_return, T_scale_fail> - digamma_beta(size_beta); - VectorBuilder, - T_partials_return, T_scale_succ, T_scale_fail> - digamma_sum(size_alpha_beta); - - if constexpr (is_any_autodiff_v) { - for (size_t i = 0; i < size_alpha; i++) { - digamma_alpha[i] = digamma(alpha_vec.val(i)); - } - for (size_t i = 0; i < size_beta; i++) { - digamma_beta[i] = digamma(beta_vec.val(i)); - } - for (size_t i = 0; i < size_alpha_beta; i++) { - digamma_sum[i] = digamma(alpha_vec.val(i) + beta_vec.val(i)); - } - } - + constexpr bool calc_alpha = is_autodiff_v; + constexpr bool calc_beta = is_autodiff_v; for (size_t n = 0; n < N; n++) { - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - const T_partials_return betafunc_dbl = beta(alpha_dbl, beta_dbl); - const T_partials_return Pn = 1.0 - inc_beta(alpha_dbl, beta_dbl, y_dbl); - const T_partials_return inv_Pn + const auto y_dbl = y_vec.val(n); + const auto alpha_dbl = alpha_vec.val(n); + const auto beta_dbl = beta_vec.val(n); + const auto betafunc_dbl = beta(alpha_dbl, beta_dbl); + const auto Pn = 1.0 - inc_beta(alpha_dbl, beta_dbl, y_dbl); + const auto inv_Pn = is_constant_all::value ? 0 : inv(Pn); ccdf_log += log(Pn); @@ -111,18 +88,24 @@ inline return_type_t beta_lccdf( / betafunc_dbl; } - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - if constexpr (is_any_autodiff_v) { - grad_reg_inc_beta(g1, g2, alpha_dbl, beta_dbl, y_dbl, digamma_alpha[n], - digamma_beta[n], digamma_sum[n], betafunc_dbl); - } - if constexpr (is_autodiff_v) { + if constexpr (calc_alpha && calc_beta) { + const auto [g1, g2] = grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0); partials<1>(ops_partials)[n] -= g1 * inv_Pn; - } - if constexpr (is_autodiff_v) { partials<2>(ops_partials)[n] -= g2 * inv_Pn; + } else if constexpr (calc_alpha) { + partials<1>(ops_partials)[n] + -= grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, + digamma(alpha_dbl), 0.0, + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; + } else if constexpr (calc_beta) { + partials<2>(ops_partials)[n] + -= grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, 0.0, + digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; } } return ops_partials.build(ccdf_log); diff --git a/stan/math/prim/prob/beta_lcdf.hpp b/stan/math/prim/prob/beta_lcdf.hpp index c2a514a766c..22517788ed4 100644 --- a/stan/math/prim/prob/beta_lcdf.hpp +++ b/stan/math/prim/prob/beta_lcdf.hpp @@ -58,40 +58,18 @@ inline return_type_t beta_lcdf( scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); - const size_t size_alpha = stan::math::size(alpha); - const size_t size_beta = stan::math::size(beta_param); - const size_t size_alpha_beta = max_size(alpha, beta_param); const size_t N = max_size(y, alpha, beta_param); - // Allocate digamma buffers only if alpha/beta contain any autodiff scalars. - constexpr bool need_digamma = is_any_autodiff_v; - VectorBuilder digamma_alpha( - size_alpha); - VectorBuilder digamma_beta( - size_beta); - VectorBuilder - digamma_sum(size_alpha_beta); - - if constexpr (is_any_autodiff_v) { - for (size_t i = 0; i < size_alpha; ++i) { - digamma_alpha[i] = digamma(alpha_vec.val(i)); - } - for (size_t i = 0; i < size_beta; ++i) { - digamma_beta[i] = digamma(beta_vec.val(i)); - } - for (size_t i = 0; i < size_alpha_beta; ++i) { - digamma_sum[i] = digamma(alpha_vec.val(i) + beta_vec.val(i)); - } - } - + constexpr bool calc_alpha = is_autodiff_v; + constexpr bool calc_beta = is_autodiff_v; for (size_t n = 0; n < N; ++n) { - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); - const T_partials_return betafunc_dbl = beta(alpha_dbl, beta_dbl); - const T_partials_return Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl); + const auto y_dbl = y_vec.val(n); + const auto alpha_dbl = alpha_vec.val(n); + const auto beta_dbl = beta_vec.val(n); + const auto betafunc_dbl = beta(alpha_dbl, beta_dbl); + const auto Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl); - const T_partials_return inv_Pn + const auto inv_Pn = is_any_autodiff_v ? inv(Pn) : 0; cdf_log += log(Pn); @@ -102,17 +80,24 @@ inline return_type_t beta_lcdf( / betafunc_dbl; } - if constexpr (is_any_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - grad_reg_inc_beta(g1, g2, alpha_dbl, beta_dbl, y_dbl, digamma_alpha[n], - digamma_beta[n], digamma_sum[n], betafunc_dbl); - if constexpr (is_any_autodiff_v) { - partials<1>(ops_partials)[n] += g1 * inv_Pn; - } - if constexpr (is_any_autodiff_v) { - partials<2>(ops_partials)[n] += g2 * inv_Pn; - } + if constexpr (calc_alpha && calc_beta) { + const auto [g1, g2] = grad_reg_inc_beta( + alpha_dbl, beta_dbl, y_dbl, digamma(alpha_dbl), digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0); + partials<1>(ops_partials)[n] += g1 * inv_Pn; + partials<2>(ops_partials)[n] += g2 * inv_Pn; + } else if constexpr (calc_alpha) { + partials<1>(ops_partials)[n] + += grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, + digamma(alpha_dbl), 0.0, + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; + } else if constexpr (calc_beta) { + partials<2>(ops_partials)[n] + += grad_reg_inc_beta(alpha_dbl, beta_dbl, y_dbl, 0.0, + digamma(beta_dbl), + digamma(alpha_dbl + beta_dbl), 0.0) + * inv_Pn; } } diff --git a/stan/math/prim/prob/beta_proportion_lccdf.hpp b/stan/math/prim/prob/beta_proportion_lccdf.hpp index c615d23b0df..04e606fd3c8 100644 --- a/stan/math/prim/prob/beta_proportion_lccdf.hpp +++ b/stan/math/prim/prob/beta_proportion_lccdf.hpp @@ -71,32 +71,8 @@ inline return_type_t beta_proportion_lccdf( scalar_seq_view y_vec(y_ref); scalar_seq_view mu_vec(mu_ref); scalar_seq_view kappa_vec(kappa_ref); - size_t size_kappa = stan::math::size(kappa); - size_t size_mu_kappa = max_size(mu, kappa); size_t N = max_size(y, mu, kappa); - VectorBuilder, T_partials_return, T_loc, - T_prec> - digamma_mukappa(size_mu_kappa); - VectorBuilder, T_partials_return, T_loc, - T_prec> - digamma_kappa_mukappa(size_mu_kappa); - VectorBuilder, T_partials_return, T_prec> - digamma_kappa(size_kappa); - - if constexpr (is_any_autodiff_v) { - for (size_t i = 0; i < size_mu_kappa; i++) { - const T_partials_return kappa_dbl = kappa_vec.val(i); - const T_partials_return mukappa_dbl = mu_vec.val(i) * kappa_dbl; - digamma_mukappa[i] = digamma(mukappa_dbl); - digamma_kappa_mukappa[i] = digamma(kappa_dbl - mukappa_dbl); - } - - for (size_t i = 0; i < size_kappa; i++) { - digamma_kappa[i] = digamma(kappa_vec.val(i)); - } - } - for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); @@ -118,20 +94,18 @@ inline return_type_t beta_proportion_lccdf( / betafunc_dbl; } - T_partials_return g1 = 0; - T_partials_return g2 = 0; - + // the mu and kappa partials each need both shape partials if constexpr (is_any_autodiff_v) { - grad_reg_inc_beta(g1, g2, mukappa_dbl, kappa_mukappa_dbl, y_dbl, - digamma_mukappa[n], digamma_kappa_mukappa[n], - digamma_kappa[n], betafunc_dbl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] -= kappa_dbl * (g1 - g2) * inv_Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - -= (g1 * mu_dbl + g2 * (1 - mu_dbl)) * inv_Pn; + const auto [g1, g2] = grad_reg_inc_beta( + mukappa_dbl, kappa_mukappa_dbl, y_dbl, digamma(mukappa_dbl), + digamma(kappa_mukappa_dbl), digamma(kappa_dbl), 0.0); + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] -= kappa_dbl * (g1 - g2) * inv_Pn; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] + -= (g1 * mu_dbl + g2 * (1 - mu_dbl)) * inv_Pn; + } } } return ops_partials.build(ccdf_log); diff --git a/stan/math/prim/prob/beta_proportion_lcdf.hpp b/stan/math/prim/prob/beta_proportion_lcdf.hpp index 5566b8e5fe2..e3e11ee3534 100644 --- a/stan/math/prim/prob/beta_proportion_lcdf.hpp +++ b/stan/math/prim/prob/beta_proportion_lcdf.hpp @@ -72,31 +72,8 @@ inline return_type_t beta_proportion_lcdf( scalar_seq_view y_vec(y_ref); scalar_seq_view mu_vec(mu_ref); scalar_seq_view kappa_vec(kappa_ref); - size_t size_kappa = stan::math::size(kappa); - size_t size_mu_kappa = max_size(mu, kappa); size_t N = max_size(y, mu, kappa); - VectorBuilder, T_partials_return, T_loc, - T_prec> - digamma_mukappa(size_mu_kappa); - VectorBuilder, T_partials_return, T_loc, - T_prec> - digamma_kappa_mukappa(size_mu_kappa); - VectorBuilder, T_partials_return, T_prec> - digamma_kappa(size_kappa); - - if constexpr (is_any_autodiff_v) { - for (size_t i = 0; i < size_mu_kappa; i++) { - const T_partials_return kappa_dbl = kappa_vec.val(i); - const T_partials_return mukappa_dbl = mu_vec.val(i) * kappa_dbl; - digamma_mukappa[i] = digamma(mukappa_dbl); - digamma_kappa_mukappa[i] = digamma(kappa_dbl - mukappa_dbl); - } - for (size_t i = 0; i < size_kappa; i++) { - digamma_kappa[i] = digamma(kappa_vec.val(i)); - } - } - for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); @@ -118,20 +95,18 @@ inline return_type_t beta_proportion_lcdf( / betafunc_dbl; } - T_partials_return g1 = 0; - T_partials_return g2 = 0; - + // the mu and kappa partials each need both shape partials if constexpr (is_any_autodiff_v) { - grad_reg_inc_beta(g1, g2, mukappa_dbl, kappa_mukappa_dbl, y_dbl, - digamma_mukappa[n], digamma_kappa_mukappa[n], - digamma_kappa[n], betafunc_dbl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += kappa_dbl * (g1 - g2) * inv_Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += (g1 * mu_dbl + g2 * (1 - mu_dbl)) * inv_Pn; + const auto [g1, g2] = grad_reg_inc_beta( + mukappa_dbl, kappa_mukappa_dbl, y_dbl, digamma(mukappa_dbl), + digamma(kappa_mukappa_dbl), digamma(kappa_dbl), 0.0); + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] += kappa_dbl * (g1 - g2) * inv_Pn; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] + += (g1 * mu_dbl + g2 * (1 - mu_dbl)) * inv_Pn; + } } } diff --git a/stan/math/prim/prob/neg_binomial_2_cdf.hpp b/stan/math/prim/prob/neg_binomial_2_cdf.hpp index 47cab11d76d..6a96393e9cb 100644 --- a/stan/math/prim/prob/neg_binomial_2_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_cdf.hpp @@ -84,30 +84,28 @@ inline return_type_t neg_binomial_2_cdf( return ops_partials.build(1.0); } - const T_partials_return n_dbl_p1 = n_vec.val(i) + 1; - const T_partials_return mu_dbl = mu_vec.val(i); - const T_partials_return phi_dbl = phi_vec.val(i); - const T_partials_return inv_mu_plus_phi = inv(mu_dbl + phi_dbl); - const T_partials_return p_dbl = phi_dbl * inv_mu_plus_phi; - const T_partials_return d_dbl = square(inv_mu_plus_phi); - const T_partials_return P_i = inc_beta(phi_dbl, n_dbl_p1, p_dbl); - const T_partials_return inc_beta_ddz_i - = is_constant_all::value - ? 0 - : inc_beta_ddz(phi_dbl, n_dbl_p1, p_dbl) * d_dbl / P_i; + const auto n_dbl_p1 = n_vec.val(i) + 1; + const auto mu_dbl = mu_vec.val(i); + const auto phi_dbl = phi_vec.val(i); + const auto inv_mu_plus_phi = inv(mu_dbl + phi_dbl); + const auto p_dbl = phi_dbl * inv_mu_plus_phi; + const auto P_i = inc_beta(phi_dbl, n_dbl_p1, p_dbl); P *= P_i; - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] -= inc_beta_ddz_i * phi_dbl; - } - - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] - += inc_beta_dda(phi_dbl, n_dbl_p1, p_dbl, digamma_phi_vec[i], - digamma_sum_vec[i]) - / P_i - + inc_beta_ddz_i * mu_dbl; + if constexpr (is_any_autodiff_v) { + const auto d_dbl = square(inv_mu_plus_phi); + const auto inc_beta_ddz_i + = inc_beta_ddz(phi_dbl, n_dbl_p1, p_dbl) * d_dbl / P_i; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] -= inc_beta_ddz_i * phi_dbl; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[i] + += inc_beta_dda(phi_dbl, n_dbl_p1, p_dbl, digamma_phi_vec[i], + digamma_sum_vec[i]) + / P_i + + inc_beta_ddz_i * mu_dbl; + } } } diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp index 31163102be1..b4e7c2b436a 100644 --- a/stan/math/prim/prob/neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_cdf.hpp @@ -3,122 +3,115 @@ #include #include +#include +#include +#include #include #include #include #include #include -#include -#include -#include +#include +#include #include #include -#include +#include #include -#include #include namespace stan { namespace math { -template +/** \ingroup prob_dists + * Returns the CDF of the negative binomial distribution. Given containers of + * matching sizes, returns the product of probabilities. + * + * With p = beta / (1 + beta), the CDF at n is the regularized incomplete beta + * function I_p(alpha, n + 1). + * + * @tparam T_n type of integer parameter + * @tparam T_shape type of shape parameter + * @tparam T_inv_scale type of inverse scale parameter + * @param n failures variable + * @param alpha shape parameter + * @param beta inverse scale parameter + * @return probability or product of probabilities + * @throw std::domain_error if alpha or beta is not positive and finite + * @throw std::invalid_argument if container sizes mismatch + */ +template * = nullptr> inline return_type_t neg_binomial_cdf( - const T_n& n, const T_shape& alpha, const T_inv_scale& beta) { + T_n&& n, T_shape&& alpha, T_inv_scale&& beta) { using T_partials_return = partials_return_t; using T_n_ref = ref_type_t; - using T_alpha_ref = ref_type_t; - using T_beta_ref = ref_type_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 = "neg_binomial_cdf"; check_consistent_sizes(function, "Failures variable", n, "Shape parameter", alpha, "Inverse scale parameter", beta); - T_n_ref n_ref = n; - T_alpha_ref alpha_ref = alpha; - T_beta_ref beta_ref = beta; - check_positive_finite(function, "Shape parameter", alpha_ref); - check_positive_finite(function, "Inverse scale parameter", beta_ref); + T_n_ref n_ref = std::forward(n); + T_alpha_ref alpha_ref = std::forward(alpha); + T_beta_ref beta_ref = std::forward(beta); - if (size_zero(n, alpha, beta)) { + decltype(auto) n_val = to_ref(as_value_column_array_or_scalar(n_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_positive_finite(function, "Shape parameter", alpha_val); + check_positive_finite(function, "Inverse scale parameter", beta_val); + + if (size_zero(n_ref, alpha_ref, beta_ref)) { return 1.0; } - T_partials_return P(1.0); auto ops_partials = make_partials_propagator(alpha_ref, beta_ref); - scalar_seq_view n_vec(n_ref); - scalar_seq_view alpha_vec(alpha_ref); - scalar_seq_view beta_vec(beta_ref); - size_t size_alpha = stan::math::size(alpha); - size_t size_n_alpha = max_size(n, alpha); - size_t max_size_seq_view = max_size(n, alpha, beta); - // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(0.0); - } - } - - VectorBuilder, T_partials_return, T_shape> - digamma_alpha_vec(size_alpha); - VectorBuilder, T_partials_return, T_n, T_shape> - digamma_sum_vec(size_n_alpha); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < size_alpha; i++) { - digamma_alpha_vec[i] = digamma(alpha_vec.val(i)); - } - for (size_t i = 0; i < size_n_alpha; i++) { - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - digamma_sum_vec[i] = digamma(n_dbl + alpha_dbl + 1); - } - } - - for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_vec.val(i) == std::numeric_limits::max()) { - return ops_partials.build(1.0); - } - - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - const T_partials_return beta_dbl = beta_vec.val(i); - const T_partials_return inv_beta_p1 = inv(beta_dbl + 1); - const T_partials_return p_dbl = beta_dbl * inv_beta_p1; - const T_partials_return d_dbl = square(inv_beta_p1); - - const T_partials_return P_i = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); - - P *= P_i; - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[i] - += inc_beta_dda(alpha_dbl, n_dbl + 1, p_dbl, digamma_alpha_vec[i], - digamma_sum_vec[i]) - / P_i; - } - - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[i] - += inc_beta_ddz(alpha_dbl, n_dbl + 1.0, p_dbl) * d_dbl / P_i; - } + if (any(n_val < 0)) { + return ops_partials.build(0.0); } - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < size_alpha; ++i) { - partials<0>(ops_partials)[i] *= P; - } + if (any(n_val == std::numeric_limits::max())) { + return ops_partials.build(1.0); } - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < stan::math::size(beta); ++i) { - partials<1>(ops_partials)[i] *= P; + if constexpr (is_autodiff_v || is_autodiff_v) { + auto inv_beta_p1 = to_ref(inv(beta_val + 1.0)); + auto p_val = to_ref(beta_val * inv_beta_p1); + auto n_p1 = to_ref(promote_scalar(n_val + 1.0)); + auto P_i = to_ref(as_array_or_scalar(inc_beta(alpha_val, n_p1, p_val))); + const T_partials_return P = prod(P_i); + // d/dx prod(P_i) = sum(dP_i/dx * P / P_i) + if constexpr (is_autodiff_v && is_autodiff_v) { + auto P_div_P_i = to_ref(P / P_i); + partials<0>(ops_partials) + = inc_beta_dda(alpha_val, n_p1, p_val, to_ref(digamma(alpha_val)), + to_ref(digamma(alpha_val + n_p1))) + * P_div_P_i; + partials<1>(ops_partials) + = as_array_or_scalar(inc_beta_ddz(alpha_val, n_p1, p_val)) + * square(inv_beta_p1) * P_div_P_i; + } else if constexpr (is_autodiff_v) { + partials<0>(ops_partials) + = inc_beta_dda(alpha_val, n_p1, p_val, to_ref(digamma(alpha_val)), + to_ref(digamma(alpha_val + n_p1))) + / P_i * P; + } else if constexpr (is_autodiff_v) { + partials<1>(ops_partials) + = as_array_or_scalar(inc_beta_ddz(alpha_val, n_p1, p_val)) + * square(inv_beta_p1) / P_i * P; } + return ops_partials.build(P); + } else { + auto inv_beta_p1 = to_ref(inv(beta_val + 1.0)); + auto p_val = to_ref(beta_val * inv_beta_p1); + auto n_p1 = to_ref(promote_scalar(n_val + 1.0)); + const T_partials_return P + = prod(as_array_or_scalar(inc_beta(alpha_val, n_p1, p_val))); + return ops_partials.build(P); } - - return ops_partials.build(P); } } // namespace math diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp index cc24a6b79a7..3de84a34589 100644 --- a/stan/math/prim/prob/neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp @@ -53,8 +53,6 @@ inline return_type_t neg_binomial_lccdf( scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); size_t size_n = stan::math::size(n); - size_t size_alpha = stan::math::size(alpha); - size_t size_n_alpha = max_size(n, alpha); size_t max_size_seq_view = max_size(n, alpha, beta_param); // Explicit return for extreme values @@ -65,27 +63,6 @@ inline return_type_t neg_binomial_lccdf( } } - VectorBuilder, T_partials_return, T_n> digammaN_vec( - size_n); - VectorBuilder, T_partials_return, T_shape> - digammaAlpha_vec(size_alpha); - VectorBuilder, T_partials_return, T_n, T_shape> - digammaSum_vec(size_n_alpha); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < size_n; i++) { - digammaN_vec[i] = digamma(n_vec.val(i) + 1); - } - for (size_t i = 0; i < size_alpha; i++) { - digammaAlpha_vec[i] = digamma(alpha_vec.val(i)); - } - for (size_t i = 0; i < size_n_alpha; i++) { - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1); - } - } - for (size_t i = 0; i < max_size_seq_view; i++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -93,25 +70,24 @@ inline return_type_t neg_binomial_lccdf( return ops_partials.build(negative_infinity()); } - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - const T_partials_return beta_dbl = beta_vec.val(i); - const T_partials_return inv_beta_p1 = inv(beta_dbl + 1); - const T_partials_return p_dbl = beta_dbl * inv_beta_p1; - const T_partials_return d_dbl = square(inv_beta_p1); - const T_partials_return Pi = 1.0 - inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); - const T_partials_return beta_func = beta(n_dbl + 1, alpha_dbl); + const auto n_dbl = n_vec.val(i); + const auto alpha_dbl = alpha_vec.val(i); + const auto beta_dbl = beta_vec.val(i); + const auto inv_beta_p1 = inv(beta_dbl + 1); + const auto p_dbl = beta_dbl * inv_beta_p1; + const auto d_dbl = square(inv_beta_p1); + const auto Pi = 1.0 - inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); + const auto beta_func = beta(n_dbl + 1, alpha_dbl); P += log(Pi); + // alpha is the first shape of inc_beta, so only its partial is needed if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, alpha_dbl, n_dbl + 1, p_dbl, - digammaAlpha_vec[i], digammaN_vec[i], digammaSum_vec[i], - beta_func); - partials<0>(ops_partials)[i] -= g1 / Pi; + partials<0>(ops_partials)[i] + -= grad_reg_inc_beta( + alpha_dbl, n_dbl + 1, p_dbl, digamma(alpha_dbl), 0.0, + digamma(n_dbl + alpha_dbl + 1), 0.0) + / Pi; } if constexpr (is_autodiff_v) { partials<1>(ops_partials)[i] -= d_dbl * pow(1 - p_dbl, n_dbl) diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp index ba85ba41d0d..7296870e58e 100644 --- a/stan/math/prim/prob/neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp @@ -53,8 +53,6 @@ inline return_type_t neg_binomial_lcdf( scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); size_t size_n = stan::math::size(n); - size_t size_alpha = stan::math::size(alpha); - size_t size_n_alpha = max_size(n, alpha); size_t max_size_seq_view = max_size(n, alpha, beta_param); // Explicit return for extreme values @@ -65,27 +63,6 @@ inline return_type_t neg_binomial_lcdf( } } - VectorBuilder, T_partials_return, T_n> digammaN_vec( - size_n); - VectorBuilder, T_partials_return, T_shape> - digammaAlpha_vec(size_alpha); - VectorBuilder, T_partials_return, T_n, T_shape> - digammaSum_vec(size_n_alpha); - - if constexpr (is_autodiff_v) { - for (size_t i = 0; i < size_n; i++) { - digammaN_vec[i] = digamma(n_vec.val(i) + 1); - } - for (size_t i = 0; i < size_alpha; i++) { - digammaAlpha_vec[i] = digamma(alpha_vec.val(i)); - } - for (size_t i = 0; i < size_n_alpha; i++) { - const T_partials_return n_dbl = n_vec.val(i); - const T_partials_return alpha_dbl = alpha_vec.val(i); - digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1); - } - } - for (size_t i = 0; i < max_size_seq_view; i++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -104,14 +81,13 @@ inline return_type_t neg_binomial_lcdf( P += log(Pi); + // alpha is the first shape of inc_beta, so only its partial is needed if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, alpha_dbl, n_dbl + 1, p_dbl, - digammaAlpha_vec[i], digammaN_vec[i], digammaSum_vec[i], - beta_func); - partials<0>(ops_partials)[i] += g1 / Pi; + partials<0>(ops_partials)[i] + += grad_reg_inc_beta( + alpha_dbl, n_dbl + 1, p_dbl, digamma(alpha_dbl), 0.0, + digamma(n_dbl + alpha_dbl + 1), 0.0) + / Pi; } if constexpr (is_autodiff_v) { partials<1>(ops_partials)[i] += d_dbl * pow(1 - p_dbl, n_dbl) diff --git a/stan/math/prim/prob/student_t_cdf.hpp b/stan/math/prim/prob/student_t_cdf.hpp index f9720c61209..78cfa8f5ee3 100644 --- a/stan/math/prim/prob/student_t_cdf.hpp +++ b/stan/math/prim/prob/student_t_cdf.hpp @@ -60,26 +60,6 @@ inline return_type_t student_t_cdf( } } - T_partials_return digammaHalf = 0; - - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digammaNu_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> - digammaNuPlusHalf_vec(math::size(nu)); - - if constexpr (is_autodiff_v) { - digammaHalf = digamma(0.5); - - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - - digammaNu_vec[i] = digamma(0.5 * nu_dbl); - digammaNuPlusHalf_vec[i] = digamma(0.5 + 0.5 * nu_dbl); - } - } - for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -87,20 +67,19 @@ inline return_type_t student_t_cdf( continue; } - const T_partials_return sigma_inv = 1.0 / sigma_vec.val(n); - const T_partials_return t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; - const T_partials_return nu_dbl = nu_vec.val(n); - const T_partials_return q = nu_dbl / (t * t); - const T_partials_return r = 1.0 / (1.0 + q); - const T_partials_return J = 2 * r * r * q / t; - const T_partials_return betaNuHalf = beta(0.5, 0.5 * nu_dbl); + const auto sigma_inv = 1.0 / sigma_vec.val(n); + const auto t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; + const auto nu_dbl = nu_vec.val(n); + const auto q = nu_dbl / (t * t); + const auto r = 1.0 / (1.0 + q); + const auto J = 2 * r * r * q / t; + const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); double zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - T_partials_return z - = inc_beta(0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r); - const T_partials_return Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - const T_partials_return d_ibeta + auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; + const auto d_ibeta = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; P *= Pn; @@ -110,12 +89,10 @@ inline return_type_t student_t_cdf( += -zJacobian * d_ibeta * J * sigma_inv / Pn; } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r, - digammaNu_vec[n], digammaHalf, - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the first shape of inc_beta here and 1 / 2 the second + const auto g1 = grad_reg_inc_beta( + 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] += zJacobian * (d_ibeta * (r / t) * (r / t) + 0.5 * g1) / Pn; @@ -131,15 +108,13 @@ inline return_type_t student_t_cdf( } } else { - T_partials_return z - = 1.0 - inc_beta((T_partials_return)0.5, 0.5 * nu_dbl, r); + auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); zJacobian *= -1; - const T_partials_return Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; + const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - T_partials_return d_ibeta - = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P *= Pn; @@ -148,12 +123,10 @@ inline return_type_t student_t_cdf( += zJacobian * d_ibeta * J * sigma_inv / Pn; } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, (T_partials_return)0.5, 0.5 * nu_dbl, r, - digammaHalf, digammaNu_vec[n], - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the second shape of inc_beta here and 1 / 2 the first + const auto g2 = grad_reg_inc_beta( + 0.5, 0.5 * nu_dbl, r, 0.0, digamma(0.5 * nu_dbl), + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] += zJacobian * (-d_ibeta * (r / t) * (r / t) + 0.5 * g2) / Pn; diff --git a/stan/math/prim/prob/student_t_lccdf.hpp b/stan/math/prim/prob/student_t_lccdf.hpp index af8cb6208b0..224a3165f50 100644 --- a/stan/math/prim/prob/student_t_lccdf.hpp +++ b/stan/math/prim/prob/student_t_lccdf.hpp @@ -62,26 +62,6 @@ inline return_type_t student_t_lccdf( } } - T_partials_return digammaHalf = 0; - - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digammaNu_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> - digammaNuPlusHalf_vec(math::size(nu)); - - if constexpr (is_autodiff_v) { - digammaHalf = digamma(0.5); - - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - - digammaNu_vec[i] = digamma(0.5 * nu_dbl); - digammaNuPlusHalf_vec[i] = digamma(0.5 + 0.5 * nu_dbl); - } - } - for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -89,20 +69,19 @@ inline return_type_t student_t_lccdf( return ops_partials.build(negative_infinity()); } - const T_partials_return sigma_inv = 1.0 / sigma_vec.val(n); - const T_partials_return t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; - const T_partials_return nu_dbl = nu_vec.val(n); - const T_partials_return q = nu_dbl / (t * t); - const T_partials_return r = 1.0 / (1.0 + q); - const T_partials_return J = 2 * r * r * q / t; - const T_partials_return betaNuHalf = beta(0.5, 0.5 * nu_dbl); - T_partials_return zJacobian = t > 0 ? -0.5 : 0.5; + const auto sigma_inv = 1.0 / sigma_vec.val(n); + const auto t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; + const auto nu_dbl = nu_vec.val(n); + const auto q = nu_dbl / (t * t); + const auto r = 1.0 / (1.0 + q); + const auto J = 2 * r * r * q / t; + const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); + auto zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - T_partials_return z - = inc_beta(0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r); - const T_partials_return Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; - const T_partials_return d_ibeta + auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + const auto Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; + const auto d_ibeta = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; P += log(Pn); @@ -113,12 +92,10 @@ inline return_type_t student_t_lccdf( } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r, - digammaNu_vec[n], digammaHalf, - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the first shape of inc_beta here and 1 / 2 the second + const auto g1 = grad_reg_inc_beta( + 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] -= zJacobian * (d_ibeta * (r / t) * (r / t) + 0.5 * g1) / Pn; @@ -134,14 +111,12 @@ inline return_type_t student_t_lccdf( } } else { - T_partials_return z - = 1.0 - inc_beta((T_partials_return)0.5, 0.5 * nu_dbl, r); + auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); zJacobian *= -1; - const T_partials_return Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; + const auto Pn = t > 0 ? 0.5 * z : 1.0 - 0.5 * z; - T_partials_return d_ibeta - = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P += log(Pn); @@ -151,12 +126,10 @@ inline return_type_t student_t_lccdf( } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, (T_partials_return)0.5, 0.5 * nu_dbl, r, - digammaHalf, digammaNu_vec[n], - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the second shape of inc_beta here and 1 / 2 the first + const auto g2 = grad_reg_inc_beta( + 0.5, 0.5 * nu_dbl, r, 0.0, digamma(0.5 * nu_dbl), + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] -= zJacobian * (-d_ibeta * (r / t) * (r / t) + 0.5 * g2) / Pn; diff --git a/stan/math/prim/prob/student_t_lcdf.hpp b/stan/math/prim/prob/student_t_lcdf.hpp index 793b23273ed..be791924977 100644 --- a/stan/math/prim/prob/student_t_lcdf.hpp +++ b/stan/math/prim/prob/student_t_lcdf.hpp @@ -62,26 +62,6 @@ inline return_type_t student_t_lcdf( } } - T_partials_return digammaHalf = 0; - - VectorBuilder, T_partials_return, T_dof> digamma_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> digammaNu_vec( - math::size(nu)); - VectorBuilder, T_partials_return, T_dof> - digammaNuPlusHalf_vec(math::size(nu)); - - if constexpr (is_autodiff_v) { - digammaHalf = digamma(0.5); - - for (size_t i = 0; i < stan::math::size(nu); i++) { - const T_partials_return nu_dbl = nu_vec.val(i); - - digammaNu_vec[i] = digamma(0.5 * nu_dbl); - digammaNuPlusHalf_vec[i] = digamma(0.5 + 0.5 * nu_dbl); - } - } - for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero @@ -89,20 +69,19 @@ inline return_type_t student_t_lcdf( continue; } - const T_partials_return sigma_inv = 1.0 / sigma_vec.val(n); - const T_partials_return t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; - const T_partials_return nu_dbl = nu_vec.val(n); - const T_partials_return q = nu_dbl / (t * t); - const T_partials_return r = 1.0 / (1.0 + q); - const T_partials_return J = 2 * r * r * q / t; - const T_partials_return betaNuHalf = beta(0.5, 0.5 * nu_dbl); - T_partials_return zJacobian = t > 0 ? -0.5 : 0.5; + const auto sigma_inv = 1.0 / sigma_vec.val(n); + const auto t = (y_vec.val(n) - mu_vec.val(n)) * sigma_inv; + const auto nu_dbl = nu_vec.val(n); + const auto q = nu_dbl / (t * t); + const auto r = 1.0 / (1.0 + q); + const auto J = 2 * r * r * q / t; + const auto betaNuHalf = beta(0.5, 0.5 * nu_dbl); + auto zJacobian = t > 0 ? -0.5 : 0.5; if (q < 2) { - T_partials_return z - = inc_beta(0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r); - const T_partials_return Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - const T_partials_return d_ibeta + auto z = inc_beta(0.5 * nu_dbl, 0.5, 1.0 - r); + const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; + const auto d_ibeta = pow(r, -0.5) * pow(1.0 - r, 0.5 * nu_dbl - 1) / betaNuHalf; P += log(Pn); @@ -113,12 +92,10 @@ inline return_type_t student_t_lcdf( } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, 0.5 * nu_dbl, (T_partials_return)0.5, 1.0 - r, - digammaNu_vec[n], digammaHalf, - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the first shape of inc_beta here and 1 / 2 the second + const auto g1 = grad_reg_inc_beta( + 0.5 * nu_dbl, 0.5, 1.0 - r, digamma(0.5 * nu_dbl), 0.0, + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] += zJacobian * (d_ibeta * (r / t) * (r / t) + 0.5 * g1) / Pn; @@ -134,14 +111,12 @@ inline return_type_t student_t_lcdf( } } else { - T_partials_return z - = 1.0 - inc_beta((T_partials_return)0.5, 0.5 * nu_dbl, r); + auto z = 1.0 - inc_beta(0.5, 0.5 * nu_dbl, r); zJacobian *= -1; - const T_partials_return Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; + const auto Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z; - T_partials_return d_ibeta - = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; + auto d_ibeta = pow(1.0 - r, 0.5 * nu_dbl - 1) * pow(r, -0.5) / betaNuHalf; P += log(Pn); @@ -151,12 +126,10 @@ inline return_type_t student_t_lcdf( } if constexpr (is_autodiff_v) { - T_partials_return g1 = 0; - T_partials_return g2 = 0; - - grad_reg_inc_beta(g1, g2, (T_partials_return)0.5, 0.5 * nu_dbl, r, - digammaHalf, digammaNu_vec[n], - digammaNuPlusHalf_vec[n], betaNuHalf); + // nu / 2 is the second shape of inc_beta here and 1 / 2 the first + const auto g2 = grad_reg_inc_beta( + 0.5, 0.5 * nu_dbl, r, 0.0, digamma(0.5 * nu_dbl), + digamma(0.5 + 0.5 * nu_dbl), 0.0); partials<1>(ops_partials)[n] += zJacobian * (-d_ibeta * (r / t) * (r / t) + 0.5 * g2) / Pn; diff --git a/test/unit/math/rev/prob/beta_cdf_log_test.cpp b/test/unit/math/rev/prob/beta_cdf_log_test.cpp new file mode 100644 index 00000000000..ec46553b456 --- /dev/null +++ b/test/unit/math/rev/prob/beta_cdf_log_test.cpp @@ -0,0 +1,46 @@ +#include +#include +#include +#include +#include + +namespace { + +/** + * Shape gradients of beta_lcdf against a high-precision reference. The + * reference is (d/da I) / I and (d/db I) / I from mpmath at 80 digits. + */ +void expect_beta_lcdf_grads(double y, double alpha, double beta, double dalpha, + double dbeta, double rtol) { + using stan::math::var; + var a = alpha; + var b = beta; + var lp = stan::math::beta_lcdf(y, a, b); + std::vector x{a, b}; + std::vector g; + lp.grad(x, g); + EXPECT_TRUE(std::isfinite(g[0])) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_TRUE(std::isfinite(g[1])) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_NEAR(dalpha, g[0], rtol * std::fabs(dalpha)) + << "y=" << y << " a=" << alpha << " b=" << beta; + EXPECT_NEAR(dbeta, g[1], rtol * std::fabs(dbeta)) + << "y=" << y << " a=" << alpha << " b=" << beta; +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsBeta_lcdf_shape_gradients) { + expect_beta_lcdf_grads(0.25, 2, 3, -0.87136898798556226, 0.40853599542368053, + 1e-12); + // large shapes: beta(alpha, beta) underflows + expect_beta_lcdf_grads(0.5, 600, 600, -0.0230409429460951, 0.0230409429460951, + 1e-12); + // y far above the mean: the cdf is 1 and the gradients are 1e-138 + expect_beta_lcdf_grads(0.758, 352, 590, -3.3772717797123241e-138, + 4.5390666074606647e-138, 1e-12); + // integer beta + expect_beta_lcdf_grads(0.5, 10, 1, -0.69314718055994531, 2.3201212683796071, + 1e-12); +} diff --git a/test/unit/math/rev/prob/beta_cdf_test.cpp b/test/unit/math/rev/prob/beta_cdf_test.cpp new file mode 100644 index 00000000000..3f842144099 --- /dev/null +++ b/test/unit/math/rev/prob/beta_cdf_test.cpp @@ -0,0 +1,39 @@ +#include +#include +#include +#include +#include + +/** + * For one observation the shape partials of beta_cdf are d/da I and d/db I. + * Reference: mpmath, 80 digits. The point has b >> a with a small y. + */ +TEST_F(AgradRev, ProbDistributionsBeta_cdf_shape_gradients) { + using stan::math::var; + var a = 1.145; + var b = 6786; + var p = stan::math::beta_cdf(1.8e-4, a, b); + std::vector x{a, b}; + std::vector g; + p.grad(x, g); + EXPECT_NEAR(0.64944133941212122, p.val(), 1e-12); + EXPECT_NEAR(-0.38621356822721375, g[0], 1e-10); + EXPECT_NEAR(5.8439109664263782e-5, g[1], 1e-10 * 5.8e-5); +} + +TEST_F(AgradRev, ProbDistributionsBeta_cdf_mixed_shape_types) { + // one autodiff shape with the other constant must give the same partial + using stan::math::var; + { + var a = 1.145; + var p = stan::math::beta_cdf(1.8e-4, a, 6786.0); + p.grad(); + EXPECT_NEAR(-0.38621356822721375, a.adj(), 1e-10); + } + { + var b = 6786; + var p = stan::math::beta_cdf(1.8e-4, 1.145, b); + p.grad(); + EXPECT_NEAR(5.8439109664263782e-5, b.adj(), 1e-10 * 5.8e-5); + } +}