diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index df349a7b483..bbab27af27c 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -313,19 +313,21 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit, opencl_kernels::inv_logit_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::log1m_device_function, opencl_kernels::logit_device_function) -ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function) +ADD_UNARY_FUNCTION_WITH_INCLUDES( + Phi, opencl_kernels::std_normal_lcdf_device_function, + opencl_kernels::phi_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx, opencl_kernels::inv_logit_device_function, opencl_kernels::phi_approx_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES( - std_normal_lcdf_scaled_impl, - opencl_kernels::std_normal_lcdf_device_function) + std_normal_lcdf_impl, opencl_kernels::std_normal_lcdf_device_function) +ADD_UNARY_FUNCTION_WITH_INCLUDES( + std_normal_lcdf_derivative, opencl_kernels::std_normal_lcdf_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES( - std_normal_lcdf_dscaled_impl, - opencl_kernels::std_normal_lcdf_device_function) -ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_Phi, opencl_kernels::log1m_device_function, - opencl_kernels::phi_device_function, - opencl_kernels::inv_phi_device_function) + inv_Phi, opencl_kernels::log1m_device_function, + opencl_kernels::std_normal_lcdf_device_function, + opencl_kernels::phi_device_function, + opencl_kernels::inv_phi_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES( log1m_inv_logit, opencl_kernels::log1p_exp_device_function, opencl_kernels::log1m_inv_logit_device_function) diff --git a/stan/math/opencl/kernels/device_functions/Phi.hpp b/stan/math/opencl/kernels/device_functions/Phi.hpp index 0db9aa2167b..bb9ff7eb331 100644 --- a/stan/math/opencl/kernels/device_functions/Phi.hpp +++ b/stan/math/opencl/kernels/device_functions/Phi.hpp @@ -20,17 +20,7 @@ static constexpr const char* phi_device_function * * @return Phi(x) */ - inline double Phi(double x) { - if (x < -37.5) { - return 0; - } else if (x < -5.0) { - return 0.5 * erfc(-M_SQRT1_2 * x); - } else if (x > 8.25) { - return 1; - } else { - return 0.5 * (1.0 + erf(M_SQRT1_2 * x)); - } - } + inline double Phi(double x) { return exp(std_normal_lcdf_impl(x)); } // \cond ) "\n#endif\n"; // NOLINT // \endcond diff --git a/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp b/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp index c48ce98acea..7180a7e504c 100644 --- a/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp +++ b/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp @@ -14,151 +14,104 @@ static constexpr const char* std_normal_lcdf_device_function "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_STD_NORMAL_LCDF\n" "#define " "STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_STD_NORMAL_LCDF\n" STRINGIFY( - /** \ingroup opencl_kernels - * Return the log standard normal cumulative distribution function - * evaluated from the scaled input `x / sqrt(2)`. - * - * @param scaled_y input scaled by `1 / sqrt(2)` - * @return log(Phi(x)) - */ - inline double std_normal_lcdf_scaled_impl(double scaled_y) { - double lcdf_n; - if (scaled_y > 0.0) { - // CDF(x) = 1/2 + 1/2 erf(x) = 1 - 1/2 erfc(x) - lcdf_n = log1p(-0.5 * erfc(scaled_y)); - if (isnan(lcdf_n)) { - lcdf_n = 0; - } - } else if (scaled_y > -4.0) { - // CDF(x) = 1/2 - 1/2 erf(-x) = 1/2 erfc(-x) - lcdf_n = log(erfc(-scaled_y)) - M_LN2; - } else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) { - // Need direct approximation once erfc(-x) underflows. - const double x2 = scaled_y * scaled_y; - const double x4 = pow(scaled_y, 4); - const double x6 = pow(scaled_y, 6); - const double x8 = pow(scaled_y, 8); - const double x10 = pow(scaled_y, 10); - const double temp_p = 0.000658749161529837803157 - + 0.0160837851487422766278 / x2 - + 0.125781726111229246204 / x4 - + 0.360344899949804439429 / x6 - + 0.305326634961232344035 / x8 - + 0.0163153871373020978498 / x10; - const double temp_q - = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 - - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6 - - 2.56852019228982242072 / x8 - 1.0 / x10; - lcdf_n = log(0.5 * M_2_SQRTPI + (temp_p / temp_q) / x2) - M_LN2 - - log(-scaled_y) - x2; - } else { - lcdf_n = -INFINITY; + // Cody (1969) rationals, sharing the CPU kernel's coefficients. + inline double std_normal_erf_small(double x) { + const double a[] = {3.16112374387056560, 1.13864154151050156e2, + 3.77485237685302021e2, 3.20937758913846947e3, + 1.85777706184603153e-1}; + const double b[] = {2.36012909523441209e1, 2.44024637934444173e2, + 1.28261652607737228e3, 2.84423683343917062e3}; + const double x2 = x * x; + double numerator = a[4] * x2; + double denominator = x2; + for (int i = 0; i < 3; ++i) { + numerator = (numerator + a[i]) * x2; + denominator = (denominator + b[i]) * x2; } - return lcdf_n; + return x * (numerator + a[3]) / (denominator + b[3]); } - /** \ingroup opencl_kernels - * Return the derivative of log standard normal cumulative - * distribution function with respect to the scaled input - * `x / sqrt(2)`. - * - * @param scaled_y input scaled by `1 / sqrt(2)` - * @return d / d(scaled_y) log(Phi(x)) - */ - inline double std_normal_lcdf_dscaled_impl(double scaled_y) { - double dnlcdf = 0.0; - double t = 0.0; - double t2 = 0.0; - double t4 = 0.0; - const double x2 = scaled_y * scaled_y; + inline double std_normal_lcdf_tail_correction(double r) { + const double p[] + = {0.000658749161529837803157, 0.0160837851487422766278, + 0.125781726111229246204, 0.360344899949804439429, + 0.305326634961232344035, 0.0163153871373020978498}; + const double q[] + = {-0.00233520497626869185443, -0.0605183413124413191178, + -0.527905102951428412248, -1.87295284992346047209, + -2.56852019228982242072, -1.0}; + double numerator = p[5] * r + p[4]; + double denominator = q[5] * r + q[4]; + for (int i = 3; i >= 0; --i) { + numerator = numerator * r + p[i]; + denominator = denominator * r + q[i]; + } + return (numerator / denominator) / (0.5 * M_2_SQRTPI); + } + + /** erfcx(x) = exp(x^2) erfc(x) for x >= 0.46875. */ + inline double std_normal_erfcx(double x) { + if (x > 4.0) { + const double r = 1.0 / (x * x); + return (1.0 + r * std_normal_lcdf_tail_correction(r)) + * (0.5 * M_2_SQRTPI) / x; + } + const double c[] = {5.64188496988670089e-1, 8.88314979438837594, + 6.61191906371416295e1, 2.98635138197400131e2, + 8.81952221241769090e2, 1.71204761263407058e3, + 2.05107837782607147e3, 1.23033935479799725e3, + 2.15311535474403846e-8}; + const double d[] = {1.57449261107098347e1, 1.17693950891312499e2, + 5.37181101862009858e2, 1.62138957456669019e3, + 3.29079923573345963e3, 4.36261909014324716e3, + 3.43936767414372164e3, 1.23033935480374942e3}; + double numerator = c[8] * x; + double denominator = x; + for (int i = 0; i < 7; ++i) { + numerator = (numerator + c[i]) * x; + denominator = (denominator + d[i]) * x; + } + return (numerator + c[7]) / (denominator + d[7]); + } - if (scaled_y > 2.9) { - t = 1.0 / (1.0 + 0.3275911 * scaled_y); - t2 = t * t; - t4 = pow(t, 4); - // A&S 7.1.26 keeps exp(-x2) in the numerator, as R's pnorm - // does; refs in stan/math/prim/prob/std_normal_lcdf.hpp - const double exp_m_x2 = exp(-x2); - dnlcdf - = 0.5 * M_2_SQRTPI * exp_m_x2 - / (1.0 - - exp_m_x2 - * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 - - 1.453152027 * t2 * t + 1.061405429 * t4)); - } else if (scaled_y > 2.5) { - t = scaled_y - 2.7; - t2 = t * t; - t4 = pow(t, 4); - dnlcdf = 0.0003849882382 - 0.002079084702 * t - + 0.005229340880 * t2 - 0.008029540137 * t2 * t - + 0.008232190507 * t4 - 0.005692364250 * t4 * t - + 0.002399496363 * pow(t, 6); - } else if (scaled_y > 2.1) { - t = scaled_y - 2.3; - t2 = t * t; - t4 = pow(t, 4); - dnlcdf = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 - - 0.03326906904 * t2 * t + 0.02482478940 * t4 - - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); - } else if (scaled_y > 1.5) { - t = scaled_y - 1.85; - t2 = t * t; - t4 = pow(t, 4); - dnlcdf = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 - - 0.09274533184 * t2 * t + 0.03543327418 * t4 - + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); - } else if (scaled_y > 0.8) { - t = scaled_y - 1.15; - t2 = t * t; - t4 = pow(t, 4); - dnlcdf = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 - - 0.09748053605 * t2 * t - 0.04347986191 * t4 - + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); - } else if (scaled_y > 0.1) { - t = scaled_y - 0.45; - t2 = t * t; - t4 = pow(t, 4); - dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 - + 0.04700850676 * t2 * t - 0.03478651979 * t4 - - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); - } else if (scaled_y < -29.0) { - // asymptotic Mills ratio, DLMF 7.12.1; grows linearly as - // -2*scaled_y, same 1/x^2 series shape as R's pnorm uses - const double inv_x2 = 1.0 / x2; - dnlcdf - = -2.0 * scaled_y - / (1.0 - + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875))); - } else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) { - t = 1.0 / (1.0 - 0.3275911 * scaled_y); - t2 = t * t; - t4 = pow(t, 4); - dnlcdf - = M_2_SQRTPI - / (0.254829592 * t - 0.284496736 * t2 + 1.421413741 * t2 * t - - 1.453152027 * t4 + 1.061405429 * t4 * t); - if (scaled_y < -17.0) { - dnlcdf += 0.0001263257217272 * x2 * scaled_y - + 0.0123586859488623 * x2 - - 0.0860505264736028 * scaled_y - 1.252783383752970; - } else if (scaled_y < -7.0) { - dnlcdf += 0.000471585349920831 * x2 * scaled_y - + 0.0296839305424034 * x2 - + 0.207402143352332 * scaled_y + 0.425316974683324; - } else if (scaled_y < -3.9) { - dnlcdf += -0.0006972280656443 * x2 * scaled_y - + 0.0068218494628567 * x2 - + 0.0585761964460277 * scaled_y + 0.1034397670201370; - } else if (scaled_y < -2.1) { - dnlcdf += -0.0018742199480885 * x2 * scaled_y - - 0.0097119598291202 * x2 - - 0.0170137970924080 * scaled_y - 0.0100428567412041; - } - } else { - dnlcdf = INFINITY; + /** Log Phi(x), with the original, unscaled argument. */ + inline double std_normal_lcdf_impl(double x) { + if (x <= -4.0 * M_SQRT2) { + const double r = 2.0 / (x * x); + return -(0.5 * x) * x - log(-x) - 0.91893853320467274178 + + log1p(r * std_normal_lcdf_tail_correction(r)); } + const double s = fabs(x) * M_SQRT1_2; + if (s < 0.46875) { + const double e = std_normal_erf_small(s); + return log1p(x < 0.0 ? -e : e) - M_LN2; + } + const double erfcx = std_normal_erfcx(s); + if (x < 0.0) { + return -(0.5 * x) * x - M_LN2 + log(erfcx); + } + return log1p(-0.5 * exp(-(0.5 * x) * x) * erfcx); + } - return dnlcdf; + /** Slope phi(x) / Phi(x) in the original units. */ + inline double std_normal_lcdf_derivative(double x) { + if (x <= -4.0 * M_SQRT2) { + const double r = 2.0 / (x * x); + return -x / (1.0 + r * std_normal_lcdf_tail_correction(r)); + } + const double s = fabs(x) * M_SQRT1_2; + if (s < 0.46875) { + const double e = std_normal_erf_small(s); + return (0.5 * M_SQRT2 * M_2_SQRTPI) * exp(-s * s) + / (x < 0.0 ? 1.0 - e : 1.0 + e); + } + const double erfcx = std_normal_erfcx(s); + if (x < 0.0) { + return (0.5 * M_SQRT2 * M_2_SQRTPI) / erfcx; + } + const double density = exp(-(0.5 * x) * x); + return (0.5 * M_SQRT1_2 * M_2_SQRTPI) * density + / (1.0 - 0.5 * density * erfcx); }) "\n#endif\n"; // NOLINT // \endcond diff --git a/stan/math/opencl/prim/exp_mod_normal_cdf.hpp b/stan/math/opencl/prim/exp_mod_normal_cdf.hpp index 9510340e2fa..9b2d4b40efb 100644 --- a/stan/math/opencl/prim/exp_mod_normal_cdf.hpp +++ b/stan/math/opencl/prim/exp_mod_normal_cdf.hpp @@ -3,30 +3,12 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include +#include namespace stan { namespace math { -/** \ingroup opencl - * Returns the double exponential cumulative density function. Given - * containers of matching sizes, returns the product of probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @tparam T_inv_scale_cl type of inverse scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @param lambda (Sequence of) inverse scale(s). - * @return The log of the product of densities. - */ template exp_mod_normal_cdf(const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma, const T_inv_scale_cl& lambda) { - static constexpr const char* function = "exp_mod_normal_cdf(OpenCL)"; - using T_partials_return - = partials_return_t; - using std::isfinite; - using std::isnan; - - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); - if (N == 0) { - return 1.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& mu_col = as_column_vector_or_scalar(mu); - const auto& sigma_col = as_column_vector_or_scalar(sigma); - const auto& lambda_col = as_column_vector_or_scalar(lambda); - - const auto& y_val = value_of(y_col); - const auto& mu_val = value_of(mu_col); - const auto& sigma_val = value_of(sigma_col); - const auto& lambda_val = value_of(lambda_col); - - auto check_y_not_nan - = check_cl(function, "Random variable", y_val, "not NaN"); - auto y_not_nan_expr = !isnan(y_val); - auto check_mu_finite - = check_cl(function, "Location parameter", mu_val, "finite"); - auto mu_finite_expr = isfinite(mu_val); - auto check_sigma_positive_finite - = check_cl(function, "Scale parameter", sigma_val, "positive finite"); - auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - auto check_lambda_positive_finite - = check_cl(function, "Inv_cale parameter", lambda_val, "positive finite"); - auto lambda_positive_finite_expr = 0 < lambda_val && isfinite(lambda_val); - - auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); - auto inv_sigma = elt_divide(1.0, sigma_val); - auto diff = y_val - mu_val; - auto v = elt_multiply(lambda_val, sigma_val); - auto scaled_diff = elt_multiply(diff, inv_sigma * INV_SQRT_TWO); - auto scaled_diff_diff = scaled_diff - v * INV_SQRT_TWO; - auto erf_calc = 0.5 * (1.0 + erf(scaled_diff_diff)); - auto exp_term = exp(0.5 * square(v) - elt_multiply(lambda_val, diff)); - auto cdf_n = 0.5 + 0.5 * erf(scaled_diff) - elt_multiply(exp_term, erf_calc); - auto cdf_expr = colwise_prod(cdf_n); - - auto exp_term_2 = exp(-square(scaled_diff_diff)); - auto deriv_1 = elt_multiply(elt_multiply(lambda_val, exp_term), erf_calc); - auto deriv_2 = INV_SQRT_TWO_PI - * elt_multiply(elt_multiply(exp_term, exp_term_2), inv_sigma); - auto deriv_3 - = INV_SQRT_TWO_PI * elt_multiply(exp(-square(scaled_diff)), inv_sigma); - auto mu_deriv1 = elt_divide(deriv_2 - deriv_1 - deriv_3, cdf_n); - auto sigma_deriv1 = elt_divide( - -elt_multiply(deriv_1 - deriv_2, v) - + elt_multiply(deriv_3 - deriv_2, scaled_diff) * SQRT_TWO, - cdf_n); - auto lambda_deriv1 = elt_divide( - elt_multiply( - exp_term, - INV_SQRT_TWO_PI * elt_multiply(sigma_val, exp_term_2) - - elt_multiply(elt_multiply(v, sigma_val) - diff, erf_calc)), - cdf_n); - - matrix_cl any_y_neg_inf_cl; - matrix_cl cdf_cl; - matrix_cl y_deriv_cl; - matrix_cl mu_deriv_cl; - matrix_cl sigma_deriv_cl; - matrix_cl lambda_deriv_cl; - - results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite, - check_lambda_positive_finite, any_y_neg_inf_cl, cdf_cl, y_deriv_cl, - mu_deriv_cl, sigma_deriv_cl, lambda_deriv_cl) - = expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr, - lambda_positive_finite_expr, any_y_neg_inf, cdf_expr, - calc_if>(cdf_n), - calc_if>(mu_deriv1), - calc_if>(sigma_deriv1), - calc_if>(lambda_deriv1)); - - if (from_matrix_cl(any_y_neg_inf_cl).maxCoeff()) { - return 0.0; - } - - T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod(); - - auto ops_partials - = make_partials_propagator(y_col, mu_col, sigma_col, lambda_col); - if constexpr (is_any_autodiff_v) { - auto mu_deriv = elt_multiply( - static_select::value>(0, mu_deriv_cl), - cdf); - auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_multiply( - static_select>(0, sigma_deriv_cl), cdf); - auto lambda_deriv = elt_multiply( - static_select>(0, lambda_deriv_cl), cdf); - - results(y_deriv_cl, mu_deriv_cl, sigma_deriv_cl, lambda_deriv_cl) - = expressions(calc_if>(y_deriv), - calc_if>(mu_deriv), - calc_if>(sigma_deriv), - calc_if>(lambda_deriv)); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(mu_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = std::move(sigma_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<3>(ops_partials) = std::move(lambda_deriv_cl); - } - } - return ops_partials.build(cdf); + return exp(internal::exp_mod_normal_lcdf_opencl_impl( + "exp_mod_normal_cdf(OpenCL)", y, mu, sigma, lambda)); } } // namespace math diff --git a/stan/math/opencl/prim/exp_mod_normal_lccdf.hpp b/stan/math/opencl/prim/exp_mod_normal_lccdf.hpp index 8f250a0d40f..a91b980e698 100644 --- a/stan/math/opencl/prim/exp_mod_normal_lccdf.hpp +++ b/stan/math/opencl/prim/exp_mod_normal_lccdf.hpp @@ -3,31 +3,11 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { -/** \ingroup opencl - * Returns the exp mod normal log complementary cumulative density - * function. Given containers of matching sizes, returns the log sum of - * probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @tparam T_inv_scale_cl type of inverse scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @param lambda (Sequence of) inverse scale(s). - * @return The log of the product of densities. - */ template exp_mod_normal_lccdf(const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma, const T_inv_scale_cl& lambda) { - static constexpr const char* function = "exp_mod_normal_lccdf(OpenCL)"; - using T_partials_return - = partials_return_t; - using std::isfinite; - using std::isnan; - - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); - if (N == 0) { - return 0.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& mu_col = as_column_vector_or_scalar(mu); - const auto& sigma_col = as_column_vector_or_scalar(sigma); - const auto& lambda_col = as_column_vector_or_scalar(lambda); - - const auto& y_val = value_of(y_col); - const auto& mu_val = value_of(mu_col); - const auto& sigma_val = value_of(sigma_col); - const auto& lambda_val = value_of(lambda_col); - - auto check_y_not_nan - = check_cl(function, "Random variable", y_val, "not NaN"); - auto y_not_nan_expr = !isnan(y_val); - auto check_mu_finite - = check_cl(function, "Location parameter", mu_val, "finite"); - auto mu_finite_expr = isfinite(mu_val); - auto check_sigma_positive_finite - = check_cl(function, "Scale parameter", sigma_val, "positive finite"); - auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - auto check_lambda_positive_finite - = check_cl(function, "Inv_cale parameter", lambda_val, "positive finite"); - auto lambda_positive_finite_expr = 0 < lambda_val && isfinite(lambda_val); - - auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); - auto any_y_pos_inf = colwise_max(cast(y_val == INFTY)); - auto inv_sigma = elt_divide(1.0, sigma_val); - auto diff = y_val - mu_val; - auto scaled_diff = elt_multiply(diff, inv_sigma * INV_SQRT_TWO); - auto v = elt_multiply(lambda_val, sigma_val); - auto scaled_diff_diff = scaled_diff - v * INV_SQRT_TWO; - auto erf_calc = 0.5 * (1.0 + erf(scaled_diff_diff)); - auto exp_term = exp(0.5 * square(v) - elt_multiply(lambda_val, diff)); - auto ccdf_n = 0.5 - 0.5 * erf(scaled_diff) + elt_multiply(exp_term, erf_calc); - auto ccdf_log_expr = colwise_sum(log(ccdf_n)); - - auto exp_term_2 = exp(-square(scaled_diff_diff)); - auto deriv_1 = elt_multiply(elt_multiply(lambda_val, exp_term), erf_calc); - auto deriv_2 = INV_SQRT_TWO_PI - * elt_multiply(elt_multiply(exp_term, exp_term_2), inv_sigma); - auto deriv_3 - = INV_SQRT_TWO_PI * elt_multiply(exp(-square(scaled_diff)), inv_sigma); - auto mu_deriv = elt_divide(deriv_1 - deriv_2 + deriv_3, ccdf_n); - auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_divide( - elt_multiply(deriv_1 - deriv_2, v) - + elt_multiply(deriv_3 - deriv_2, scaled_diff) * SQRT_TWO, - ccdf_n); - auto lambda_deriv = elt_divide( - elt_multiply(exp_term, - elt_multiply(elt_multiply(v, sigma_val) - diff, erf_calc) - - INV_SQRT_TWO_PI * elt_multiply(sigma_val, exp_term_2)), - ccdf_n); - - matrix_cl any_y_neg_inf_cl; - matrix_cl any_y_pos_inf_cl; - matrix_cl ccdf_log_cl; - matrix_cl mu_deriv_cl; - matrix_cl y_deriv_cl; - matrix_cl sigma_deriv_cl; - matrix_cl lambda_deriv_cl; - - results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite, - check_lambda_positive_finite, any_y_neg_inf_cl, any_y_pos_inf_cl, - ccdf_log_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl, lambda_deriv_cl) - = expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr, - lambda_positive_finite_expr, any_y_neg_inf, any_y_pos_inf, - ccdf_log_expr, calc_if>(y_deriv), - calc_if>(mu_deriv), - calc_if>(sigma_deriv), - calc_if>(lambda_deriv)); - - if (from_matrix_cl(any_y_pos_inf_cl).maxCoeff()) { - return NEGATIVE_INFTY; - } - - if (from_matrix_cl(any_y_neg_inf_cl).maxCoeff()) { - return 0.0; - } - - T_partials_return ccdf_log = (from_matrix_cl(ccdf_log_cl)).sum(); - - auto ops_partials - = make_partials_propagator(y_col, mu_col, sigma_col, lambda_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(mu_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = std::move(sigma_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<3>(ops_partials) = std::move(lambda_deriv_cl); - } - return ops_partials.build(ccdf_log); + return internal::exp_mod_normal_lcdf_opencl_impl( + "exp_mod_normal_lccdf(OpenCL)", y, mu, sigma, lambda); } } // namespace math diff --git a/stan/math/opencl/prim/exp_mod_normal_lcdf.hpp b/stan/math/opencl/prim/exp_mod_normal_lcdf.hpp index 14f67b9e5f3..4242c2dfaf9 100644 --- a/stan/math/opencl/prim/exp_mod_normal_lcdf.hpp +++ b/stan/math/opencl/prim/exp_mod_normal_lcdf.hpp @@ -8,45 +8,27 @@ #include #include #include -#include #include namespace stan { namespace math { +namespace internal { -/** \ingroup opencl - * Returns the exp mod normal log cumulative density - * function. Given containers of matching sizes, returns the log sum of - * probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @tparam T_inv_scale_cl type of inverse scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @param lambda (Sequence of) inverse scale(s). - * @return The log of the product of densities. - */ -template * = nullptr, - require_any_not_stan_scalar_t* = nullptr> +/** Same log-space formulation as the prim exp_mod_normal_lcdf_impl. */ +template inline return_type_t -exp_mod_normal_lcdf(const T_y_cl& y, const T_loc_cl& mu, - const T_scale_cl& sigma, const T_inv_scale_cl& lambda) { - static constexpr const char* function = "exp_mod_normal_lcdf(OpenCL)"; - using T_partials_return - = partials_return_t; +exp_mod_normal_lcdf_opencl_impl(const char* function, const T_y_cl& y, + const T_loc_cl& mu, const T_scale_cl& sigma, + const T_inv_scale_cl& lambda) { + constexpr double sign = upper ? -1.0 : 1.0; using std::isfinite; using std::isnan; check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); + mu, "Scale parameter", sigma, "Inv_scale parameter", + lambda); + const size_t N = max_size(y, mu, sigma, lambda); if (N == 0) { return 0.0; } @@ -70,89 +52,47 @@ exp_mod_normal_lcdf(const T_y_cl& y, const T_loc_cl& mu, auto check_sigma_positive_finite = check_cl(function, "Scale parameter", sigma_val, "positive finite"); auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - auto check_lambda_positive_finite - = check_cl(function, "Inv_cale parameter", lambda_val, "positive finite"); + auto check_lambda_positive_finite = check_cl(function, "Inv_scale parameter", + lambda_val, "positive finite"); auto lambda_positive_finite_expr = 0 < lambda_val && isfinite(lambda_val); auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto any_y_pos_inf = colwise_max(cast(y_val == INFTY)); - auto sigma_inv = elt_divide(1.0, sigma_val); auto diff = y_val - mu_val; - auto scaled_diff = elt_multiply(diff * INV_SQRT_TWO, sigma_inv); + auto z = elt_divide(diff, sigma_val); auto v = elt_multiply(lambda_val, sigma_val); - auto scaled_diff_diff = scaled_diff - v * INV_SQRT_TWO; - auto cdf_term_1 = 0.5 + 0.5 * erf(scaled_diff); - auto cdf_term_2_phi = 0.5 * (1.0 + erf(scaled_diff_diff)); - auto log_exp_term = 0.5 * square(v) - elt_multiply(lambda_val, diff); - auto exp_term = exp(log_exp_term); - auto cdf_term_2 = elt_multiply(exp_term, cdf_term_2_phi); - auto cdf_n = cdf_term_1 - cdf_term_2; - auto use_stable = cdf_n <= 0.0 || !isfinite(cdf_n); - - auto exp_term_2 = exp(-square(scaled_diff_diff)); - auto deriv_1 - = elt_multiply(elt_multiply(lambda_val, exp_term), cdf_term_2_phi); - auto deriv_2 = INV_SQRT_TWO_PI - * elt_multiply(elt_multiply(exp_term, exp_term_2), sigma_inv); - auto deriv_3 - = INV_SQRT_TWO_PI * elt_multiply(exp(-square(scaled_diff)), sigma_inv); - auto direct_cdf_log = log(cdf_n); - auto direct_y_deriv = elt_divide(deriv_1 - deriv_2 + deriv_3, cdf_n); - auto direct_mu_deriv = -direct_y_deriv; - auto direct_sigma_deriv = -elt_divide( - elt_multiply(deriv_1 - deriv_2, v) - + elt_multiply(deriv_3 - deriv_2, scaled_diff) * SQRT_TWO, - cdf_n); - auto direct_lambda_deriv = elt_divide( - elt_multiply(exp_term, - INV_SQRT_TWO_PI * elt_multiply(sigma_val, exp_term_2) - - elt_multiply(elt_multiply(v, sigma_val) - diff, - cdf_term_2_phi)), - cdf_n); - - auto log_cdf_term_1 = std_normal_lcdf_scaled_impl(scaled_diff); - auto dlog_cdf_term_1 = std_normal_lcdf_dscaled_impl(scaled_diff); - auto log_cdf_term_2_phi = std_normal_lcdf_scaled_impl(scaled_diff_diff); - auto dlog_cdf_term_2_phi = std_normal_lcdf_dscaled_impl(scaled_diff_diff); - auto log_cdf_term_2 = log_exp_term + log_cdf_term_2_phi; - auto log_cdf_n = log_diff_exp(log_cdf_term_1, log_cdf_term_2); - auto cdf_term_1_weight = exp(log_cdf_term_1 - log_cdf_n); - auto cdf_term_2_weight = exp(log_cdf_term_2 - log_cdf_n); - auto scaled_diff_deriv - = elt_multiply(dlog_cdf_term_1, sigma_inv * INV_SQRT_TWO); - auto scaled_diff_diff_deriv - = elt_multiply(dlog_cdf_term_2_phi, sigma_inv * INV_SQRT_TWO); - auto stable_y_deriv - = elt_multiply(cdf_term_1_weight, scaled_diff_deriv) - - elt_multiply(cdf_term_2_weight, -lambda_val + scaled_diff_diff_deriv); - auto stable_mu_deriv = -stable_y_deriv; - auto stable_sigma_deriv - = elt_multiply(cdf_term_1_weight, - -elt_multiply(dlog_cdf_term_1, - elt_multiply(scaled_diff, sigma_inv))) - - elt_multiply( - cdf_term_2_weight, - elt_multiply(lambda_val, v) - - elt_multiply( - dlog_cdf_term_2_phi, - elt_multiply(scaled_diff + v * INV_SQRT_TWO, sigma_inv))); - auto stable_lambda_deriv = -elt_multiply( - cdf_term_2_weight, - elt_multiply(v, sigma_val) - diff - - elt_multiply(dlog_cdf_term_2_phi, sigma_val * INV_SQRT_TWO)); - auto cdf_log_expr - = colwise_sum(select(use_stable, log_cdf_n, direct_cdf_log)); - auto y_deriv = select(use_stable, stable_y_deriv, direct_y_deriv); - auto mu_deriv = select(use_stable, stable_mu_deriv, direct_mu_deriv); - auto sigma_deriv = select(use_stable, stable_sigma_deriv, direct_sigma_deriv); - auto lambda_deriv - = select(use_stable, stable_lambda_deriv, direct_lambda_deriv); + auto z_a = sign * z; + auto z_b = z - v; + auto log_a = math::std_normal_lcdf_impl(z_a); + auto log_b = 0.5 * square(v) - elt_multiply(lambda_val, diff) + + math::std_normal_lcdf_impl(z_b); + auto lp = [&]() { + if constexpr (upper) { + return fmax(log_a, log_b) + log1p_exp(-fabs(log_a - log_b)); + } else { + return log_diff_exp(log_a, log_b); + } + }(); + auto cdf_log_expr = colwise_sum(lp); + + auto w_b = -sign * exp(log_b - lp); + auto s_b = elt_multiply(w_b, std_normal_lcdf_derivative(z_b)); + auto s = elt_divide( + sign * elt_multiply(exp(log_a - lp), std_normal_lcdf_derivative(z_a)) + + s_b, + sigma_val); + auto q = elt_multiply(w_b, v) - s_b; + auto y_deriv = s - elt_multiply(w_b, lambda_val); + auto mu_deriv = -y_deriv; + auto sigma_deriv + = select(s == 0, 0.0, -elt_multiply(s, z)) + elt_multiply(lambda_val, q); + auto lambda_deriv = elt_multiply(sigma_val, q) - elt_multiply(w_b, diff); matrix_cl any_y_neg_inf_cl; matrix_cl any_y_pos_inf_cl; matrix_cl cdf_log_cl; - matrix_cl mu_deriv_cl; matrix_cl y_deriv_cl; + matrix_cl mu_deriv_cl; matrix_cl sigma_deriv_cl; matrix_cl lambda_deriv_cl; @@ -166,15 +106,14 @@ exp_mod_normal_lcdf(const T_y_cl& y, const T_loc_cl& mu, calc_if>(sigma_deriv), calc_if>(lambda_deriv)); - if (from_matrix_cl(any_y_pos_inf_cl).maxCoeff()) { - return 0.0; - } - if (from_matrix_cl(any_y_neg_inf_cl).maxCoeff()) { - return NEGATIVE_INFTY; + return upper ? 0.0 : NEGATIVE_INFTY; + } + if (from_matrix_cl(any_y_pos_inf_cl).maxCoeff()) { + return upper ? NEGATIVE_INFTY : 0.0; } - T_partials_return cdf_log = (from_matrix_cl(cdf_log_cl)).sum(); + double cdf_log = sum(from_matrix_cl(cdf_log_cl)); auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col, lambda_col); @@ -194,6 +133,21 @@ exp_mod_normal_lcdf(const T_y_cl& y, const T_loc_cl& mu, return ops_partials.build(cdf_log); } +} // namespace internal + +template * = nullptr, + require_any_not_stan_scalar_t* = nullptr> +inline return_type_t +exp_mod_normal_lcdf(const T_y_cl& y, const T_loc_cl& mu, + const T_scale_cl& sigma, const T_inv_scale_cl& lambda) { + return internal::exp_mod_normal_lcdf_opencl_impl( + "exp_mod_normal_lcdf(OpenCL)", y, mu, sigma, lambda); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/opencl/prim/exp_mod_normal_lpdf.hpp b/stan/math/opencl/prim/exp_mod_normal_lpdf.hpp index 00bb0440cf8..7fd57590222 100644 --- a/stan/math/opencl/prim/exp_mod_normal_lpdf.hpp +++ b/stan/math/opencl/prim/exp_mod_normal_lpdf.hpp @@ -81,33 +81,29 @@ exp_mod_normal_lpdf(const T_y_cl& y, const T_loc_cl& mu, auto lambda_positive_finite_expr = isfinite(lambda_val) && lambda_val > 0; auto inv_sigma_expr = elt_divide(1.0, sigma_val); - auto sigma_sq_expr = elt_multiply(sigma_val, sigma_val); + auto sigma_sq_expr = square(sigma_val); auto lambda_sigma_sq_expr = elt_multiply(lambda_val, sigma_sq_expr); auto mu_minus_y_expr = mu_val - y_val; - auto inner_term_expr = elt_multiply(mu_minus_y_expr + lambda_sigma_sq_expr, - INV_SQRT_TWO * inv_sigma_expr); - auto erfc_calc_expr = erfc(inner_term_expr); + // log(erfc(t) / 2) = log Phi(-sqrt(2) t) cancels the log(1/2) constant. + auto z_expr + = -elt_multiply(mu_minus_y_expr + lambda_sigma_sq_expr, inv_sigma_expr); auto logp1_expr = elt_multiply(lambda_val, mu_minus_y_expr + 0.5 * lambda_sigma_sq_expr) - + log(erfc_calc_expr); + + math::std_normal_lcdf_impl(z_expr); auto logp_expr = colwise_sum( static_select::value>( logp1_expr + log(lambda_val), logp1_expr)); - auto deriv_logerfc_expr - = elt_divide(-SQRT_TWO_OVER_SQRT_PI - * exp(-elt_multiply(inner_term_expr, inner_term_expr)), - erfc_calc_expr); - auto deriv_expr - = lambda_val + elt_multiply(deriv_logerfc_expr, inv_sigma_expr); + auto slope_expr = std_normal_lcdf_derivative(z_expr); + auto deriv_expr = lambda_val - elt_multiply(slope_expr, inv_sigma_expr); auto deriv_sigma_expr - = elt_multiply(sigma_val, elt_multiply(lambda_val, lambda_val)) - + elt_multiply( - deriv_logerfc_expr, - (lambda_val - elt_divide(mu_minus_y_expr, sigma_sq_expr))); + = elt_multiply(sigma_val, square(lambda_val)) + - elt_multiply( + slope_expr, + lambda_val - elt_multiply(mu_minus_y_expr, square(inv_sigma_expr))); auto deriv_lambda_expr = elt_divide(1.0, lambda_val) + lambda_sigma_sq_expr + mu_minus_y_expr - + elt_multiply(deriv_logerfc_expr, sigma_val); + - elt_multiply(slope_expr, sigma_val); matrix_cl logp_cl; matrix_cl y_deriv_cl; @@ -126,9 +122,6 @@ exp_mod_normal_lpdf(const T_y_cl& y, const T_loc_cl& mu, calc_if>(deriv_lambda_expr)); T_partials_return logp = sum(from_matrix_cl(logp_cl)); - if constexpr (include_summand::value) { - logp -= LOG_TWO * N; - } auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col, lambda_col); diff --git a/stan/math/opencl/prim/lognormal_cdf.hpp b/stan/math/opencl/prim/lognormal_cdf.hpp index 19781c96639..6217b1d17da 100644 --- a/stan/math/opencl/prim/lognormal_cdf.hpp +++ b/stan/math/opencl/prim/lognormal_cdf.hpp @@ -3,29 +3,12 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include +#include namespace stan { namespace math { -/** \ingroup opencl - * Returns the loghormal cumulative distribution function for the given - * location, and scale. If given containers of matching sizes - * returns the product of probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @return The log of the product of densities. - */ template < typename T_y_cl, typename T_loc_cl, typename T_scale_cl, require_all_prim_or_rev_kernel_expression_t* = nullptr> inline return_type_t lognormal_cdf( const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = "lognormal_cdf(OpenCL)"; - using T_partials_return = partials_return_t; - using std::isfinite; - using std::isnan; - - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); - if (N == 0) { - return 1.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& mu_col = as_column_vector_or_scalar(mu); - const auto& sigma_col = as_column_vector_or_scalar(sigma); - - const auto& y_val = value_of(y_col); - const auto& mu_val = value_of(mu_col); - const auto& sigma_val = value_of(sigma_col); - - auto check_y_nonnegative - = check_cl(function, "Random variable", y_val, "nonnegative"); - auto y_nonnegative_expr = 0.0 <= y_val; - auto check_mu_finite - = check_cl(function, "Location parameter", mu_val, "finite"); - auto mu_finite_expr = isfinite(mu_val); - auto check_sigma_positive_finite - = check_cl(function, "Scale parameter", sigma_val, "positive finite"); - auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - - auto any_y_zero = colwise_max(cast(y_val == 0.0)); - auto log_y = log(y_val); - auto scaled_diff = elt_divide(log_y - mu_val, sigma_val * SQRT_TWO); - auto erfc_m_diff = erfc(-scaled_diff); - auto cdf_n = 0.5 * erfc_m_diff; - auto cdf_expr = colwise_prod(cdf_n); - auto mu_deriv_tmp - = -INV_SQRT_TWO_PI - * elt_divide(exp(-square(scaled_diff)), elt_multiply(sigma_val, cdf_n)); - auto y_deriv_tmp = elt_divide(-mu_deriv_tmp, y_val); - auto sigma_deriv_tmp = elt_multiply(mu_deriv_tmp, scaled_diff * SQRT_TWO); - - matrix_cl any_y_zero_cl; - matrix_cl cdf_cl; - matrix_cl mu_deriv_cl; - matrix_cl y_deriv_cl; - matrix_cl sigma_deriv_cl; - - results(check_y_nonnegative, check_mu_finite, check_sigma_positive_finite, - any_y_zero_cl, cdf_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl) - = expressions(y_nonnegative_expr, mu_finite_expr, - sigma_positive_finite_expr, any_y_zero, cdf_expr, - calc_if>(y_deriv_tmp), - calc_if>(mu_deriv_tmp), - calc_if>(sigma_deriv_tmp)); - - if (from_matrix_cl(any_y_zero_cl).maxCoeff()) { - return 0.0; - } - - T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod(); - - auto mu_deriv = mu_deriv_cl * cdf; - auto y_deriv = y_deriv_cl * cdf; - auto sigma_deriv = sigma_deriv_cl * cdf; - - results(mu_deriv_cl, y_deriv_cl, sigma_deriv_cl) - = expressions(calc_if>(mu_deriv), - calc_if>(y_deriv), - calc_if>(sigma_deriv)); - - auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(mu_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = std::move(sigma_deriv_cl); - } - return ops_partials.build(cdf); + return exp(internal::lognormal_lcdf_opencl_impl( + "lognormal_cdf(OpenCL)", y, mu, sigma)); } } // namespace math diff --git a/stan/math/opencl/prim/lognormal_lccdf.hpp b/stan/math/opencl/prim/lognormal_lccdf.hpp index af407a42967..76a0a1d58dd 100644 --- a/stan/math/opencl/prim/lognormal_lccdf.hpp +++ b/stan/math/opencl/prim/lognormal_lccdf.hpp @@ -3,29 +3,11 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { -/** \ingroup opencl - * Returns the lognormal log complementary cumulative distribution function - * for the given location, and scale. If given containers of matching sizes - * returns the log sum of probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @return The log of the product of densities. - */ template < typename T_y_cl, typename T_loc_cl, typename T_scale_cl, require_all_prim_or_rev_kernel_expression_t* = nullptr> inline return_type_t lognormal_lccdf( const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = "lognormal_lccdf(OpenCL)"; - using T_partials_return = partials_return_t; - using std::isfinite; - using std::isnan; - - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); - if (N == 0) { - return 0.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& mu_col = as_column_vector_or_scalar(mu); - const auto& sigma_col = as_column_vector_or_scalar(sigma); - - const auto& y_val = value_of(y_col); - const auto& mu_val = value_of(mu_col); - const auto& sigma_val = value_of(sigma_col); - - auto check_y_nonnegative - = check_cl(function, "Random variable", y_val, "nonnegative"); - auto y_nonnegative = 0.0 <= y_val; - auto check_mu_finite - = check_cl(function, "Location parameter", mu_val, "finite"); - auto mu_finite_expr = isfinite(mu_val); - auto check_sigma_positive_finite - = check_cl(function, "Scale parameter", sigma_val, "positive finite"); - auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - - auto any_y_zero = colwise_max(cast(y_val == 0.0)); - auto log_y = log(y_val); - auto scaled_diff = elt_divide(log_y - mu_val, sigma_val * SQRT_TWO); - auto erfc_calc = erfc(scaled_diff); - auto lccdf_expr = colwise_sum(log(erfc_calc)); - auto mu_deriv = elt_divide(SQRT_TWO_OVER_SQRT_PI * exp(-square(scaled_diff)), - elt_multiply(sigma_val, erfc_calc)); - auto y_deriv = elt_divide(mu_deriv, -y_val); - auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff * SQRT_TWO); - - matrix_cl any_y_zero_cl; - matrix_cl lccdf_cl; - matrix_cl y_deriv_cl; - matrix_cl mu_deriv_cl; - matrix_cl sigma_deriv_cl; - - results(check_y_nonnegative, check_mu_finite, check_sigma_positive_finite, - any_y_zero_cl, lccdf_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl) - = expressions(y_nonnegative, mu_finite_expr, sigma_positive_finite_expr, - any_y_zero, lccdf_expr, - calc_if>(y_deriv), - calc_if>(mu_deriv), - calc_if>(sigma_deriv)); - - if (from_matrix_cl(any_y_zero_cl).maxCoeff()) { - return 0.0; - } - - T_partials_return lccdf = N * LOG_HALF + sum(from_matrix_cl(lccdf_cl)); - - auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(mu_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = std::move(sigma_deriv_cl); - } - return ops_partials.build(lccdf); + return internal::lognormal_lcdf_opencl_impl("lognormal_lccdf(OpenCL)", + y, mu, sigma); } } // namespace math diff --git a/stan/math/opencl/prim/lognormal_lcdf.hpp b/stan/math/opencl/prim/lognormal_lcdf.hpp index 3f977a24e0a..31d79939311 100644 --- a/stan/math/opencl/prim/lognormal_lcdf.hpp +++ b/stan/math/opencl/prim/lognormal_lcdf.hpp @@ -12,29 +12,13 @@ namespace stan { namespace math { +namespace internal { -/** \ingroup opencl - * Returns the lognormal log cumulative distribution function - * for the given location, and scale. If given containers of matching sizes - * returns the log sum of probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @return The log of the product of densities. - */ -template < - typename T_y_cl, typename T_loc_cl, typename T_scale_cl, - require_all_prim_or_rev_kernel_expression_t* = nullptr, - require_any_not_stan_scalar_t* = nullptr> -inline return_type_t lognormal_lcdf( - const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = "lognormal_lcdf(OpenCL)"; - using T_partials_return = partials_return_t; +template +inline return_type_t lognormal_lcdf_opencl_impl( + const char* function, const T_y_cl& y, const T_loc_cl& mu, + const T_scale_cl& sigma) { + constexpr double sign = reflect ? -1.0 : 1.0; using std::isfinite; using std::isnan; @@ -64,14 +48,13 @@ inline return_type_t lognormal_lcdf( auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); auto any_y_zero = colwise_max(cast(y_val == 0.0)); - auto log_y = log(y_val); - auto scaled_diff = elt_divide(log_y - mu_val, sigma_val * SQRT_TWO); - auto erfc_calc = erfc(-scaled_diff); - auto lcdf_expr = colwise_sum(log(erfc_calc)); - auto mu_deriv = elt_divide(-SQRT_TWO_OVER_SQRT_PI * exp(-square(scaled_diff)), - elt_multiply(sigma_val, erfc_calc)); - auto y_deriv = elt_divide(mu_deriv, -y_val); - auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff * SQRT_TWO); + auto z = sign * elt_divide(log(y_val) - mu_val, sigma_val); + auto lcdf_expr = colwise_sum(math::std_normal_lcdf_impl(z)); + auto slope = std_normal_lcdf_derivative(z); + auto scaled_slope = elt_divide(slope, sigma_val); + auto y_deriv = sign * elt_divide(scaled_slope, y_val); + auto mu_deriv = -sign * scaled_slope; + auto sigma_deriv = select(slope == 0, 0.0, -elt_multiply(scaled_slope, z)); matrix_cl any_y_zero_cl; matrix_cl lcdf_cl; @@ -88,10 +71,10 @@ inline return_type_t lognormal_lcdf( calc_if>(sigma_deriv)); if (from_matrix_cl(any_y_zero_cl).maxCoeff()) { - return NEGATIVE_INFTY; + return reflect ? 0.0 : NEGATIVE_INFTY; } - T_partials_return lcdf = N * LOG_HALF + sum(from_matrix_cl(lcdf_cl)); + double lcdf = sum(from_matrix_cl(lcdf_cl)); auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col); @@ -107,6 +90,19 @@ inline return_type_t lognormal_lcdf( return ops_partials.build(lcdf); } +} // namespace internal + +template < + typename T_y_cl, typename T_loc_cl, typename T_scale_cl, + require_all_prim_or_rev_kernel_expression_t* = nullptr, + require_any_not_stan_scalar_t* = nullptr> +inline return_type_t lognormal_lcdf( + const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { + return internal::lognormal_lcdf_opencl_impl("lognormal_lcdf(OpenCL)", + y, mu, sigma); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/opencl/prim/normal_cdf.hpp b/stan/math/opencl/prim/normal_cdf.hpp index 2e764e8eb81..f5e9d06cc88 100644 --- a/stan/math/opencl/prim/normal_cdf.hpp +++ b/stan/math/opencl/prim/normal_cdf.hpp @@ -3,20 +3,16 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include +#include namespace stan { namespace math { /** \ingroup opencl * Returns the normal cumulative distribution function for the given - * location, and scale. If given containers of matching sizes - * returns the product of probabilities. + * location and scale. If given containers of matching sizes, returns the + * product of probabilities. * * @tparam T_y_cl type of scalar outcome * @tparam T_loc_cl type of location @@ -24,7 +20,7 @@ namespace math { * @param y (Sequence of) scalar(s). * @param mu (Sequence of) location(s). * @param sigma (Sequence of) scale(s). - * @return The log of the product of densities. + * @return The product of cumulative probabilities. */ template < typename T_y_cl, typename T_loc_cl, typename T_scale_cl, @@ -33,84 +29,8 @@ template < require_any_not_stan_scalar_t* = nullptr> inline return_type_t normal_cdf( const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = "normal_cdf(OpenCL)"; - using T_partials_return = partials_return_t; - using std::isfinite; - using std::isnan; - - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - const size_t N = max_size(y, mu, sigma); - if (N == 0) { - return 1.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& mu_col = as_column_vector_or_scalar(mu); - const auto& sigma_col = as_column_vector_or_scalar(sigma); - - const auto& y_val = value_of(y_col); - const auto& mu_val = value_of(mu_col); - const auto& sigma_val = value_of(sigma_col); - - auto check_y_not_nan - = check_cl(function, "Random variable", y_val, "not NaN"); - auto y_not_nan_expr = !isnan(y_val); - auto check_mu_finite - = check_cl(function, "Location parameter", mu_val, "finite"); - auto mu_finite_expr = isfinite(mu_val); - auto check_sigma_positive - = check_cl(function, "Scale parameter", sigma_val, "positive"); - auto sigma_positive_expr = 0 < sigma_val; - - auto scaled_diff = elt_divide(y_val - mu_val, sigma_val * SQRT_TWO); - auto cdf_n = select( - scaled_diff < -37.5 * INV_SQRT_TWO, 0.0, - select(scaled_diff < -5.0 * INV_SQRT_TWO, 0.5 * erfc(-scaled_diff), - select(scaled_diff > 8.25 * INV_SQRT_TWO, 1.0, - 0.5 * (1.0 + erf(scaled_diff))))); - auto cdf_expr = colwise_prod(cdf_n); - auto mu_deriv_tmp = select(scaled_diff < -37.5 * INV_SQRT_TWO, 0.0, - INV_SQRT_TWO_PI - * elt_divide(exp(-square(scaled_diff)), - elt_multiply(cdf_n, sigma_val))); - auto sigma_deriv_tmp = elt_multiply(mu_deriv_tmp, scaled_diff); - - matrix_cl cdf_cl; - matrix_cl mu_deriv_cl; - matrix_cl y_deriv_cl; - matrix_cl sigma_deriv_cl; - - results(check_y_not_nan, check_mu_finite, check_sigma_positive, cdf_cl, - mu_deriv_cl, sigma_deriv_cl) - = expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_expr, - cdf_expr, - calc_if>(mu_deriv_tmp), - calc_if>(sigma_deriv_tmp)); - - T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod(); - - auto y_deriv = elt_multiply(mu_deriv_cl, cdf); - auto mu_deriv = -y_deriv; - auto sigma_deriv = elt_multiply(sigma_deriv_cl, -SQRT_TWO * cdf); - - results(mu_deriv_cl, y_deriv_cl, sigma_deriv_cl) - = expressions(calc_if>(mu_deriv), - calc_if>(y_deriv), - calc_if>(sigma_deriv)); - - auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(mu_deriv_cl); - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = std::move(sigma_deriv_cl); - } - return ops_partials.build(cdf); + return exp(internal::normal_lcdf_opencl_impl("normal_cdf(OpenCL)", y, + mu, sigma)); } } // namespace math diff --git a/stan/math/opencl/prim/normal_lccdf.hpp b/stan/math/opencl/prim/normal_lccdf.hpp index 25f4f6ed0d8..ff539c6ac66 100644 --- a/stan/math/opencl/prim/normal_lccdf.hpp +++ b/stan/math/opencl/prim/normal_lccdf.hpp @@ -6,23 +6,7 @@ namespace stan { namespace math { -namespace internal { -constexpr char normal_lccdf_opencl_func[] = "normal_lccdf(OpenCL)"; -} // namespace internal -/** \ingroup opencl - * Returns the normal log complementary cumulative distribution function - * for the given location, and scale. If given containers of matching sizes - * returns the log sum of probabilities. - * - * @tparam T_y_cl type of scalar outcome - * @tparam T_loc_cl type of location - * @tparam T_scale_cl type of scale - * @param y (Sequence of) scalar(s). - * @param mu (Sequence of) location(s). - * @param sigma (Sequence of) scale(s). - * @return The log of the product of densities. - */ template < typename T_y_cl, typename T_loc_cl, typename T_scale_cl, require_all_prim_or_rev_kernel_expression_t* = nullptr> inline return_type_t normal_lccdf( const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - return normal_lcdf(-y, -mu, sigma); + return internal::normal_lcdf_opencl_impl("normal_lccdf(OpenCL)", y, mu, + sigma); } } // namespace math diff --git a/stan/math/opencl/prim/normal_lcdf.hpp b/stan/math/opencl/prim/normal_lcdf.hpp index 187eda05866..bc7361276c3 100644 --- a/stan/math/opencl/prim/normal_lcdf.hpp +++ b/stan/math/opencl/prim/normal_lcdf.hpp @@ -13,189 +13,12 @@ namespace stan { namespace math { namespace internal { -constexpr char normal_lcdf_opencl_func[] = "normal_lcdf(OpenCL)"; -const char opencl_normal_lcdf_impl[] = STRINGIFY( - double x2 = normal_lcdf_scaled_diff * normal_lcdf_scaled_diff; - double normal_lcdf_n = 0; - // Rigorous numerical approximations are applied here to deal with values - // of |normal_lcdf_scaled_diff|>>0. This is needed to deal with rare - // base-rate logistic regression problems where it is useful to use an - // alternative link function instead. - // - // use erfc() instead of erf() in order to retain precision - // since for x>0 erfc()->0 - if (normal_lcdf_scaled_diff > 0.0) { - // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) - normal_lcdf_n = log1p(-0.5 * erfc(normal_lcdf_scaled_diff)); - if (isnan(normal_lcdf_n)) { - normal_lcdf_n = 0; - } - } else if (normal_lcdf_scaled_diff > -4.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x) - normal_lcdf_n = log(erfc(-normal_lcdf_scaled_diff)) - M_LN2; - } else if (10.0 * log(fabs(normal_lcdf_scaled_diff)) < log(DBL_MAX)) { - // entering territory where erfc(-x)~0 - // need to use direct numerical approximation of normal_lcdf_n instead - // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) - // CDF(x) = 1/2erfc(-x) - double x4 = pow(normal_lcdf_scaled_diff, 4); - double x6 = pow(normal_lcdf_scaled_diff, 6); - double x8 = pow(normal_lcdf_scaled_diff, 8); - double x10 = pow(normal_lcdf_scaled_diff, 10); - double temp_p - = 0.000658749161529837803157 + 0.0160837851487422766278 / x2 - + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6 - + 0.305326634961232344035 / x8 + 0.0163153871373020978498 / x10; - double temp_q = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 - - 0.527905102951428412248 / x4 - - 1.87295284992346047209 / x6 - - 2.56852019228982242072 / x8 - 1.0 / x10; - normal_lcdf_n = -M_LN2 + log(0.5 * M_2_SQRTPI + (temp_p / temp_q) / x2) - - log(-normal_lcdf_scaled_diff) - x2; - } else { - // normal_lcdf_scaled_diff^10 term will overflow - normal_lcdf_n = -INFINITY; - }); -// NOLINTBEGIN -const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY( - double normal_ldncdf = 0.0; double t = 0.0; double t2 = 0.0; - double t4 = 0.0; double normal_lcdf_exp_m_x2 = 0.0; - double normal_lcdf_inv_x2 = 0.0; - // calculate using piecewise function - // (due to instability / inaccuracy in the various approximations) - if (normal_lcdf_deriv_scaled_diff > 2.9) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - t = 1.0 / (1.0 + 0.3275911 * normal_lcdf_deriv_scaled_diff); - t2 = t * t; - t4 = pow(t, 4); - // A&S 7.1.26 keeps exp(-x2) in the numerator, as R's pnorm do_del - // does; refs in stan/math/prim/prob/normal_lcdf.hpp - normal_lcdf_exp_m_x2 = exp(-x2); - normal_ldncdf - = 0.5 * M_2_SQRTPI * normal_lcdf_exp_m_x2 - / (1.0 - - normal_lcdf_exp_m_x2 - * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 - - 1.453152027 * t2 * t + 1.061405429 * t4)); - } else if (normal_lcdf_deriv_scaled_diff > 2.5) { - // in the trouble area where all of the standard numerical - // approximations are unstable - bridge the gap using Taylor - // expansions of the analytic function - // use Taylor expansion centred around x=2.7 - t = normal_lcdf_deriv_scaled_diff - 2.7; - t2 = t * t; - t4 = pow(t, 4); - normal_ldncdf = 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2 - - 0.008029540137 * t2 * t + 0.008232190507 * t4 - - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6); - } else if (normal_lcdf_deriv_scaled_diff > 2.1) { - // use Taylor expansion centred around x=2.3 - t = normal_lcdf_deriv_scaled_diff - 2.3; - t2 = t * t; - t4 = pow(t, 4); - normal_ldncdf = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 - - 0.03326906904 * t2 * t + 0.02482478940 * t4 - - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); - } else if (normal_lcdf_deriv_scaled_diff > 1.5) { - // use Taylor expansion centred around x=1.85 - t = normal_lcdf_deriv_scaled_diff - 1.85; - t2 = t * t; - t4 = pow(t, 4); - normal_ldncdf = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 - - 0.09274533184 * t2 * t + 0.03543327418 * t4 - + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); - } else if (normal_lcdf_deriv_scaled_diff > 0.8) { - // use Taylor expansion centred around x=1.15 - t = normal_lcdf_deriv_scaled_diff - 1.15; - t2 = t * t; - t4 = pow(t, 4); - normal_ldncdf = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 - - 0.09748053605 * t2 * t - 0.04347986191 * t4 - + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); - } else if (normal_lcdf_deriv_scaled_diff > 0.1) { - // use Taylor expansion centred around x=0.45 - t = normal_lcdf_deriv_scaled_diff - 0.45; - t2 = t * t; - t4 = pow(t, 4); - normal_ldncdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 - + 0.04700850676 * t2 * t - 0.03478651979 * t4 - - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); - } else if (normal_lcdf_deriv_scaled_diff < -29.0) { - // asymptotic Mills ratio, DLMF 7.12.1; grows linearly as -2*scaled_diff, - // so no quadratic fit can track it. Same 1/x^2 series shape as R's pnorm - normal_lcdf_inv_x2 = 1.0 / x2; - normal_ldncdf - = -2.0 * normal_lcdf_deriv_scaled_diff - / (1.0 - + normal_lcdf_inv_x2 - * (-0.5 - + normal_lcdf_inv_x2 - * (0.75 + normal_lcdf_inv_x2 * -1.875))); - } else if (10.0 * log(fabs(normal_lcdf_deriv_scaled_diff)) < log(DBL_MAX)) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - // use fact that erf(x)=-erf(-x) - // Abramowitz and Stegun define this for -inf* = nullptr, - require_any_not_stan_scalar_t* = nullptr> -inline return_type_t normal_lcdf( - const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = func; +template +inline return_type_t normal_lcdf_opencl_impl( + const char* function, const T_y_cl& y, const T_loc_cl& mu, + const T_scale_cl& sigma) { + constexpr double sign = reflect ? -1.0 : 1.0; using std::isfinite; using std::isnan; @@ -224,22 +47,13 @@ inline return_type_t normal_lcdf( = check_cl(function, "Scale parameter", sigma_val, "positive"); auto sigma_positive_expr = 0 < sigma_val; - auto scaled_diff = elt_divide(y_val - mu_val, sigma_val * SQRT_TWO); - - auto sigma_sqrt2 = sigma_val * SQRT_TWO; - - auto lcdf_n = opencl_code( - std::make_tuple("normal_lcdf_scaled_diff"), scaled_diff) - .template output("normal_lcdf_n"); - auto lcdf_expr = colwise_sum(lcdf_n); - - auto ldncdf - = opencl_code( - std::make_tuple("normal_lcdf_deriv_scaled_diff"), scaled_diff) - .template output("normal_ldncdf"); - auto y_deriv = elt_divide(ldncdf, sigma_sqrt2); + auto z = sign * elt_divide(y_val - mu_val, sigma_val); + auto lcdf_expr = colwise_sum(math::std_normal_lcdf_impl(z)); + auto slope = std_normal_lcdf_derivative(z); + auto scaled_slope = elt_divide(slope, sigma_val); + auto y_deriv = sign * scaled_slope; auto mu_deriv = -y_deriv; - auto sigma_deriv = -elt_divide(elt_multiply(ldncdf, scaled_diff), sigma_val); + auto sigma_deriv = select(slope == 0, 0.0, -elt_multiply(scaled_slope, z)); matrix_cl lcdf_cl; matrix_cl y_deriv_cl; @@ -269,6 +83,32 @@ inline return_type_t normal_lcdf( return ops_partials.build(lcdf); } +} // namespace internal + +/** \ingroup opencl + * Returns the normal log cumulative distribution function for the given + * location and scale. If given containers of matching sizes, returns the + * sum of log probabilities. + * + * @tparam T_y_cl type of scalar outcome + * @tparam T_loc_cl type of location + * @tparam T_scale_cl type of scale + * @param y (Sequence of) scalar(s). + * @param mu (Sequence of) location(s). + * @param sigma (Sequence of) scale(s). + * @return The log of the product of cumulative probabilities. + */ +template < + typename T_y_cl, typename T_loc_cl, typename T_scale_cl, + require_all_prim_or_rev_kernel_expression_t* = nullptr, + require_any_not_stan_scalar_t* = nullptr> +inline return_type_t normal_lcdf( + const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { + return internal::normal_lcdf_opencl_impl("normal_lcdf(OpenCL)", y, mu, + sigma); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/opencl/prim/skew_normal_lpdf.hpp b/stan/math/opencl/prim/skew_normal_lpdf.hpp index 8a02bd379be..64c98fefaa3 100644 --- a/stan/math/opencl/prim/skew_normal_lpdf.hpp +++ b/stan/math/opencl/prim/skew_normal_lpdf.hpp @@ -85,8 +85,8 @@ inline return_type_t skew_normal_lpdf( auto inv_sigma = elt_divide(1., sigma_val); auto y_minus_mu_over_sigma = elt_multiply((y_val - mu_val), inv_sigma); - auto log_erfc_alpha_z = log( - erfc(elt_multiply(alpha_val, y_minus_mu_over_sigma) * -INV_SQRT_TWO)); + auto alpha_z = elt_multiply(alpha_val, y_minus_mu_over_sigma); + auto log_erfc_alpha_z = LOG_TWO + std_normal_lcdf_impl(alpha_z); auto logp1 = log_erfc_alpha_z; auto logp2 = static_select::value>( @@ -94,14 +94,9 @@ inline return_type_t skew_normal_lpdf( auto logp_expr = colwise_sum( static_select< include_summand::value>( - logp2 - - elt_multiply(y_minus_mu_over_sigma, y_minus_mu_over_sigma) - * 0.5, - logp2)); - - auto scaled = elt_multiply(alpha_val, y_minus_mu_over_sigma) * INV_SQRT_TWO; - auto deriv_logerf = SQRT_TWO_OVER_SQRT_PI - * exp(-elt_multiply(scaled, scaled) - log_erfc_alpha_z); + logp2 - 0.5 * square(y_minus_mu_over_sigma), logp2)); + + auto deriv_logerf = std_normal_lcdf_derivative(alpha_z); auto y_loc_deriv = elt_multiply( y_minus_mu_over_sigma - elt_multiply(deriv_logerf, alpha_val), inv_sigma); auto sigma_deriv diff --git a/stan/math/opencl/prim/std_normal_cdf.hpp b/stan/math/opencl/prim/std_normal_cdf.hpp index 5a0cd0c1316..d3b19792cdb 100644 --- a/stan/math/opencl/prim/std_normal_cdf.hpp +++ b/stan/math/opencl/prim/std_normal_cdf.hpp @@ -3,68 +3,26 @@ #ifdef STAN_OPENCL #include -#include -#include -#include -#include -#include -#include +#include +#include namespace stan { namespace math { /** \ingroup opencl - * Returns the standard normal cumulative distribution function. + * Returns the standard normal cumulative distribution function. If given a + * container, returns the product of probabilities. * * @tparam T_y_cl type of scalar outcome * @param y (Sequence of) scalar(s). - * @return The log of the product of densities. + * @return The product of cumulative probabilities. */ template * = nullptr, require_any_not_stan_scalar_t* = nullptr> inline return_type_t std_normal_cdf(const T_y_cl& y) { - static constexpr const char* function = "std_normal_cdf(OpenCL)"; - using T_partials_return = partials_return_t; - using std::isfinite; - using std::isnan; - - const size_t N = math::size(y); - if (N == 0) { - return 1.0; - } - - const auto& y_col = as_column_vector_or_scalar(y); - const auto& y_val = value_of(y_col); - - auto check_y_not_nan - = check_cl(function, "Random variable", y_val, "not NaN"); - auto y_not_nan_expr = !isnan(y_val); - - auto scaled_y = y_val * INV_SQRT_TWO; - auto cdf_n - = select(y_val < -37.5, 0.0, - select(y_val < -5.0, 0.5 * erfc(-scaled_y), - select(y_val > 8.25, 1.0, 0.5 * (1.0 + erf(scaled_y))))); - auto cdf_expr = colwise_prod(cdf_n); - auto y_deriv1 - = select(y_val < -37.5, 0.0, - INV_SQRT_TWO_PI * elt_divide(exp(-square(scaled_y)), cdf_n)); - - matrix_cl cdf_cl; - matrix_cl y_deriv_cl; - - results(check_y_not_nan, cdf_cl, y_deriv_cl) = expressions( - y_not_nan_expr, cdf_expr, calc_if>(y_deriv1)); - - T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod(); - - auto ops_partials = make_partials_propagator(y_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = y_deriv_cl * cdf; - } - return ops_partials.build(cdf); + return exp(internal::std_normal_lcdf_opencl_impl( + "std_normal_cdf(OpenCL)", y)); } } // namespace math diff --git a/stan/math/opencl/prim/std_normal_lccdf.hpp b/stan/math/opencl/prim/std_normal_lccdf.hpp index 0804d5d9ad6..5eb3a1cf439 100644 --- a/stan/math/opencl/prim/std_normal_lccdf.hpp +++ b/stan/math/opencl/prim/std_normal_lccdf.hpp @@ -6,23 +6,13 @@ namespace stan { namespace math { -namespace internal { -constexpr char std_normal_lccdf_opencl_func[] = "std_normal_lccdf(OpenCL)"; -} // namespace internal -/** \ingroup opencl - * Returns the log standard normal complementary cumulative distribution - * function. - * - * @tparam T_y_cl type of scalar outcome - * @param y (Sequence of) scalar(s). - * @return The log of the product of densities. - */ template * = nullptr, require_any_not_stan_scalar_t* = nullptr> inline return_type_t std_normal_lccdf(const T_y_cl& y) { - return std_normal_lcdf(-y); + return internal::std_normal_lcdf_opencl_impl("std_normal_lccdf(OpenCL)", + y); } } // namespace math diff --git a/stan/math/opencl/prim/std_normal_lcdf.hpp b/stan/math/opencl/prim/std_normal_lcdf.hpp index 463b59b6f49..86c182950cb 100644 --- a/stan/math/opencl/prim/std_normal_lcdf.hpp +++ b/stan/math/opencl/prim/std_normal_lcdf.hpp @@ -14,28 +14,16 @@ namespace stan { namespace math { namespace internal { -constexpr char std_normal_lcdf_opencl_func[] = "std_normal_lcdf(OpenCL)"; -} // namespace internal -/** \ingroup opencl - * Returns the log standard normal complementary cumulative distribution - * function. - * - * @tparam T_y_cl type of scalar outcome - * @param y (Sequence of) scalar(s). - * @return The log of the product of densities. - */ -template * = nullptr, - require_any_not_stan_scalar_t* = nullptr> -inline return_type_t std_normal_lcdf(const T_y_cl& y) { - static constexpr const char* function = func; +template +inline return_type_t std_normal_lcdf_opencl_impl(const char* function, + const T_y_cl& y) { + constexpr double sign = reflect ? -1.0 : 1.0; using std::isfinite; using std::isnan; const size_t N = math::size(y); if (N == 0) { - return 1.0; + return 0.0; } const auto& y_col = as_column_vector_or_scalar(y); @@ -45,10 +33,9 @@ inline return_type_t std_normal_lcdf(const T_y_cl& y) { = check_cl(function, "Random variable", y_val, "not NaN"); auto y_not_nan_expr = !isnan(y_val); - auto scaled_y = y_val * INV_SQRT_TWO; - auto lcdf_expr = colwise_sum(std_normal_lcdf_scaled_impl(scaled_y)); - auto dnlcdf = std_normal_lcdf_dscaled_impl(scaled_y); - auto y_deriv = dnlcdf * INV_SQRT_TWO; + auto z = sign * y_val; + auto lcdf_expr = colwise_sum(math::std_normal_lcdf_impl(z)); + auto y_deriv = sign * std_normal_lcdf_derivative(z); matrix_cl lcdf_cl; matrix_cl y_deriv_cl; @@ -66,6 +53,24 @@ inline return_type_t std_normal_lcdf(const T_y_cl& y) { return ops_partials.build(lcdf); } +} // namespace internal + +/** \ingroup opencl + * Returns the log standard normal cumulative distribution + * function. + * + * @tparam T_y_cl type of scalar outcome + * @param y (Sequence of) scalar(s). + * @return The log of the product of cumulative probabilities. + */ +template * = nullptr, + require_any_not_stan_scalar_t* = nullptr> +inline return_type_t std_normal_lcdf(const T_y_cl& y) { + return internal::std_normal_lcdf_opencl_impl("std_normal_lcdf(OpenCL)", + y); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/fun/Phi.hpp b/stan/math/prim/fun/Phi.hpp index 0ae82a61755..fc51df8f448 100644 --- a/stan/math/prim/fun/Phi.hpp +++ b/stan/math/prim/fun/Phi.hpp @@ -3,10 +3,8 @@ #include #include -#include -#include -#include -#include +#include +#include #include namespace stan { @@ -24,22 +22,12 @@ namespace math { * This function can be used to implement the inverse link function * for probit regression. * - * Phi will underflow to 0 below -37.5 and overflow to 1 above 8 - * * @param x Argument. * @return Probability random sample is less than or equal to argument. */ inline double Phi(double x) { check_not_nan("Phi", "x", x); - if (x < -37.5) { - return 0; - } else if (x < -5.0) { - return 0.5 * erfc(-INV_SQRT_TWO * x); - } else if (x > 8.25) { - return 1; - } else { - return 0.5 * (1.0 + erf(INV_SQRT_TWO * x)); - } + return exp(internal::std_normal_lcdf_value_grad(x).first); } /** diff --git a/stan/math/prim/fun/std_normal_lcdf_impl.hpp b/stan/math/prim/fun/std_normal_lcdf_impl.hpp new file mode 100644 index 00000000000..3149f2e17d6 --- /dev/null +++ b/stan/math/prim/fun/std_normal_lcdf_impl.hpp @@ -0,0 +1,164 @@ +#ifndef STAN_MATH_PRIM_FUN_STD_NORMAL_LCDF_IMPL_HPP +#define STAN_MATH_PRIM_FUN_STD_NORMAL_LCDF_IMPL_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { +namespace internal { + +// Rational approximations from Cody (1969), Math. Comp. 23:631-637 (SPECFUN +// CALERF). Polynomial arithmetic only: no libm erf or erfc on the hot path. + +/** erf(x) for |x| < 0.46875. */ +template +inline T std_normal_erf_small(const T& x) { + static constexpr double a[] + = {3.16112374387056560, 1.13864154151050156e2, 3.77485237685302021e2, + 3.20937758913846947e3, 1.85777706184603153e-1}; + static constexpr double b[] = {2.36012909523441209e1, 2.44024637934444173e2, + 1.28261652607737228e3, 2.84423683343917062e3}; + const T x2 = square(x); + T numerator = a[4] * x2; + T denominator = x2; + for (int i = 0; i < 3; ++i) { + numerator = (numerator + a[i]) * x2; + denominator = (denominator + b[i]) * x2; + } + return x * (numerator + a[3]) / (denominator + b[3]); +} + +/** Tail factor in r = 1/x^2: erfcx(x) = (1 + r * correction) / (sqrt(pi) x). */ +template +inline T std_normal_tail_correction(const T& r) { + static constexpr double p[] + = {0.000658749161529837803157, 0.0160837851487422766278, + 0.125781726111229246204, 0.360344899949804439429, + 0.305326634961232344035, 0.0163153871373020978498}; + static constexpr double q[] + = {-0.00233520497626869185443, -0.0605183413124413191178, + -0.527905102951428412248, -1.87295284992346047209, + -2.56852019228982242072, -1.0}; + T numerator = p[5] * r + p[4]; + T denominator = q[5] * r + q[4]; + for (int i = 3; i >= 0; --i) { + numerator = numerator * r + p[i]; + denominator = denominator * r + q[i]; + } + return (numerator / denominator) / INV_SQRT_PI; +} + +/** erfcx(x) = exp(x^2) erfc(x) for x >= 0.46875. */ +template +inline T std_normal_erfcx(const T& x) { + if (x > 4) { + const T r = inv_square(x); + return (1 + r * std_normal_tail_correction(r)) * INV_SQRT_PI / x; + } + static constexpr double c[] + = {5.64188496988670089e-1, 8.88314979438837594, 6.61191906371416295e1, + 2.98635138197400131e2, 8.81952221241769090e2, 1.71204761263407058e3, + 2.05107837782607147e3, 1.23033935479799725e3, 2.15311535474403846e-8}; + static constexpr double d[] + = {1.57449261107098347e1, 1.17693950891312499e2, 5.37181101862009858e2, + 1.62138957456669019e3, 3.29079923573345963e3, 4.36261909014324716e3, + 3.43936767414372164e3, 1.23033935480374942e3}; + T numerator = c[8] * x; + T denominator = x; + for (int i = 0; i < 7; ++i) { + numerator = (numerator + c[i]) * x; + denominator = (denominator + d[i]) * x; + } + return (numerator + c[7]) / (denominator + d[7]); +} + +/** Scalar log Phi(z) and its slope phi(z) / Phi(z). + * For z < 0 both come from erfcx with no exp; z > 0 needs one exp for the + * complement. Infinite z gives -inf/0 values and inf/0 slopes. + * The lower tail keeps log1p, r = 2 (1/a)^2 and a / (1 + r c) so nested + * autodiff neither overflows nor loses the 2 / a^3 third derivative; the + * z > 40 return keeps derivatives finite when x^2 overflows. + */ +template * = nullptr> +inline std::pair, return_type_t> std_normal_lcdf_value_grad( + const T& z_in) { + using R = return_type_t; + const R z = z_in; + if (z > 40) { + return {0, 0}; + } + if (z <= -4 * SQRT_TWO) { + const R a = -z; + const R r = 2 * square(inv(a)); + const R rc = r * std_normal_tail_correction(r); + const R value = -(0.5 * z) * z - HALF_LOG_TWO_PI - log(a) + log1p(rc); + if constexpr (calc_grad) { + return {value, a / (1 + rc)}; + } else { + return {value, 0}; + } + } + // Not abs(z): its autodiff tangent is 0 at z == 0, which the slope needs. + const R x = (z < 0 ? R(-z) : z) * INV_SQRT_TWO; + if (x < 0.46875) { + const R e = std_normal_erf_small(x); + const R value = LOG_HALF + log1p(z < 0 ? R(-e) : e); + if constexpr (calc_grad) { + return {value, SQRT_TWO_OVER_SQRT_PI * exp(-square(x)) + / (z < 0 ? R(1 - e) : R(1 + e))}; + } else { + return {value, 0}; + } + } + const R erfcx = std_normal_erfcx(x); + if (z < 0) { + const R value = -(0.5 * z) * z + LOG_HALF + log(erfcx); + if constexpr (calc_grad) { + return {value, SQRT_TWO_OVER_SQRT_PI / erfcx}; + } else { + return {value, 0}; + } + } + const R density = exp(-(0.5 * z) * z); + const R tail = 0.5 * density * erfcx; + // Not log1m(tail): its domain check costs ~3 ns per element here. + const R value = log1p(-tail); + if constexpr (calc_grad) { + return {value, INV_SQRT_TWO_PI * density / (1 - tail)}; + } else { + return {value, 0}; + } +} + +/** Elementwise values and slopes; empty slopes when the gradient is off. */ +template * = nullptr> +inline auto std_normal_lcdf_value_grad(const T& z) { + using R = return_type_t>; + using Array = Eigen::Array; + const auto& z_ref = to_ref(z); + Array values(z_ref.size()); + Array slopes(calc_grad ? z_ref.size() : 0); + for (Eigen::Index i = 0; i < z_ref.size(); ++i) { + const auto result + = std_normal_lcdf_value_grad(R(z_ref.coeff(i))); + values[i] = result.first; + if constexpr (calc_grad) { + slopes[i] = result.second; + } + } + return std::make_pair(std::move(values), std::move(slopes)); +} + +} // namespace internal +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/meta/ref_type.hpp b/stan/math/prim/meta/ref_type.hpp index 538370c997b..0dc1da27295 100644 --- a/stan/math/prim/meta/ref_type.hpp +++ b/stan/math/prim/meta/ref_type.hpp @@ -48,8 +48,11 @@ struct ref_type_if< template struct ref_type_if> { - using type = - typename ref_type_if::Base>::type; + using T_base = typename std::decay_t::Base; + // Keep rvalues by value so a temporary is not referenced after it dies. + using type = typename ref_type_if< + Condition, std::conditional_t::value, + T_base&&, T_base>>::type; }; template diff --git a/stan/math/prim/prob/exp_mod_normal_ccdf_log.hpp b/stan/math/prim/prob/exp_mod_normal_ccdf_log.hpp index c5a84309029..b41a39b3825 100644 --- a/stan/math/prim/prob/exp_mod_normal_ccdf_log.hpp +++ b/stan/math/prim/prob/exp_mod_normal_ccdf_log.hpp @@ -14,8 +14,7 @@ template inline return_type_t exp_mod_normal_ccdf_log( const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda) { - return exp_mod_normal_lccdf(y, mu, sigma, - lambda); + return exp_mod_normal_lccdf(y, mu, sigma, lambda); } } // namespace math diff --git a/stan/math/prim/prob/exp_mod_normal_cdf.hpp b/stan/math/prim/prob/exp_mod_normal_cdf.hpp index baf9a5fd25f..8fcd392f7bb 100644 --- a/stan/math/prim/prob/exp_mod_normal_cdf.hpp +++ b/stan/math/prim/prob/exp_mod_normal_cdf.hpp @@ -2,23 +2,8 @@ #define STAN_MATH_PRIM_PROB_EXP_MOD_NORMAL_CDF_HPP #include -#include -#include -#include -#include -#include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { @@ -27,116 +12,10 @@ template * = nullptr> inline return_type_t exp_mod_normal_cdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, - const T_inv_scale& lambda) { - 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_sigma_ref = ref_type_if_not_constant_t; - using T_lambda_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "exp_mod_normal_cdf"; - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma, "Inv_scale parameter", - lambda); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - T_lambda_ref lambda_ref = lambda; - - 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); - decltype(auto) lambda_val - = to_ref(as_value_column_array_or_scalar(lambda_ref)); - - check_not_nan(function, "Random variable", y_val); - check_finite(function, "Location parameter", mu_val); - check_positive_finite(function, "Scale parameter", sigma_val); - check_positive_finite(function, "Inv_scale parameter", lambda_val); - - if (size_zero(y, mu, sigma, lambda)) { - return 1.0; - } - - auto ops_partials - = make_partials_propagator(y_ref, mu_ref, sigma_ref, lambda_ref); - - if constexpr (is_vector::value) { - if ((y_val == NEGATIVE_INFTY).any()) { - return ops_partials.build(0.0); - } - } else { - if (y_val == NEGATIVE_INFTY) { - return ops_partials.build(0.0); - } - } - - const auto& inv_sigma - = to_ref_if>(inv(sigma_val)); - const auto& diff = to_ref(y_val - mu_val); - const auto& v = to_ref(lambda_val * sigma_val); - const auto& scaled_diff = to_ref(diff * INV_SQRT_TWO * inv_sigma); - const auto& scaled_diff_diff - = to_ref_if>( - scaled_diff - v * INV_SQRT_TWO); - const auto& erf_calc = to_ref(0.5 * (1 + erf(scaled_diff_diff))); - - const auto& exp_term - = to_ref_if>( - exp(0.5 * square(v) - lambda_val * diff)); - const auto& cdf_n - = to_ref(0.5 + 0.5 * erf(scaled_diff) - exp_term * erf_calc); - - T_partials_return cdf(1.0); - if constexpr (is_vector::value) { - cdf = cdf_n.prod(); - } else { - cdf = cdf_n; - } - - if constexpr (is_any_autodiff_v) { - const auto& exp_term_2 = to_ref_if<( - is_any_autodiff_v && is_autodiff_v)>( - exp(-square(scaled_diff_diff))); - if constexpr (is_any_autodiff_v) { - constexpr bool need_deriv_refs - = is_any_autodiff_v && is_autodiff_v; - const auto& deriv_1 - = to_ref_if(lambda_val * exp_term * erf_calc); - const auto& deriv_2 = to_ref_if( - INV_SQRT_TWO_PI * exp_term * exp_term_2 * inv_sigma); - const auto& sq_scaled_diff = square(scaled_diff); - const auto& exp_m_sq_scaled_diff = exp(-sq_scaled_diff); - const auto& deriv_3 = to_ref_if( - INV_SQRT_TWO_PI * exp_m_sq_scaled_diff * inv_sigma); - if constexpr (is_any_autodiff_v) { - const auto& deriv - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - cdf * (deriv_1 - deriv_2 + deriv_3) / cdf_n); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = deriv; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = -deriv; - } - } - if constexpr (is_autodiff_v) { - edge<2>(ops_partials).partials_ - = -cdf - * ((deriv_1 - deriv_2) * v - + (deriv_3 - deriv_2) * scaled_diff * SQRT_TWO) - / cdf_n; - } - } - if constexpr (is_autodiff_v) { - edge<3>(ops_partials).partials_ - = cdf * exp_term - * (INV_SQRT_TWO_PI * sigma_val * exp_term_2 - - (v * sigma_val - diff) * erf_calc) - / cdf_n; - } - } - return ops_partials.build(cdf); + T_y&& y, T_loc&& mu, T_scale&& sigma, T_inv_scale&& lambda) { + return exp(internal::exp_mod_normal_lcdf_impl( + "exp_mod_normal_cdf", std::forward(y), std::forward(mu), + std::forward(sigma), std::forward(lambda))); } } // namespace math diff --git a/stan/math/prim/prob/exp_mod_normal_cdf_log.hpp b/stan/math/prim/prob/exp_mod_normal_cdf_log.hpp index 4908971ad37..1c59046db7e 100644 --- a/stan/math/prim/prob/exp_mod_normal_cdf_log.hpp +++ b/stan/math/prim/prob/exp_mod_normal_cdf_log.hpp @@ -14,8 +14,7 @@ template inline return_type_t exp_mod_normal_cdf_log( const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda) { - return exp_mod_normal_lcdf(y, mu, sigma, - lambda); + return exp_mod_normal_lcdf(y, mu, sigma, lambda); } } // namespace math diff --git a/stan/math/prim/prob/exp_mod_normal_lccdf.hpp b/stan/math/prim/prob/exp_mod_normal_lccdf.hpp index 2e5901704e0..28c542098fd 100644 --- a/stan/math/prim/prob/exp_mod_normal_lccdf.hpp +++ b/stan/math/prim/prob/exp_mod_normal_lccdf.hpp @@ -2,25 +2,7 @@ #define STAN_MATH_PRIM_PROB_EXP_MOD_NORMAL_LCCDF_HPP #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 { @@ -29,108 +11,10 @@ template * = nullptr> inline return_type_t exp_mod_normal_lccdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, - const T_inv_scale& lambda) { - 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_sigma_ref = ref_type_if_not_constant_t; - using T_lambda_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "exp_mod_normal_lccdf"; - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma, "Inv_scale parameter", - lambda); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - T_lambda_ref lambda_ref = lambda; - - 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); - decltype(auto) lambda_val - = to_ref(as_value_column_array_or_scalar(lambda_ref)); - - check_not_nan(function, "Random variable", y_val); - check_finite(function, "Location parameter", mu_val); - check_positive_finite(function, "Scale parameter", sigma_val); - check_positive_finite(function, "Inv_scale parameter", lambda_val); - - if (size_zero(y, mu, sigma, lambda)) { - return 0; - } - - auto ops_partials - = make_partials_propagator(y_ref, mu_ref, sigma_ref, lambda_ref); - - scalar_seq_view y_vec(y_val); - for (size_t n = 0, size_y = stan::math::size(y); n < size_y; n++) { - if (is_inf(y_vec[n])) { - return ops_partials.build(y_vec[n] > 0 ? negative_infinity() : 0); - } - } - - const auto& inv_sigma - = to_ref_if>(inv(sigma_val)); - const auto& diff = to_ref(y_val - mu_val); - const auto& v = to_ref(lambda_val * sigma_val); - const auto& scaled_diff = to_ref(diff * INV_SQRT_TWO * inv_sigma); - const auto& scaled_diff_diff - = to_ref_if>( - scaled_diff - v * INV_SQRT_TWO); - const auto& erf_calc = to_ref(0.5 * (1 + erf(scaled_diff_diff))); - - const auto& exp_term - = to_ref_if>( - exp(0.5 * square(v) - lambda_val * diff)); - const auto& ccdf_n - = to_ref(0.5 - 0.5 * erf(scaled_diff) + exp_term * erf_calc); - - T_partials_return ccdf_log = sum(log(ccdf_n)); - - if constexpr (is_any_autodiff_v) { - const auto& exp_term_2 = to_ref_if<( - is_any_autodiff_v && is_autodiff_v)>( - exp(-square(scaled_diff_diff))); - if constexpr (is_any_autodiff_v) { - constexpr bool need_deriv_refs - = is_any_autodiff_v && is_autodiff_v; - const auto& deriv_1 - = to_ref_if(lambda_val * exp_term * erf_calc); - const auto& deriv_2 = to_ref_if( - INV_SQRT_TWO_PI * exp_term * exp_term_2 * inv_sigma); - const auto& sq_scaled_diff = square(scaled_diff); - const auto& exp_m_sq_scaled_diff = exp(-sq_scaled_diff); - const auto& deriv_3 = to_ref_if( - INV_SQRT_TWO_PI * exp_m_sq_scaled_diff * inv_sigma); - if constexpr (is_any_autodiff_v) { - const auto& deriv - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - (deriv_1 - deriv_2 + deriv_3) / ccdf_n); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -deriv; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = deriv; - } - } - if constexpr (is_autodiff_v) { - edge<2>(ops_partials).partials_ - = ((deriv_1 - deriv_2) * v - + (deriv_3 - deriv_2) * scaled_diff * SQRT_TWO) - / ccdf_n; - } - } - if constexpr (is_autodiff_v) { - edge<3>(ops_partials).partials_ - = exp_term - * ((v * sigma_val - diff) * erf_calc - - INV_SQRT_TWO_PI * sigma_val * exp_term_2) - / ccdf_n; - } - } - - return ops_partials.build(ccdf_log); + T_y&& y, T_loc&& mu, T_scale&& sigma, T_inv_scale&& lambda) { + return internal::exp_mod_normal_lcdf_impl( + "exp_mod_normal_lccdf", std::forward(y), std::forward(mu), + std::forward(sigma), std::forward(lambda)); } } // namespace math diff --git a/stan/math/prim/prob/exp_mod_normal_lcdf.hpp b/stan/math/prim/prob/exp_mod_normal_lcdf.hpp index d5222619469..ac556cf4eda 100644 --- a/stan/math/prim/prob/exp_mod_normal_lcdf.hpp +++ b/stan/math/prim/prob/exp_mod_normal_lcdf.hpp @@ -3,135 +3,130 @@ #include #include -#include -#include +#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 { +namespace internal { -template * = nullptr> -inline return_type_t exp_mod_normal_lcdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, - const T_inv_scale& lambda) { +template +inline auto exp_mod_normal_log_combine(const T_a& a, const T_b& b) { + if constexpr (upper) { + return log_sum_exp(a, b); + } else { + return log_diff_exp(a, b); + } +} + +/** log F = log_diff_exp(a, b) with a = log Phi(z), b = v^2/2 - lambda (y - mu) + * + log Phi(z - v), v = lambda sigma; log (1 - F) = log_sum_exp(log Phi(-z), + * b). + */ +template +inline return_type_t exp_mod_normal_lcdf_impl( + const char* function, T_y&& y, T_loc&& mu, T_scale&& sigma, + T_inv_scale&& lambda) { 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_sigma_ref = ref_type_if_not_constant_t; using T_lambda_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "exp_mod_normal_lcdf"; + constexpr bool any_autodiff + = is_any_autodiff_v; + constexpr double sign = upper ? -1.0 : 1.0; check_consistent_sizes(function, "Random variable", y, "Location parameter", mu, "Scale parameter", sigma, "Inv_scale parameter", lambda); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - T_lambda_ref lambda_ref = lambda; - + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_sigma_ref sigma_ref = std::forward(sigma); + T_lambda_ref lambda_ref = std::forward(lambda); 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); decltype(auto) lambda_val = to_ref(as_value_column_array_or_scalar(lambda_ref)); - check_not_nan(function, "Random variable", y_val); check_finite(function, "Location parameter", mu_val); check_positive_finite(function, "Scale parameter", sigma_val); check_positive_finite(function, "Inv_scale parameter", lambda_val); - if (size_zero(y, mu, sigma, lambda)) { + if (size_zero(y_ref, mu_ref, sigma_ref, lambda_ref)) { return 0; } auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref, lambda_ref); - - scalar_seq_view y_vec(y_val); - for (size_t n = 0, size_y = stan::math::size(y); n < size_y; n++) { - if (is_inf(y_vec[n])) { - return ops_partials.build(y_vec[n] < 0 ? negative_infinity() : 0); - } + if (any(y_val == NEGATIVE_INFTY)) { + return ops_partials.build(upper ? 0.0 : NEGATIVE_INFTY); + } + if (any(y_val == INFTY)) { + return ops_partials.build(upper ? NEGATIVE_INFTY : 0.0); } - const auto& inv_sigma - = to_ref_if>(inv(sigma_val)); const auto& diff = to_ref(y_val - mu_val); + const auto& z = to_ref(diff / sigma_val); const auto& v = to_ref(lambda_val * sigma_val); - const auto& scaled_diff = to_ref(diff * INV_SQRT_TWO * inv_sigma); - const auto& scaled_diff_diff - = to_ref_if>( - scaled_diff - v * INV_SQRT_TWO); - const auto& erf_calc = to_ref(0.5 * (1 + erf(scaled_diff_diff))); + const auto [log_a, slope_a] + = internal::std_normal_lcdf_value_grad(sign * z); + const auto [log_phi_b, slope_b] + = internal::std_normal_lcdf_value_grad(z - v); + const auto& log_b = to_ref_if(0.5 * square(v) + - lambda_val * diff + log_phi_b); + const auto& lp = to_ref_if( + exp_mod_normal_log_combine(log_a, log_b)); + const T_partials_return cdf_log = sum(lp); - const auto& exp_term - = to_ref_if>( - exp(0.5 * square(v) - lambda_val * diff)); - const auto& cdf_n - = to_ref(0.5 + 0.5 * erf(scaled_diff) - exp_term * erf_calc); - - T_partials_return cdf_log = sum(log(cdf_n)); - - if constexpr (is_any_autodiff_v) { - const auto& exp_term_2 = to_ref_if<( - is_any_autodiff_v && is_autodiff_v)>( - exp(-square(scaled_diff_diff))); - if constexpr (is_any_autodiff_v) { - constexpr bool need_deriv_refs - = is_any_autodiff_v && is_autodiff_v; - const auto& deriv_1 - = to_ref_if(lambda_val * exp_term * erf_calc); - const auto& deriv_2 = to_ref_if( - INV_SQRT_TWO_PI * exp_term * exp_term_2 * inv_sigma); - const auto& sq_scaled_diff = square(scaled_diff); - const auto& exp_m_sq_scaled_diff = exp(-sq_scaled_diff); - const auto& deriv_3 = to_ref_if( - INV_SQRT_TWO_PI * exp_m_sq_scaled_diff * inv_sigma); - if constexpr (is_any_autodiff_v) { - const auto& deriv - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - (deriv_1 - deriv_2 + deriv_3) / cdf_n); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = deriv; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = -deriv; - } - } - if constexpr (is_autodiff_v) { - edge<2>(ops_partials).partials_ - = -((deriv_1 - deriv_2) * v - + (deriv_3 - deriv_2) * scaled_diff * SQRT_TWO) - / cdf_n; - } + if constexpr (any_autodiff) { + // Weights of the two terms in the total; the second is signed. + const auto& w_b = to_ref(-sign * exp(log_b - lp)); + const auto& s_b = to_ref(w_b * slope_b); + const auto& s = to_ref_if< + (is_autodiff_v + is_autodiff_v + is_autodiff_v) + >= 2>((sign * exp(log_a - lp) * slope_a + s_b) / sigma_val); + const auto& q + = to_ref_if>(w_b * v - s_b); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = s - w_b * lambda_val; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = w_b * lambda_val - s; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) = select(s == 0, 0.0, -s * z) + lambda_val * q; } if constexpr (is_autodiff_v) { - edge<3>(ops_partials).partials_ - = exp_term - * (INV_SQRT_TWO_PI * sigma_val * exp_term_2 - - (v * sigma_val - diff) * erf_calc) - / cdf_n; + partials<3>(ops_partials) = sigma_val * q - w_b * diff; } } return ops_partials.build(cdf_log); } +} // namespace internal + +template * = nullptr> +inline return_type_t exp_mod_normal_lcdf( + T_y&& y, T_loc&& mu, T_scale&& sigma, T_inv_scale&& lambda) { + return internal::exp_mod_normal_lcdf_impl( + "exp_mod_normal_lcdf", std::forward(y), std::forward(mu), + std::forward(sigma), std::forward(lambda)); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp index cb3374ddc1b..e9002c65cd0 100644 --- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp +++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp @@ -7,7 +7,6 @@ #include #include #include -#include #include #include #include @@ -17,6 +16,7 @@ #include #include #include +#include #include namespace stan { @@ -27,8 +27,7 @@ template * = nullptr> inline return_type_t exp_mod_normal_lpdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, - const T_inv_scale& lambda) { + T_y&& y, T_loc&& mu, T_scale&& sigma, T_inv_scale&& lambda) { 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; @@ -38,10 +37,10 @@ inline return_type_t exp_mod_normal_lpdf( check_consistent_sizes(function, "Random variable", y, "Location parameter", mu, "Scale parameter", sigma, "Inv_scale parameter", lambda); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - T_lambda_ref lambda_ref = lambda; + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_sigma_ref sigma_ref = std::forward(sigma); + T_lambda_ref lambda_ref = std::forward(lambda); 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)); @@ -54,7 +53,7 @@ inline return_type_t exp_mod_normal_lpdf( check_positive_finite(function, "Scale parameter", sigma_val); check_positive_finite(function, "Inv_scale parameter", lambda_val); - if (size_zero(y, mu, sigma, lambda)) { + if (size_zero(y_ref, mu_ref, sigma_ref, lambda_ref)) { return 0.0; } if constexpr (!include_summand exp_mod_normal_lpdf( const auto& sigma_sq = to_ref_if>(square(sigma_val)); const auto& lambda_sigma_sq = to_ref(lambda_val * sigma_sq); const auto& mu_minus_y = to_ref(mu_val - y_val); - const auto& inner_term - = to_ref_if>( - (mu_minus_y + lambda_sigma_sq) * INV_SQRT_TWO * inv_sigma); - const auto& erfc_calc = to_ref(erfc(inner_term)); + // log(erfc(t) / 2) = log Phi(-sqrt(2) t) cancels the log(1/2) constant. + const auto& z = to_ref(-(mu_minus_y + lambda_sigma_sq) * inv_sigma); + const auto [values, slopes] = internal::std_normal_lcdf_value_grad< + is_any_autodiff_v>(z); - size_t N = max_size(y, mu, sigma, lambda); + size_t N = max_size(y_ref, mu_ref, sigma_ref, lambda_ref); T_partials_return logp(0.0); - if constexpr (include_summand::value) { - logp -= LOG_TWO * N; - } if constexpr (include_summand::value) { - logp += sum(log(lambda_val)) * N / math::size(lambda); + logp += sum(log(lambda_val)) * N / math::size(lambda_ref); } - const auto& log_erfc_calc = log(erfc_calc); - logp - += sum(lambda_val * (mu_minus_y + 0.5 * lambda_sigma_sq) + log_erfc_calc); + logp += sum(lambda_val * (mu_minus_y + 0.5 * lambda_sigma_sq) + values); auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref, lambda_ref); if constexpr (is_any_autodiff_v) { - const auto& exp_m_sq_inner_term = exp(-square(inner_term)); - const auto& deriv_logerfc = to_ref_if< - is_any_autodiff_v< - T_y, - T_loc> + is_autodiff_v + is_autodiff_v >= 2>( - -SQRT_TWO_OVER_SQRT_PI * exp_m_sq_inner_term / erfc_calc); if constexpr (is_any_autodiff_v) { const auto& deriv = to_ref_if>( - lambda_val + deriv_logerfc * inv_sigma); + lambda_val - slopes * inv_sigma); if constexpr (is_autodiff_v) { partials<0>(ops_partials) = -deriv; } @@ -107,11 +95,11 @@ inline return_type_t exp_mod_normal_lpdf( if constexpr (is_autodiff_v) { edge<2>(ops_partials).partials_ = sigma_val * square(lambda_val) - + deriv_logerfc * (lambda_val - mu_minus_y / sigma_sq); + - slopes * (lambda_val - mu_minus_y * square(inv_sigma)); } if constexpr (is_autodiff_v) { - partials<3>(ops_partials) = inv(lambda_val) + lambda_sigma_sq + mu_minus_y - + deriv_logerfc * sigma_val; + partials<3>(ops_partials) + = inv(lambda_val) + lambda_sigma_sq + mu_minus_y - slopes * sigma_val; } } @@ -120,9 +108,10 @@ inline return_type_t exp_mod_normal_lpdf( template inline return_type_t exp_mod_normal_lpdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, - const T_inv_scale& lambda) { - return exp_mod_normal_lpdf(y, mu, sigma, lambda); + T_y&& y, T_loc&& mu, T_scale&& sigma, T_inv_scale&& lambda) { + return exp_mod_normal_lpdf( + std::forward(y), std::forward(mu), + std::forward(sigma), std::forward(lambda)); } } // namespace math diff --git a/stan/math/prim/prob/lognormal_ccdf_log.hpp b/stan/math/prim/prob/lognormal_ccdf_log.hpp index 04288340b21..0c8e0cf6ec8 100644 --- a/stan/math/prim/prob/lognormal_ccdf_log.hpp +++ b/stan/math/prim/prob/lognormal_ccdf_log.hpp @@ -13,7 +13,7 @@ namespace math { template inline return_type_t lognormal_ccdf_log( const T_y& y, const T_loc& mu, const T_scale& sigma) { - return lognormal_lccdf(y, mu, sigma); + return lognormal_lccdf(y, mu, sigma); } } // namespace math diff --git a/stan/math/prim/prob/lognormal_cdf.hpp b/stan/math/prim/prob/lognormal_cdf.hpp index f739cc98acb..fba45f3254d 100644 --- a/stan/math/prim/prob/lognormal_cdf.hpp +++ b/stan/math/prim/prob/lognormal_cdf.hpp @@ -2,22 +2,8 @@ #define STAN_MATH_PRIM_PROB_LOGNORMAL_CDF_HPP #include -#include -#include -#include -#include -#include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { @@ -25,63 +11,11 @@ namespace math { template * = nullptr> -inline return_type_t lognormal_cdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { - 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_sigma_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "lognormal_cdf"; - - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - - 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); - - check_nonnegative(function, "Random variable", y_val); - check_finite(function, "Location parameter", mu_val); - check_positive_finite(function, "Scale parameter", sigma_val); - - if (size_zero(y, mu, sigma)) { - return 1.0; - } - - auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); - - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(0.0); - } - - const auto& log_y = log(y_val); - const auto& scaled_diff = to_ref_if>( - (log_y - mu_val) / (sigma_val * SQRT_TWO)); - const auto& erfc_m_diff = erfc(-scaled_diff); - const auto& cdf_n - = to_ref_if>(0.5 * erfc_m_diff); - - T_partials_return cdf = prod(cdf_n); - - if constexpr (is_any_autodiff_v) { - const auto& exp_m_sq_diff = exp(-scaled_diff * scaled_diff); - const auto& rep_deriv = to_ref_if< - is_autodiff_v< - T_y> + is_autodiff_v + is_autodiff_v >= 2>( - -cdf * INV_SQRT_TWO_PI * exp_m_sq_diff / (sigma_val * cdf_n)); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -rep_deriv / y_val; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = rep_deriv; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = rep_deriv * scaled_diff * SQRT_TWO; - } - } - return ops_partials.build(cdf); +inline return_type_t lognormal_cdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return exp(internal::lognormal_lcdf_impl( + "lognormal_cdf", std::forward(y), std::forward(mu), + std::forward(sigma))); } } // namespace math diff --git a/stan/math/prim/prob/lognormal_cdf_log.hpp b/stan/math/prim/prob/lognormal_cdf_log.hpp index effbc38c8ec..4eff99e367c 100644 --- a/stan/math/prim/prob/lognormal_cdf_log.hpp +++ b/stan/math/prim/prob/lognormal_cdf_log.hpp @@ -13,7 +13,7 @@ namespace math { template inline return_type_t lognormal_cdf_log( const T_y& y, const T_loc& mu, const T_scale& sigma) { - return lognormal_lcdf(y, mu, sigma); + return lognormal_lcdf(y, mu, sigma); } } // namespace math diff --git a/stan/math/prim/prob/lognormal_lccdf.hpp b/stan/math/prim/prob/lognormal_lccdf.hpp index 39391bedb0f..4b10aa06d12 100644 --- a/stan/math/prim/prob/lognormal_lccdf.hpp +++ b/stan/math/prim/prob/lognormal_lccdf.hpp @@ -2,22 +2,7 @@ #define STAN_MATH_PRIM_PROB_LOGNORMAL_LCCDF_HPP #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { @@ -25,62 +10,11 @@ namespace math { template * = nullptr> -inline return_type_t lognormal_lccdf( - const T_y& y, const T_loc& mu, const T_scale& sigma) { - 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_sigma_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "lognormal_lccdf"; - - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - - 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); - - check_nonnegative(function, "Random variable", y_val); - check_finite(function, "Location parameter", mu_val); - check_positive_finite(function, "Scale parameter", sigma_val); - - if (size_zero(y, mu, sigma)) { - return 0; - } - - auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); - - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(0.0); - } - - const auto& log_y = log(y_val); - const auto& scaled_diff = to_ref_if>( - (log_y - mu_val) / (sigma_val * SQRT_TWO)); - const auto& erfc_calc - = to_ref_if>(erfc(scaled_diff)); - - size_t N = max_size(y, mu, sigma); - T_partials_return ccdf_log = N * LOG_HALF + sum(log(erfc_calc)); - - if constexpr (is_any_autodiff_v) { - const auto& exp_m_sq_diff = exp(-scaled_diff * scaled_diff); - const auto& rep_deriv = to_ref_if< - is_autodiff_v< - T_y> + is_autodiff_v + is_autodiff_v >= 2>( - SQRT_TWO_OVER_SQRT_PI * exp_m_sq_diff / (sigma_val * erfc_calc)); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -rep_deriv / y_val; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = rep_deriv; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = rep_deriv * scaled_diff * SQRT_TWO; - } - } - return ops_partials.build(ccdf_log); +inline return_type_t lognormal_lccdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return internal::lognormal_lcdf_impl( + "lognormal_lccdf", std::forward(y), std::forward(mu), + std::forward(sigma)); } } // namespace math diff --git a/stan/math/prim/prob/lognormal_lcdf.hpp b/stan/math/prim/prob/lognormal_lcdf.hpp index 482f79885bd..8223f5f8472 100644 --- a/stan/math/prim/prob/lognormal_lcdf.hpp +++ b/stan/math/prim/prob/lognormal_lcdf.hpp @@ -3,86 +3,84 @@ #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 { +namespace internal { -template * = nullptr> -inline return_type_t lognormal_lcdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { +template +inline return_type_t lognormal_lcdf_impl( + const char* function, T_y&& y, T_loc&& mu, T_scale&& sigma) { 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_sigma_ref = ref_type_if_not_constant_t; - static constexpr const char* function = "lognormal_lcdf"; - - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - + constexpr double sign = reflect ? -1.0 : 1.0; + check_consistent_sizes(function, "Random variable", y, "Location parameter", + mu, "Scale parameter", sigma); + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_sigma_ref sigma_ref = std::forward(sigma); 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); - check_nonnegative(function, "Random variable", y_val); check_finite(function, "Location parameter", mu_val); check_positive_finite(function, "Scale parameter", sigma_val); - if (size_zero(y, mu, sigma)) { + if (size_zero(y_ref, mu_ref, sigma_ref)) { return 0; } auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); - - if (sum(promote_scalar(y_val == 0))) { - return ops_partials.build(NEGATIVE_INFTY); + if (any(y_val == 0)) { + return ops_partials.build(reflect ? 0.0 : NEGATIVE_INFTY); } - const auto& log_y = log(y_val); - const auto& scaled_diff = to_ref_if>( - (log_y - mu_val) / (sigma_val * SQRT_TWO)); - const auto& erfc_calc - = to_ref_if>(erfc(-scaled_diff)); - size_t N = max_size(y, mu, sigma); - T_partials_return cdf_log = N * LOG_HALF + sum(log(erfc_calc)); - + const auto& z = to_ref(sign * (log(y_val) - mu_val) / sigma_val); + const auto [values, slopes] = internal::std_normal_lcdf_value_grad< + is_any_autodiff_v>(z); + const T_partials_return cdf_log = sum(values); if constexpr (is_any_autodiff_v) { - const auto& exp_m_sq_diff = exp(-scaled_diff * scaled_diff); - const auto& rep_deriv = to_ref_if< - is_autodiff_v< - T_y> + is_autodiff_v + is_autodiff_v >= 2>( - -SQRT_TWO_OVER_SQRT_PI * exp_m_sq_diff / (sigma_val * erfc_calc)); + const auto& scaled_slope = to_ref_if< + (is_autodiff_v + is_autodiff_v + is_autodiff_v) + >= 2>(slopes / sigma_val); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -rep_deriv / y_val; + partials<0>(ops_partials) = sign * scaled_slope / y_val; } if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = rep_deriv; + partials<1>(ops_partials) = -sign * scaled_slope; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = rep_deriv * scaled_diff * SQRT_TWO; + partials<2>(ops_partials) + = select(scaled_slope == 0, 0.0, -scaled_slope * z); } } return ops_partials.build(cdf_log); } +} // namespace internal + +template * = nullptr> +inline return_type_t lognormal_lcdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return internal::lognormal_lcdf_impl( + "lognormal_lcdf", std::forward(y), std::forward(mu), + std::forward(sigma)); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/prim/prob/normal_ccdf_log.hpp b/stan/math/prim/prob/normal_ccdf_log.hpp index d7df99b01dd..72ce810fec8 100644 --- a/stan/math/prim/prob/normal_ccdf_log.hpp +++ b/stan/math/prim/prob/normal_ccdf_log.hpp @@ -13,7 +13,7 @@ namespace math { template inline return_type_t normal_ccdf_log( const T_y& y, const T_loc& mu, const T_scale& sigma) { - return normal_lccdf(y, mu, sigma); + return normal_lccdf(y, mu, sigma); } } // namespace math diff --git a/stan/math/prim/prob/normal_cdf.hpp b/stan/math/prim/prob/normal_cdf.hpp index 5cb9bce6e0d..e115bf3c7ce 100644 --- a/stan/math/prim/prob/normal_cdf.hpp +++ b/stan/math/prim/prob/normal_cdf.hpp @@ -2,18 +2,8 @@ #define STAN_MATH_PRIM_PROB_NORMAL_CDF_HPP #include -#include -#include -#include -#include #include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { @@ -24,6 +14,9 @@ namespace math { * * \f$\Phi(x) = \frac{1}{\sqrt{2 \pi}} \int_{-\inf}^x e^{-t^2/2} dt\f$. * + * Evaluated as the exponential of the log cdf so the tails and gradients + * share the log kernel instead of hard cutoffs. + * * @tparam T_y type of y * @tparam T_loc type of mean parameter * @tparam T_scale type of standard deviation parameter @@ -35,89 +28,11 @@ namespace math { template * = nullptr> -inline return_type_t normal_cdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { - using T_partials_return = partials_return_t; - using std::exp; - using T_y_ref = ref_type_t; - using T_mu_ref = ref_type_t; - using T_sigma_ref = ref_type_t; - static constexpr const char* function = "normal_cdf"; - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - check_not_nan(function, "Random variable", y_ref); - check_finite(function, "Location parameter", mu_ref); - check_positive(function, "Scale parameter", sigma_ref); - - if (size_zero(y, mu, sigma)) { - return 1.0; - } - - T_partials_return cdf(1.0); - auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); - - scalar_seq_view y_vec(y_ref); - scalar_seq_view mu_vec(mu_ref); - scalar_seq_view sigma_vec(sigma_ref); - size_t N = max_size(y, mu, sigma); - - 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); - const T_partials_return sigma_dbl = sigma_vec.val(n); - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) / (sigma_dbl * SQRT_TWO); - T_partials_return cdf_n; - if (scaled_diff < -37.5 * INV_SQRT_TWO) { - cdf_n = 0.0; - } else if (scaled_diff < -5.0 * INV_SQRT_TWO) { - cdf_n = 0.5 * erfc(-scaled_diff); - } else if (scaled_diff > 8.25 * INV_SQRT_TWO) { - cdf_n = 1; - } else { - cdf_n = 0.5 * (1.0 + erf(scaled_diff)); - } - - cdf *= cdf_n; - - if constexpr (is_any_autodiff_v) { - const T_partials_return rep_deriv - = (scaled_diff < -37.5 * INV_SQRT_TWO) - ? 0.0 - : INV_SQRT_TWO_PI * exp(-scaled_diff * scaled_diff) - / (cdf_n * sigma_dbl); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += rep_deriv; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] -= rep_deriv; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= rep_deriv * scaled_diff * SQRT_TWO; - } - } - } - - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < stan::math::size(y); ++n) { - partials<0>(ops_partials)[n] *= cdf; - } - } - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < stan::math::size(mu); ++n) { - partials<1>(ops_partials)[n] *= cdf; - } - } - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < stan::math::size(sigma); ++n) { - partials<2>(ops_partials)[n] *= cdf; - } - } - return ops_partials.build(cdf); +inline return_type_t normal_cdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return exp(internal::normal_lcdf_impl( + "normal_cdf", std::forward(y), std::forward(mu), + std::forward(sigma))); } } // namespace math diff --git a/stan/math/prim/prob/normal_lccdf.hpp b/stan/math/prim/prob/normal_lccdf.hpp index cbdd074418c..b55cff0c53a 100644 --- a/stan/math/prim/prob/normal_lccdf.hpp +++ b/stan/math/prim/prob/normal_lccdf.hpp @@ -5,18 +5,15 @@ namespace stan { namespace math { -namespace internal { -constexpr char normal_lccdf_func[] = "normal_lccdf"; -} // namespace internal template * = nullptr> -inline return_type_t normal_lccdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { - return normal_lcdf( - -as_array_or_scalar(y), -as_array_or_scalar(mu), sigma); +inline return_type_t normal_lccdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return internal::normal_lcdf_impl("normal_lccdf", std::forward(y), + std::forward(mu), + std::forward(sigma)); } } // namespace math diff --git a/stan/math/prim/prob/normal_lcdf.hpp b/stan/math/prim/prob/normal_lcdf.hpp index b7446c9fd82..5fcfe214ab1 100644 --- a/stan/math/prim/prob/normal_lcdf.hpp +++ b/stan/math/prim/prob/normal_lcdf.hpp @@ -3,130 +3,78 @@ #include #include -#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 { namespace internal { -constexpr char normal_lcdf_func[] = "normal_lcdf"; + +/** Log of the normal cdf, or of its complement when `reflect` is set: the + * standardized value is negated so no negated autodiff operands are built. + */ +template +inline return_type_t normal_lcdf_impl(const char* function, + T_y&& y, T_loc&& mu, + T_scale&& sigma) { + 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_sigma_ref = ref_type_if_not_constant_t; + constexpr double sign = reflect ? -1.0 : 1.0; + check_consistent_sizes(function, "Random variable", y, "Location parameter", + mu, "Scale parameter", sigma); + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_sigma_ref sigma_ref = std::forward(sigma); + 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) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref)); + check_not_nan(function, "Random variable", y_val); + check_finite(function, "Location parameter", mu_val); + check_positive(function, "Scale parameter", sigma_val); + + if (size_zero(y_ref, mu_ref, sigma_ref)) { + return 0; + } + + const auto& z = to_ref(sign * (y_val - mu_val) / sigma_val); + const auto [values, slopes] = internal::std_normal_lcdf_value_grad< + is_any_autodiff_v>(z); + const T_partials_return cdf_log = sum(values); + auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); + if constexpr (is_any_autodiff_v) { + const auto& scaled_slope = to_ref_if< + (is_autodiff_v + is_autodiff_v + is_autodiff_v) + >= 2>(slopes / sigma_val); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = sign * scaled_slope; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = -sign * scaled_slope; + } + if constexpr (is_autodiff_v) { + // The positive infinite endpoint has a zero slope, not 0 * infinity. + partials<2>(ops_partials) + = select(scaled_slope == 0, 0.0, -scaled_slope * z); + } + } + return ops_partials.build(cdf_log); +} + } // namespace internal /** \ingroup prob_dists * @brief Calculates the log of the cdf of the normal distribution * - * Tail branching follows three published results, and matches what other - * libraries do: - * - Abramowitz & Stegun (1964) 7.1.26 Page 299 gives - * `erf(x) = 1 - P(t) exp(-x^2)`, `t = 1/(1 + p x)`, with `exp(-x^2)` as a - * numerator factor: - * https://archive.org/details/handbookofmathem1964abra/page/298/mode/2up - * - Cody (1969), Math. Comp. 23:631-637, is the rational approximation used - * for the tail value: https://doi.org/10.1090/S0025-5718-1969-0247736-4 - * - Digital Library of Mathematical Functions (DLMF) 7.12.1 - * is the erfc asymptotic behind the far-tail Mills ratio: - * https://dlmf.nist.gov/7.12.E1 - * - * R's `pnorm` only ever forms the Gaussian factor in the numerator - * (the `do_del` macro), and switches to the Cody tail form at `y > M_SQRT_32`, - * i.e. `|x| > sqrt(32) ~= 5.657`. Since `scaled_diff = x / sqrt(2)`, that same - * crossover is `|scaled_diff| > sqrt(32)/sqrt(2) = 4` exactly, which is the - * cutoff used here for the cdf value. Measured by evaluating the - * `temp_p`/`temp_q` expression below against `log(erfc(-scaled_diff)/2)` at 60 - * significant digits, our Cody set holds to 8.6e-17 relative down to - * scaled_diff = -4 and degrades past about -3.52, so 4 sits just inside its - * range. This is the same as R's impl. - * https://github.com/wch/r-source/blob/trunk/src/nmath/pnorm.c - * SciPy's `log_ndtr` uses the identical `log1p(-erfc(t)/2)` upper branch with - * the same `t = x/sqrt(2)`, and needs no rational approximation or Taylor - * patches at all below `x = -1`, because `erfcx` never forms `exp(+t^2)`: - * https://github.com/scipy/xsf/blob/main/include/xsf/stats.h - * - * The interior cutoffs for the gradient are not from the literature. - * They were derived in stan-dev/math#1411 - * (Phil Clemson, Univ. of Liverpool, Nov 2019), fixing - * stan-dev/math#1284, which describes them as "original Taylor expansions - * that have been derived to bridge the gap where the numerical - * approximations are unstable". The author's account of how they were - * placed, from the review thread, was: "After playing around with the - * autodiff tester I found some regions where it was failing some of the - * tests, so I added some new approximations (Taylor expansions and fits of - * the residuals) to improve the accuracy." - * https://github.com/stan-dev/math/pull/1411 - * https://github.com/stan-dev/math/issues/1284 - * - * Each row of the chart below is one Taylor expansion around `scaled_diff`. - * `interval` is the range of `scaled_diff` it covers. - * `centre` is the point the series is expanded about. - * - * | interval | centre | half-width | - * |---|---|---| - * | (2.5, 2.9] | 2.7 | 0.200 | - * | (2.1, 2.5] | 2.3 | 0.200 | - * | (1.5, 2.1] | 1.85 | 0.300 | - * | (0.8, 1.5] | 1.15 | 0.350 | - * | (0.1, 0.8] | 0.45 | 0.350 | - * - * The cut points are not arbitrary. Each interval of `scaled_diff` is centred - * on its own Taylor point, which minimises the largest `|t|` the series has to - * cover. - * - * The second table below shows why each interval ends where it does: the - * series is accurate across its own range and falls apart just outside it, - * so the intervals cannot be widened to use fewer of them. - * - * Worst in-range comes from evaluating each Taylor polynomial as written - * against `(2/sqrt(pi)) * exp(-scaled_diff^2) / erfc(-scaled_diff)`, the exact - * derivative, at 60 significant digits. + * Uses the shared standard-normal kernel after standardization. See + * std_normal_lcdf_impl.hpp for the Cody approximations and their crossovers. * - * Columns below: - * `worst in-range`: The largest relative error of that branch's series over - * its own interval; - * `0.2 below lo`: The relative error the same series would give at `lo - 0.2`; - * `0.2 above hi`: The relative error the same series would give at `hi + 0.2`, - * i.e. what widening the interval either way would cost. - * - * | interval | worst in-range | 0.2 below lo | 0.2 above hi | - * |---|---|---|---| - * | (2.5, 2.9] | 3.27e-05 | 2.32e-05 | 1.61e-02 | - * | (2.1, 2.5] | 2.80e-05 | 3.25e-04 | 9.15e-03 | - * | (1.5, 2.1] | 2.30e-05 | 8.11e-05 | 3.77e-03 | - * | (0.8, 1.5] | 6.09e-05 | 9.29e-05 | 2.93e-03 | - * | (0.1, 0.8] | 7.61e-06 | 3.62e-05 | 2.82e-04 | - * - * The worst in-range error is uniformly 1e-5 to 6e-5, just inside the 1e-4 - * relative gradient tolerance `expect_ad` applies by default - * (`gradient_grad_` in test/unit/math/ad_tolerances.hpp; the 1e-8 there is - * `gradient_val_`, which bounds the value rather than the gradient), while - * widening an interval by 0.2 costs one to three orders of magnitude. That - * uniformity, not any published result, is the placement criterion. - * - * The `worst in-range` column is enforced: the branch_accuracy test in - * mix/prob/normal_cdf_log_test.cpp asserts every branch against - * 60-digit references at that error plus a small margin. - * - * The negative-tail residual corrections at `scaled_diff` of -2.1, -3.9, -7 and - * -17 are cubic fits of residuals from the same PR and likewise have no - * external source. The -29 cutoff below is justified by DLMF 7.12.1 - * - * @tparam func name reported by the error checks. Reflected distributions - * such as `normal_lccdf` delegate here and pass their own name so that - * exceptions name the function the user actually called. * @tparam T_y A vector or scalar type for the random variable. * @tparam T_loc A vector or scalar type for the location parameter. * @tparam T_scale A vector or scalar type for the scale parameter. @@ -136,216 +84,14 @@ constexpr char normal_lcdf_func[] = "normal_lcdf"; * @return The log of the normal cdf evaluated at the specified arguments. If * given containers, the log of the product of the cdfs. */ -template * = nullptr> -inline return_type_t normal_lcdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { - using T_partials_return = partials_return_t; - using T_y_ref = ref_type_t; - using T_mu_ref = ref_type_t; - using T_sigma_ref = ref_type_t; - static constexpr const char* function = func; - check_consistent_sizes(function, "Random variable", y, "Location parameter", - mu, "Scale parameter", sigma); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - check_not_nan(function, "Random variable", y_ref); - check_finite(function, "Location parameter", mu_ref); - check_positive(function, "Scale parameter", sigma_ref); - - if (size_zero(y, mu, sigma)) { - return 0; - } - - T_partials_return cdf_log(0.0); - auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref); - - scalar_seq_view y_vec(y_ref); - scalar_seq_view mu_vec(mu_ref); - scalar_seq_view sigma_vec(sigma_ref); - size_t N = max_size(y, mu, sigma); - - 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); - const T_partials_return sigma_dbl = sigma_vec.val(n); - - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) / (sigma_dbl * SQRT_TWO); - - const T_partials_return x2 = square(scaled_diff); - - // Rigorous numerical approximations are applied here to deal with values - // of |scaled_diff|>>0. This is needed to deal with rare base-rate - // logistic regression problems where it is useful to use an alternative - // link function instead. - // - // use erfc() instead of erf() in order to retain precision - // since for x>0 erfc()->0 - if (scaled_diff > 0.0) { - // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) - cdf_log += log1p(-0.5 * erfc(scaled_diff)); - if (!is_not_nan(cdf_log)) { - cdf_log = 0; - } - } else if (scaled_diff > -4.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x); -4 is R pnorm's M_SQRT_32 - // Since we scale by sqrt(2), we use sqrt(32)/sqrt(2) = 4 - cdf_log += log(erfc(-scaled_diff)) + LOG_HALF; - } else if (10.0 * log(fabs(scaled_diff)) - < log(std::numeric_limits::max())) { - // entering territory where erfc(-x)~0 - // need to use direct numerical approximation of cdf_log instead - // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) - // CDF(x) = 1/2erfc(-x) - const T_partials_return x4 = pow(scaled_diff, 4); - const T_partials_return x6 = pow(scaled_diff, 6); - const T_partials_return x8 = pow(scaled_diff, 8); - const T_partials_return x10 = pow(scaled_diff, 10); - const T_partials_return temp_p - = 0.000658749161529837803157 + 0.0160837851487422766278 / x2 - + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6 - + 0.305326634961232344035 / x8 + 0.0163153871373020978498 / x10; - const T_partials_return temp_q - = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 - - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6 - - 2.56852019228982242072 / x8 - 1.0 / x10; - cdf_log += LOG_HALF + log(INV_SQRT_PI + (temp_p / temp_q) / x2) - - log(-scaled_diff) - x2; - } else { - // scaled_diff^10 term will overflow - cdf_log = stan::math::negative_infinity(); - } - - if constexpr (is_any_autodiff_v) { - // compute partial derivatives - // based on analytic form given by: - // dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2+erf(x)/2) - T_partials_return dncdf_log = 0.0; - T_partials_return t = 0.0; - T_partials_return t2 = 0.0; - T_partials_return t4 = 0.0; - - // calculate using piecewise function - // (due to instability / inaccuracy in the various approximations) - if (scaled_diff > 2.9) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - t = 1.0 / (1.0 + 0.3275911 * scaled_diff); - t2 = square(t); - t4 = pow(t, 4); - // A&S 7.1.26 puts exp(-x2) in the numerator; keep it there so it - // underflows to zero instead of overflowing inside a denominator - const T_partials_return exp_m_x2 = exp(-x2); - dncdf_log - = exp_m_x2 - / (SQRT_PI - * (1.0 - - exp_m_x2 - * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 - - 1.453152027 * t2 * t + 1.061405429 * t4))); - } else if (scaled_diff > 2.5) { - // in the trouble area where all of the standard numerical - // approximations are unstable - bridge the gap using Taylor - // expansions of the analytic function - // use Taylor expansion centred around x=2.7 - t = scaled_diff - 2.7; - t2 = square(t); - t4 = pow(t, 4); - dncdf_log = 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2 - - 0.008029540137 * t2 * t + 0.008232190507 * t4 - - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6); - } else if (scaled_diff > 2.1) { - // use Taylor expansion centred around x=2.3 - t = scaled_diff - 2.3; - t2 = square(t); - t4 = pow(t, 4); - dncdf_log = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 - - 0.03326906904 * t2 * t + 0.02482478940 * t4 - - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); - } else if (scaled_diff > 1.5) { - // use Taylor expansion centred around x=1.85 - t = scaled_diff - 1.85; - t2 = square(t); - t4 = pow(t, 4); - dncdf_log = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 - - 0.09274533184 * t2 * t + 0.03543327418 * t4 - + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); - } else if (scaled_diff > 0.8) { - // use Taylor expansion centred around x=1.15 - t = scaled_diff - 1.15; - t2 = square(t); - t4 = pow(t, 4); - dncdf_log = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 - - 0.09748053605 * t2 * t - 0.04347986191 * t4 - + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); - } else if (scaled_diff > 0.1) { - // use Taylor expansion centred around x=0.45 - t = scaled_diff - 0.45; - t2 = square(t); - t4 = pow(t, 4); - dncdf_log = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 - + 0.04700850676 * t2 * t - 0.03478651979 * t4 - - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); - } else if (scaled_diff < -29.0) { - // asymptotic Mills ratio, DLMF 7.12.1: dncdf_log grows linearly as - // -2*scaled_diff, so no quadratic residual fit can track it - const T_partials_return inv_x2 = 1.0 / x2; - dncdf_log - = -2.0 * scaled_diff - / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875))); - } else if (10.0 * log(fabs(scaled_diff)) - < log(std::numeric_limits::max())) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - // use fact that erf(x)=-erf(-x) - // Abramowitz and Stegun define this for -inf) { - partials<0>(ops_partials)[n] += dncdf_log / sigma_sqrt2; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] -= dncdf_log / sigma_sqrt2; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= dncdf_log * scaled_diff / sigma_dbl; - } - } - } - return ops_partials.build(cdf_log); +inline return_type_t normal_lcdf(T_y&& y, T_loc&& mu, + T_scale&& sigma) { + return internal::normal_lcdf_impl("normal_lcdf", std::forward(y), + std::forward(mu), + std::forward(sigma)); } } // namespace math diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp index 50ed5217858..a1cede591f1 100644 --- a/stan/math/prim/prob/skew_normal_lpdf.hpp +++ b/stan/math/prim/prob/skew_normal_lpdf.hpp @@ -7,16 +7,16 @@ #include #include #include -#include -#include -#include +#include #include #include #include #include +#include #include #include #include +#include #include namespace stan { @@ -27,7 +27,7 @@ template * = nullptr> inline return_type_t skew_normal_lpdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha) { + T_y&& y, T_loc&& mu, T_scale&& sigma, T_shape&& alpha) { 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; @@ -37,10 +37,10 @@ inline return_type_t skew_normal_lpdf( check_consistent_sizes(function, "Random variable", y, "Location parameter", mu, "Scale parameter", sigma, "Shape parameter", alpha); - T_y_ref y_ref = y; - T_mu_ref mu_ref = mu; - T_sigma_ref sigma_ref = sigma; - T_alpha_ref alpha_ref = alpha; + T_y_ref y_ref = std::forward(y); + T_mu_ref mu_ref = std::forward(mu); + T_sigma_ref sigma_ref = std::forward(sigma); + T_alpha_ref alpha_ref = std::forward(alpha); 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)); @@ -52,7 +52,7 @@ inline return_type_t skew_normal_lpdf( check_finite(function, "Shape parameter", alpha_val); check_positive(function, "Scale parameter", sigma_val); - if (size_zero(y, mu, sigma, alpha)) { + if (size_zero(y_ref, mu_ref, sigma_ref, alpha_ref)) { return 0.0; } if constexpr (!include_summand::value) { @@ -64,62 +64,48 @@ inline return_type_t skew_normal_lpdf( const auto& inv_sigma = to_ref_if>(inv(sigma_val)); - const auto& y_minus_mu_over_sigma - = to_ref_if::value>( - (y_val - mu_val) * inv_sigma); - const auto& log_erfc_alpha_z - = to_ref_if>( - log(erfc(-alpha_val * y_minus_mu_over_sigma * INV_SQRT_TWO))); + const auto& z = to_ref((y_val - mu_val) * inv_sigma); + const auto& az = to_ref(alpha_val * z); + const auto [values, slopes] = internal::std_normal_lcdf_value_grad< + is_any_autodiff_v>(az); - size_t N = max_size(y, mu, sigma, alpha); - T_partials_return logp = sum(log_erfc_alpha_z); + size_t N = max_size(y_ref, mu_ref, sigma_ref, alpha_ref); + T_partials_return logp = N * LOG_TWO + sum(values); if constexpr (include_summand::value) { logp -= HALF_LOG_TWO_PI * N; } if constexpr (include_summand::value) { - logp -= sum(log(sigma_val)) * N / math::size(sigma); + logp -= sum(log(sigma_val)) * N / math::size(sigma_ref); } if constexpr (include_summand::value) { - logp -= sum(square(y_minus_mu_over_sigma)) * 0.5 * N - / max_size(y, mu, sigma); + logp -= sum(square(z)) * 0.5 * N / max_size(y_ref, mu_ref, sigma_ref); } - - if constexpr (is_any_autodiff_v) { - const auto& sq = square(alpha_val * y_minus_mu_over_sigma * INV_SQRT_TWO); - const auto& ex = exp(-sq - log_erfc_alpha_z); - auto deriv_logerf = to_ref_if< - is_any_autodiff_v< - T_y, T_loc> + is_autodiff_v + is_autodiff_v >= 2>( - SQRT_TWO_OVER_SQRT_PI * ex); - if constexpr (is_any_autodiff_v) { - auto deriv_y_loc - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - (y_minus_mu_over_sigma - deriv_logerf * alpha_val) * inv_sigma); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = -deriv_y_loc; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials) = std::move(deriv_y_loc); - } + if constexpr (is_any_autodiff_v) { + const auto& score = to_ref_if< + (is_autodiff_v + is_autodiff_v + is_autodiff_v) + >= 2>((slopes * alpha_val - z) * inv_sigma); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = score; } - if constexpr (is_autodiff_v) { - edge<2>(ops_partials).partials_ - = ((y_minus_mu_over_sigma - deriv_logerf * alpha_val) - * y_minus_mu_over_sigma - - 1) - * inv_sigma; + if constexpr (is_autodiff_v) { + partials<1>(ops_partials) = -score; } - if constexpr (is_autodiff_v) { - partials<3>(ops_partials) = deriv_logerf * y_minus_mu_over_sigma; + if constexpr (is_autodiff_v) { + partials<2>(ops_partials) = -score * z - inv_sigma; } } + if constexpr (is_autodiff_v) { + partials<3>(ops_partials) = slopes * z; + } return ops_partials.build(logp); } template inline return_type_t skew_normal_lpdf( - const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha) { - return skew_normal_lpdf(y, mu, sigma, alpha); + T_y&& y, T_loc&& mu, T_scale&& sigma, T_shape&& alpha) { + return skew_normal_lpdf(std::forward(y), std::forward(mu), + std::forward(sigma), + std::forward(alpha)); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_ccdf_log.hpp b/stan/math/prim/prob/std_normal_ccdf_log.hpp index 70db6c52fc3..3d8e1b05f8e 100644 --- a/stan/math/prim/prob/std_normal_ccdf_log.hpp +++ b/stan/math/prim/prob/std_normal_ccdf_log.hpp @@ -12,7 +12,7 @@ namespace math { */ template inline return_type_t std_normal_ccdf_log(const T_y& y) { - return std_normal_lccdf(y); + return std_normal_lccdf(y); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_cdf.hpp b/stan/math/prim/prob/std_normal_cdf.hpp index 2fc86dc1d9f..487a461bb4a 100644 --- a/stan/math/prim/prob/std_normal_cdf.hpp +++ b/stan/math/prim/prob/std_normal_cdf.hpp @@ -2,17 +2,8 @@ #define STAN_MATH_PRIM_PROB_STD_NORMAL_CDF_HPP #include -#include -#include -#include -#include #include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { @@ -23,6 +14,9 @@ namespace math { * * \f$\Phi(x) = \frac{1}{\sqrt{2 \pi}} \int_{-\inf}^x e^{-t^2/2} dt\f$. * + * Evaluated as the exponential of the log cdf so the tails and gradients + * share the log kernel instead of hard cutoffs. + * * @tparam T_y type of y * @param y scalar variate * @return The standard normal cdf evaluated at the specified argument. @@ -30,58 +24,9 @@ namespace math { template < typename T_y, require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> -inline return_type_t std_normal_cdf(const T_y& y) { - using T_partials_return = partials_return_t; - using std::exp; - using T_y_ref = ref_type_t; - static constexpr const char* function = "std_normal_cdf"; - T_y_ref y_ref = y; - check_not_nan(function, "Random variable", y_ref); - - if (size_zero(y)) { - return 1.0; - } - - T_partials_return cdf(1.0); - auto ops_partials = make_partials_propagator(y_ref); - - scalar_seq_view y_vec(y_ref); - size_t N = stan::math::size(y); - - for (size_t n = 0; n < N; n++) { - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return scaled_y = y_dbl * INV_SQRT_TWO; - T_partials_return cdf_n; - if (y_dbl < -37.5) { - cdf_n = 0.0; - } else if (y_dbl < -5.0) { - cdf_n = 0.5 * erfc(-scaled_y); - } else if (y_dbl > 8.25) { - cdf_n = 1; - } else { - cdf_n = 0.5 * (1.0 + erf(scaled_y)); - } - - cdf *= cdf_n; - - if constexpr (is_autodiff_v) { - const T_partials_return rep_deriv - = (y_dbl < -37.5) - ? 0.0 - : INV_SQRT_TWO_PI * exp(-scaled_y * scaled_y) / cdf_n; - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] += rep_deriv; - } - } - } - - if constexpr (is_autodiff_v) { - for (size_t n = 0; n < N; ++n) { - partials<0>(ops_partials)[n] *= cdf; - } - } - - return ops_partials.build(cdf); +inline return_type_t std_normal_cdf(T_y&& y) { + return exp(internal::std_normal_lcdf_impl("std_normal_cdf", + std::forward(y))); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_lccdf.hpp b/stan/math/prim/prob/std_normal_lccdf.hpp index baf52f18156..2c51ea2b787 100644 --- a/stan/math/prim/prob/std_normal_lccdf.hpp +++ b/stan/math/prim/prob/std_normal_lccdf.hpp @@ -5,16 +5,13 @@ namespace stan { namespace math { -namespace internal { -constexpr char std_normal_lccdf_func[] = "std_normal_lccdf"; -} // namespace internal template < typename T_y, require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> -inline return_type_t std_normal_lccdf(const T_y& y) { - return std_normal_lcdf( - -as_array_or_scalar(y)); +inline return_type_t std_normal_lccdf(T_y&& y) { + return internal::std_normal_lcdf_impl("std_normal_lccdf", + std::forward(y)); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_lcdf.hpp b/stan/math/prim/prob/std_normal_lcdf.hpp index 738fe948765..ea517fe659e 100644 --- a/stan/math/prim/prob/std_normal_lcdf.hpp +++ b/stan/math/prim/prob/std_normal_lcdf.hpp @@ -4,237 +4,68 @@ #include #include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include #include #include -#include #include #include -#include -#include +#include namespace stan { namespace math { namespace internal { -constexpr char std_normal_lcdf_func[] = "std_normal_lcdf"; + +/** Log of the standard normal cdf, or of its complement when `reflect` is + * set: the input is negated so no negated autodiff operands are built. + */ +template +inline return_type_t std_normal_lcdf_impl(const char* function, T_y&& y) { + using T_y_ref = ref_type_if_not_constant_t; + constexpr double sign = reflect ? -1.0 : 1.0; + T_y_ref y_ref = std::forward(y); + decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref)); + check_not_nan(function, "Random variable", y_val); + + if (size_zero(y_ref)) { + return 0; + } + + // Branching avoids materialising sign * y_val when not reflecting. + const auto [values, slopes] = [&]() { + if constexpr (reflect) { + return std_normal_lcdf_value_grad>(-y_val); + } else { + return std_normal_lcdf_value_grad>(y_val); + } + }(); + auto ops_partials = make_partials_propagator(y_ref); + if constexpr (is_autodiff_v) { + partials<0>(ops_partials) = sign * slopes; + } + return ops_partials.build(sum(values)); +} + } // namespace internal /** \ingroup prob_dists * @brief Calculates the log of the cdf of the standard normal distribution * - * The piecewise structure here is the same as `normal_lcdf`, and the - * cutoffs, their provenance and the measurements behind them are documented - * once, on that function. See prim/prob/normal_lcdf.hpp. That covers the - * A&S 7.1.26, Cody (1969) and DLMF 7.12.1 references, the R `pnorm` and - * SciPy `log_ndtr` cross-references, why the erfc/Cody crossover sits at 4, - * the stan-dev/math#1411 origin of the interior Taylor cutoffs, and the two - * cutoff tables. + * Shares the scalar value and slope calculation with normal_lcdf through + * std_normal_lcdf_impl.hpp. * - * Two differences apply when reading it here. The scaled variable is - * `scaled_y = y * INV_SQRT_TWO`, not `scaled_diff = (y - mu) / (sigma * - * SQRT_TWO)`; since `mu = 0` and `sigma = 1` the two coincide, so every - * cutoff value transfers unchanged. And the test that enforces the - * `worst in-range` column for this function is the branch_accuracy test in - * mix/prob/std_normal_cdf_log_test.cpp. - * - * @tparam func name reported by the error checks. Reflected distributions - * such as `std_normal_lccdf` delegate here and pass their own name so that - * exceptions name the function the user actually called. * @tparam T_y A vector or scalar type for the random variable. * @param y (Sequence of) scalar(s). * @return The log of the standard normal cdf evaluated at the specified * argument. If given a container, the log of the product of the cdfs. */ template < - const char* func = internal::std_normal_lcdf_func, typename T_y, + typename T_y, require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> -inline return_type_t std_normal_lcdf(const T_y& y) { - using T_partials_return = partials_return_t; - using std::exp; - using std::fabs; - using std::log; - using std::pow; - using T_y_ref = ref_type_t; - static constexpr const char* function = func; - T_y_ref y_ref = y; - check_not_nan(function, "Random variable", y_ref); - - if (size_zero(y)) { - return 0; - } - - T_partials_return lcdf(0.0); - auto ops_partials = make_partials_propagator(y_ref); - - scalar_seq_view y_vec(y_ref); - size_t N = stan::math::size(y); - - for (size_t n = 0; n < N; n++) { - const T_partials_return y_dbl = y_vec.val(n); - const T_partials_return scaled_y = y_dbl * INV_SQRT_TWO; - const T_partials_return x2 = square(scaled_y); - - // Rigorous numerical approximations are applied here to deal with values - // of |scaled_y|>>0. This is needed to deal with rare base-rate - // logistic regression problems where it is useful to use an alternative - // link function instead. - // - // use erfc() instead of erf() in order to retain precision - // since for x>0 erfc()->0 - if (scaled_y > 0.0) { - // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) - lcdf += log1p(-0.5 * erfc(scaled_y)); - if (!is_not_nan(lcdf)) { - lcdf = 0; - } - } else if (scaled_y > -4.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x); -4 is R pnorm's M_SQRT_32 - // crossover expressed in scaled_y, sqrt(32)/sqrt(2) - lcdf += log(erfc(-scaled_y)) + LOG_HALF; - } else if (10.0 * log(fabs(scaled_y)) - < log(std::numeric_limits::max())) { - // entering territory where erfc(-x)~0 - // need to use direct numerical approximation of lcdf instead - // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) - // CDF(x) = 1/2erfc(-x) - const T_partials_return x4 = pow(scaled_y, 4); - const T_partials_return x6 = pow(scaled_y, 6); - const T_partials_return x8 = pow(scaled_y, 8); - const T_partials_return x10 = pow(scaled_y, 10); - const T_partials_return temp_p - = 0.000658749161529837803157 + 0.0160837851487422766278 / x2 - + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6 - + 0.305326634961232344035 / x8 + 0.0163153871373020978498 / x10; - const T_partials_return temp_q - = -0.00233520497626869185443 - 0.0605183413124413191178 / x2 - - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6 - - 2.56852019228982242072 / x8 - 1.0 / x10; - lcdf += LOG_HALF + log(INV_SQRT_PI + (temp_p / temp_q) / x2) - - log(-scaled_y) - x2; - } else { - // scaled_y^10 term will overflow - lcdf = stan::math::negative_infinity(); - } - - if constexpr (is_autodiff_v) { - // compute partial derivatives - // based on analytic form given by: - // dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2+erf(x)/2) - T_partials_return dnlcdf = 0.0; - T_partials_return t = 0.0; - T_partials_return t2 = 0.0; - T_partials_return t4 = 0.0; - - // calculate using piecewise function - // (due to instability / inaccuracy in the various approximations) - if (scaled_y > 2.9) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - t = 1.0 / (1.0 + 0.3275911 * scaled_y); - t2 = square(t); - t4 = pow(t, 4); - // A&S 7.1.26 puts exp(-x2) in the numerator; keep it there so it - // underflows to zero instead of overflowing inside a denominator - const T_partials_return exp_m_x2 = exp(-x2); - dnlcdf = INV_SQRT_PI * exp_m_x2 - / (1.0 - - exp_m_x2 - * (0.254829592 - 0.284496736 * t + 1.421413741 * t2 - - 1.453152027 * t2 * t + 1.061405429 * t4)); - } else if (scaled_y > 2.5) { - // in the trouble area where all of the standard numerical - // approximations are unstable - bridge the gap using Taylor - // expansions of the analytic function - // use Taylor expansion centred around x=2.7 - t = scaled_y - 2.7; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2 - - 0.008029540137 * t2 * t + 0.008232190507 * t4 - - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6); - } else if (scaled_y > 2.1) { - // use Taylor expansion centred around x=2.3 - t = scaled_y - 2.3; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2 - - 0.03326906904 * t2 * t + 0.02482478940 * t4 - - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6); - } else if (scaled_y > 1.5) { - // use Taylor expansion centred around x=1.85 - t = scaled_y - 1.85; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2 - - 0.09274533184 * t2 * t + 0.03543327418 * t4 - + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6); - } else if (scaled_y > 0.8) { - // use Taylor expansion centred around x=1.15 - t = scaled_y - 1.15; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2 - - 0.09748053605 * t2 * t - 0.04347986191 * t4 - + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6); - } else if (scaled_y > 0.1) { - // use Taylor expansion centred around x=0.45 - t = scaled_y - 0.45; - t2 = square(t); - t4 = pow(t, 4); - dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2 - + 0.04700850676 * t2 * t - 0.03478651979 * t4 - - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6); - } else if (scaled_y < -29.0) { - // asymptotic Mills ratio, DLMF 7.12.1: dnlcdf grows linearly as - // -2*scaled_y, so no quadratic residual fit can track it - const T_partials_return inv_x2 = 1.0 / x2; - dnlcdf = -2.0 * scaled_y - / (1.0 + inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875))); - } else if (10.0 * log(fabs(scaled_y)) - < log(std::numeric_limits::max())) { - // approximation derived from Abramowitz and Stegun (1964) 7.1.26 - // use fact that erf(x)=-erf(-x) - // Abramowitz and Stegun define this for -inf) { - partials<0>(ops_partials)[n] += dnlcdf * INV_SQRT_TWO; - } - } - } - - return ops_partials.build(lcdf); +inline return_type_t std_normal_lcdf(T_y&& y) { + return internal::std_normal_lcdf_impl("std_normal_lcdf", + std::forward(y)); } } // namespace math diff --git a/stan/math/rev/fun/Phi.hpp b/stan/math/rev/fun/Phi.hpp index 06f62cb885a..da9d473eb6a 100644 --- a/stan/math/rev/fun/Phi.hpp +++ b/stan/math/rev/fun/Phi.hpp @@ -23,9 +23,7 @@ namespace math { \f[ \mbox{Phi}(x) = \begin{cases} - 0 & \mbox{if } x < -37.5 \\ - \Phi(x) & \mbox{if } -37.5 \leq x \leq 8.25 \\ - 1 & \mbox{if } x > 8.25 \\[6pt] + \Phi(x) & \mbox{if } -\infty \leq x \leq \infty \\[6pt] \textrm{error} & \mbox{if } x = \textrm{NaN} \end{cases} \f] diff --git a/stan/math/rev/fun/as_array_or_scalar.hpp b/stan/math/rev/fun/as_array_or_scalar.hpp index fddebcad31d..7fea607dca3 100644 --- a/stan/math/rev/fun/as_array_or_scalar.hpp +++ b/stan/math/rev/fun/as_array_or_scalar.hpp @@ -18,7 +18,11 @@ namespace math { */ template * = nullptr> inline auto as_array_or_scalar(T&& v) { - return v.array(); + if constexpr (is_eigen_array>::value) { + return std::forward(v); + } else { + return v.array(); + } } } // namespace math diff --git a/test/unit/math/mix/prob/exp_mod_normal_ccdf_log_test.cpp b/test/unit/math/mix/prob/exp_mod_normal_ccdf_log_test.cpp new file mode 100644 index 00000000000..197f2bf4c13 --- /dev/null +++ b/test/unit/math/mix/prob/exp_mod_normal_ccdf_log_test.cpp @@ -0,0 +1,30 @@ +#include +#include + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lccdf) { + auto f = [](const auto& p) { + return stan::math::exp_mod_normal_lccdf(p[0], p[1], p[2], p[3]); + }; + + for (double y : {-40.0, -10.0, -3.0, 0.3, 4.0, 60.0}) { + stan::test::expect_ad(f, Eigen::Vector4d(y, 0.1, 1.3, 0.7)); + } + stan::test::expect_ad( + [](const auto& mu, const auto& sigma, const auto& lambda) { + return stan::math::exp_mod_normal_lccdf(Eigen::Vector3d(-40, 0.3, 60), + mu, sigma, lambda); + }, + 0.1, 1.3, 0.7); +} + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lccdf_tails_and_endpoints) { + using stan::math::exp_mod_normal_lccdf; + using stan::math::INFTY; + using stan::math::NEGATIVE_INFTY; + // MPFR references for mu = 0, sigma = 1, lambda = 1. + EXPECT_NEAR(-1.353310396e-10, exp_mod_normal_lccdf(-6.0, 0.0, 1.0, 1.0), + 1e-18); + EXPECT_FLOAT_EQ(-799.5, exp_mod_normal_lccdf(800.0, 0.0, 1.0, 1.0)); + EXPECT_EQ(0, exp_mod_normal_lccdf(NEGATIVE_INFTY, 0.0, 1.0, 1.0)); + EXPECT_EQ(NEGATIVE_INFTY, exp_mod_normal_lccdf(INFTY, 0.0, 1.0, 1.0)); +} diff --git a/test/unit/math/mix/prob/exp_mod_normal_cdf_log_test.cpp b/test/unit/math/mix/prob/exp_mod_normal_cdf_log_test.cpp new file mode 100644 index 00000000000..cd5f17cc747 --- /dev/null +++ b/test/unit/math/mix/prob/exp_mod_normal_cdf_log_test.cpp @@ -0,0 +1,38 @@ +#include +#include + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lcdf) { + auto f = [](const auto& p) { + return stan::math::exp_mod_normal_lcdf(p[0], p[1], p[2], p[3]); + }; + + for (double y : {-10.0, -3.0, 0.3, 4.0, 60.0}) { + stan::test::expect_ad(f, Eigen::Vector4d(y, 0.1, 1.3, 0.7)); + } + // Deep lower tail: log_diff_exp of two ~-480 terms is ill-conditioned for + // the finite-difference third derivatives. + stan::test::ad_tolerances tols; + tols.grad_hessian_grad_hessian_ = 5e-2; + stan::test::expect_ad(tols, f, Eigen::Vector4d(-40, 0.1, 1.3, 0.7)); + stan::test::expect_ad( + tols, + [](const auto& mu, const auto& sigma, const auto& lambda) { + return stan::math::exp_mod_normal_lcdf(Eigen::Vector3d(-40, 0.3, 60), + mu, sigma, lambda); + }, + 0.1, 1.3, 0.7); +} + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lcdf_tails_and_endpoints) { + using stan::math::exp_mod_normal_lcdf; + using stan::math::INFTY; + using stan::math::NEGATIVE_INFTY; + // MPFR references for mu = 0, sigma = 1, lambda = 1. + EXPECT_NEAR(-808.3232158, exp_mod_normal_lcdf(-40.0, 0.0, 1.0, 1.0), 1e-6); + EXPECT_NEAR(-55.64594635, exp_mod_normal_lcdf(-10.0, 0.0, 1.0, 1.0), 1e-7); + EXPECT_NEAR(-22.72329719, exp_mod_normal_lcdf(-6.0, 0.0, 1.0, 1.0), 1e-7); + EXPECT_NEAR(-1.443704555e-26, exp_mod_normal_lcdf(60.0, 0.0, 1.0, 1.0), + 1e-34); + EXPECT_EQ(NEGATIVE_INFTY, exp_mod_normal_lcdf(NEGATIVE_INFTY, 0.0, 1.0, 1.0)); + EXPECT_EQ(0, exp_mod_normal_lcdf(INFTY, 0.0, 1.0, 1.0)); +} diff --git a/test/unit/math/mix/prob/exp_mod_normal_cdf_test.cpp b/test/unit/math/mix/prob/exp_mod_normal_cdf_test.cpp new file mode 100644 index 00000000000..4b51b692f6c --- /dev/null +++ b/test/unit/math/mix/prob/exp_mod_normal_cdf_test.cpp @@ -0,0 +1,20 @@ +#include +#include + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_cdf) { + auto f = [](const auto& p) { + return stan::math::exp_mod_normal_cdf(p[0], p[1], p[2], p[3]); + }; + + for (double y : {-40.0, -10.0, -3.0, 0.3, 4.0, 60.0}) { + stan::test::expect_ad(f, Eigen::Vector4d(y, 0.1, 1.3, 0.7)); + } +} + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_cdf_endpoints) { + using stan::math::exp_mod_normal_cdf; + using stan::math::INFTY; + using stan::math::NEGATIVE_INFTY; + EXPECT_EQ(0, exp_mod_normal_cdf(NEGATIVE_INFTY, 0.0, 1.0, 1.0)); + EXPECT_EQ(1, exp_mod_normal_cdf(INFTY, 0.0, 1.0, 1.0)); +} diff --git a/test/unit/math/mix/prob/exp_mod_normal_test.cpp b/test/unit/math/mix/prob/exp_mod_normal_test.cpp new file mode 100644 index 00000000000..dd3dabade74 --- /dev/null +++ b/test/unit/math/mix/prob/exp_mod_normal_test.cpp @@ -0,0 +1,25 @@ +#include +#include + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lpdf) { + auto f = [](const auto& p) { + return stan::math::exp_mod_normal_lpdf(p[0], p[1], p[2], p[3]); + }; + + for (double y : {-40.0, -10.0, -3.0, 0.3, 4.0, 60.0}) { + stan::test::expect_ad(f, Eigen::Vector4d(y, 0.1, 1.3, 0.7)); + } + stan::test::expect_ad(f, Eigen::Vector4d(0, 0, 2, 20)); + stan::test::expect_ad( + [](const auto& mu, const auto& sigma, const auto& lambda) { + return stan::math::exp_mod_normal_lpdf(Eigen::Vector3d(-40, 0.3, 60), + mu, sigma, lambda); + }, + 0.1, 1.3, 0.7); +} + +TEST_F(AgradRev, mathMixScalFun_exp_mod_normal_lpdf_large_lambda) { + // MPFR reference; erfc underflows here. + EXPECT_NEAR(-1.6127097401997972, + stan::math::exp_mod_normal_lpdf(0.0, 0.0, 2.0, 20.0), 1e-12); +} diff --git a/test/unit/math/mix/prob/lognormal_ccdf_log_test.cpp b/test/unit/math/mix/prob/lognormal_ccdf_log_test.cpp new file mode 100644 index 00000000000..fee332939c1 --- /dev/null +++ b/test/unit/math/mix/prob/lognormal_ccdf_log_test.cpp @@ -0,0 +1,24 @@ +#include +#include +#include + +TEST_F(AgradRev, mathMixScalFun_lognormal_lccdf) { + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::lognormal_lccdf(y, mu, sigma); + }; + + stan::test::expect_ad(f, 2.0, 0.5, 1.5); + const Eigen::Vector3d y(1, std::exp(-50.0), std::exp(50.0)); + stan::test::expect_ad( + [&](const auto& mu, const auto& sigma) { + return stan::math::lognormal_lccdf(y, mu, sigma); + }, + 0.0, 1.0); +} + +TEST_F(AgradRev, mathMixScalFun_lognormal_lccdf_tail_and_endpoints) { + using stan::math::lognormal_lccdf; + EXPECT_NEAR(-1254.8313611394199, lognormal_lccdf(std::exp(50.0), 0.0, 1.0), + 2e-12); + EXPECT_EQ(0, lognormal_lccdf(0.0, 0.0, 1.0)); +} diff --git a/test/unit/math/mix/prob/lognormal_cdf_log_test.cpp b/test/unit/math/mix/prob/lognormal_cdf_log_test.cpp new file mode 100644 index 00000000000..b19b38f6dd3 --- /dev/null +++ b/test/unit/math/mix/prob/lognormal_cdf_log_test.cpp @@ -0,0 +1,36 @@ +#include +#include +#include + +TEST_F(AgradRev, mathMixScalFun_lognormal_lcdf) { + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::lognormal_lcdf(y, mu, sigma); + }; + + stan::test::expect_ad(f, 2.0, 0.5, 1.5); + stan::test::expect_ad(f, 1.0, 50.0, 1.0); + const Eigen::Vector3d y(1, std::exp(-50.0), 2); + stan::test::expect_ad( + [&](const auto& mu, const auto& sigma) { + return stan::math::lognormal_lcdf(y, mu, sigma); + }, + 0.0, 1.0); +} + +TEST_F(AgradRev, mathMixScalFun_lognormal_lcdf_tail_and_endpoints) { + using stan::math::INFTY; + using stan::math::lognormal_lcdf; + using stan::math::NEGATIVE_INFTY; + using stan::math::var; + EXPECT_NEAR(-1254.8313611394199, lognormal_lcdf(std::exp(-50.0), 0.0, 1.0), + 2e-12); + EXPECT_EQ(NEGATIVE_INFTY, lognormal_lcdf(0.0, 0.0, 1.0)); + var sigma = 1; + auto lp = lognormal_lcdf(INFTY, 0.0, sigma); + lp.grad(); + EXPECT_EQ(0, lp.val()); + EXPECT_EQ(0, sigma.adj()); + EXPECT_THROW( + lognormal_lcdf(Eigen::VectorXd::Ones(2), Eigen::VectorXd::Ones(3), 1), + std::invalid_argument); +} diff --git a/test/unit/math/mix/prob/lognormal_cdf_test.cpp b/test/unit/math/mix/prob/lognormal_cdf_test.cpp new file mode 100644 index 00000000000..949129df6d9 --- /dev/null +++ b/test/unit/math/mix/prob/lognormal_cdf_test.cpp @@ -0,0 +1,25 @@ +#include +#include +#include + +TEST_F(AgradRev, mathMixScalFun_lognormal_cdf) { + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::lognormal_cdf(y, mu, sigma); + }; + + stan::test::expect_ad(f, 2.0, 0.5, 1.5); + const Eigen::Vector3d y(1, std::exp(-50.0), std::exp(50.0)); + stan::test::expect_ad( + [&](const auto& mu, const auto& sigma) { + return stan::math::lognormal_cdf(y, mu, sigma); + }, + 0.0, 1.0); +} + +TEST_F(AgradRev, mathMixScalFun_lognormal_cdf_tail_and_endpoints) { + using stan::math::lognormal_cdf; + EXPECT_NEAR(1, + lognormal_cdf(std::exp(-20.0), 0.0, 1.0) / 2.7536241186062337e-89, + 1e-12); + EXPECT_EQ(0, lognormal_cdf(0.0, 0.0, 1.0)); +} diff --git a/test/unit/math/mix/prob/normal_ccdf_log_test.cpp b/test/unit/math/mix/prob/normal_ccdf_log_test.cpp index f0fce7c3630..ab3d40a44e1 100644 --- a/test/unit/math/mix/prob/normal_ccdf_log_test.cpp +++ b/test/unit/math/mix/prob/normal_ccdf_log_test.cpp @@ -52,3 +52,26 @@ TEST_F(AgradRev, mathMixScalFun_normal_lccdf_branch_accuracy) { normal_lcdf_tail_test::expect_branch_accuracy(normal_lccdf_mix_test::fn, normal_lccdf_mix_test::dir); } + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_varmat) { + using stan::math::var; + using stan::math::var_value; + Eigen::VectorXd y(3); + y << -1.5, 0.2, 3.0; + Eigen::Matrix y_mat = y; + var mu_mat = 0.5, sigma_mat = 1.2; + var lp_mat = stan::math::normal_lccdf(y_mat, mu_mat, sigma_mat); + lp_mat.grad(); + const Eigen::VectorXd y_adj = y_mat.adj(); + const double mu_adj = mu_mat.adj(), sigma_adj = sigma_mat.adj(); + stan::math::set_zero_all_adjoints(); + + var_value y_var(y); + var mu = 0.5, sigma = 1.2; + var lp = stan::math::normal_lccdf(y_var, mu, sigma); + lp.grad(); + EXPECT_DOUBLE_EQ(lp_mat.val(), lp.val()); + EXPECT_MATRIX_NEAR(y_adj, y_var.adj(), 1e-14); + EXPECT_DOUBLE_EQ(mu_adj, mu.adj()); + EXPECT_DOUBLE_EQ(sigma_adj, sigma.adj()); +} diff --git a/test/unit/math/mix/prob/normal_cdf_log_test.cpp b/test/unit/math/mix/prob/normal_cdf_log_test.cpp index b17a93685d5..3b7cd0f88e6 100644 --- a/test/unit/math/mix/prob/normal_cdf_log_test.cpp +++ b/test/unit/math/mix/prob/normal_cdf_log_test.cpp @@ -52,3 +52,16 @@ TEST_F(AgradRev, mathMixScalFun_normal_lcdf_branch_accuracy) { normal_lcdf_tail_test::expect_branch_accuracy(normal_lcdf_mix_test::fn, normal_lcdf_mix_test::dir); } + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_infinite_endpoints) { + using stan::math::INFTY; + using stan::math::NEGATIVE_INFTY; + using stan::math::normal_lcdf; + using stan::math::var; + EXPECT_EQ(0, normal_lcdf(INFTY, 0.0, 1.0)); + EXPECT_EQ(NEGATIVE_INFTY, normal_lcdf(NEGATIVE_INFTY, 0.0, 1.0)); + var sigma = 1; + auto lp = normal_lcdf(INFTY, 0.0, sigma); + lp.grad(); + EXPECT_EQ(0, sigma.adj()); +} diff --git a/test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp b/test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp index 4cadb205dd5..09b3c99a7d6 100644 --- a/test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp +++ b/test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp @@ -5,6 +5,7 @@ #include #include #include +#include /** * Shared checks for the tails of normal_lcdf, std_normal_lcdf and their @@ -18,11 +19,50 @@ * * References were evaluated at 60 significant digits from * `Phi(y) = erfc(-y/sqrt(2))/2` and `phi(y) = exp(-y^2/2)/sqrt(2 pi)`, then - * rounded to double. The cutoff tables they enforce are documented on - * stan::math::normal_lcdf in prim/prob/normal_lcdf.hpp. + * rounded to double. The historical approximation cutoffs are retained as + * regression inputs; the shared kernel now uses erfc and a Cody tail instead + * of the former Taylor expansions and residual fits. */ namespace normal_lcdf_tail_test { +template +void check_tail_derivatives(const F& f, double z, + const std::array& expected) { + using namespace stan::math; + const auto check = [&](double actual, size_t order) { + EXPECT_NEAR(expected[order], actual, 1e-12 * std::abs(expected[order])); + }; + { + nested_rev_autodiff nested; + fvar> x; + x.val_.val_ = z; + x.val_.d_ = 1; + x.d_.val_ = 1; + auto y = f(x); + y.d_.d_.grad(); + check(y.d_.val_.val(), 0); + check(y.d_.d_.val(), 1); + check(x.val_.val_.adj(), 2); + } + { + nested_rev_autodiff nested; + fvar x(z, 1); + auto y = f(x); + y.d_.grad(); + check(y.d_.val(), 0); + check(x.val_.adj(), 1); + } + fvar>> x; + x.val_.val_.val_ = z; + x.val_.val_.d_ = 1; + x.val_.d_.val_ = 1; + x.d_.val_.val_ = 1; + const auto y = f(x); + check(y.d_.val_.val_, 0); + check(y.d_.d_.val_, 1); + check(y.d_.d_.d_, 2); +} + /** Sign relating the function under test to normal_lcdf. */ struct orientation { static constexpr double lcdf = 1.0; @@ -69,7 +109,7 @@ void value_and_grad(const F& f, double y, double* value, double* grad) { /** * Five-term DLMF 7.12.1 asymptotic for d/dy log Phi(y), valid in the far lower - * tail. Independent of the four-term form the implementation uses. + * tail. Independent of the Cody rational approximation used by the kernel. */ inline double mills_dlogphi_dy(double y) { const double s = y / std::sqrt(2.0); @@ -80,18 +120,14 @@ inline double mills_dlogphi_dy(double y) { } /** - * expect_ad on both sides of every interior cutoff of the piecewise value and - * derivative, given in scaled units so they line up with the table in - * prim/prob/normal_lcdf.hpp. + * expect_ad on both sides of the historical value and derivative cutoffs, + * including the current Cody crossover at scaled input -4. * * The 0.01 offset keeps clear of `finite_diff_grad_hessian_auto`'s stencil, * which spans about `1.2e-5 * max(1, |y|)`, so neither side straddles a seam. * - * This pins that every branch is executed and is accurate where used. It - * cannot pin the cutoff location: adjacent branches agree to well within 1e-4 - * for about 0.1 either side of a seam, by construction. Mutation-checked -- - * re-centring the 1.85 series on 1.80 fails, shifting the 2.5 boundary to 2.4 - * does not, 2.75 does. + * Keep the former cutoffs as regression inputs even though the Taylor and + * residual-fit branches have been replaced by the shared normal kernel. */ template void expect_ad_across_cutoffs(const F& f, double dir) { @@ -118,45 +154,42 @@ void expect_ad_at_defect_inputs(const F& f, double dir) { } /** - * Per-branch derivative accuracy, each row at that branch's own measured worst - * in-range error plus a small margin. Per-branch rather than blanket because - * expect_ad's 1e-4 `gradient_grad_` leaves the (0.8, 1.5] branch only 1.6x. + * Derivative references at inputs covered by the former Taylor branches. + * The shared kernel should agree with these references to 12 digits. */ template void expect_branch_accuracy(const F& f, double dir) { struct branch_ref { double y; double d1; - double tol; }; - const branch_ref cases[] - = {// (2.5, 2.9], Taylor centre 2.7, measured worst 3.27e-05 - {3.5496760415564683, 0.0007326476540693806, 5e-5}, - {3.818376618407357, 0.00027222779390121885, 5e-5}, - {4.1012193308819755, 8.881828792625139e-05, 5e-5}, - // (2.1, 2.5], Taylor centre 2.3, measured worst 2.80e-05 - {2.9839906166072305, 0.004655923761682003, 4e-5}, - {3.2526911934581184, 0.00201252166907532, 4e-5}, - {3.5355339059327378, 0.0007702965121768906, 4e-5}, - // (1.5, 2.1], Taylor centre 1.85, measured worst 2.30e-05 - {2.135462479183374, 0.04148009614764338, 3e-5}, - {2.5455844122715714, 0.015709826784999336, 3e-5}, - {2.9698484809835, 0.004856449376114489, 3e-5}, - // (0.8, 1.5], Taylor centre 1.15, measured worst 6.09e-05 - {1.145512985522207, 0.236841175835376, 9e-5}, - {1.6263455967290592, 0.1121292480767062, 9e-5}, - {2.121320343559643, 0.042773100995777136, 9e-5}, - // (0.1, 0.8], Taylor centre 0.45, measured worst 7.61e-06 - {0.15556349186104046, 0.7015595130271106, 1e-5}, - {0.6363961030678928, 0.4416330793820557, 1e-5}, - {1.1313708498984762, 0.24150063210766093, 1e-5}}; + const branch_ref cases[] = {// (2.5, 2.9], Taylor centre 2.7 + {3.5496760415564683, 0.0007326476540693806}, + {3.818376618407357, 0.00027222779390121885}, + {4.1012193308819755, 8.881828792625139e-05}, + // (2.1, 2.5], Taylor centre 2.3 + {2.9839906166072305, 0.004655923761682003}, + {3.2526911934581184, 0.00201252166907532}, + {3.5355339059327378, 0.0007702965121768906}, + // (1.5, 2.1], Taylor centre 1.85 + {2.135462479183374, 0.04148009614764338}, + {2.5455844122715714, 0.015709826784999336}, + {2.9698484809835, 0.004856449376114489}, + // (0.8, 1.5], Taylor centre 1.15 + {1.145512985522207, 0.236841175835376}, + {1.6263455967290592, 0.1121292480767062}, + {2.121320343559643, 0.042773100995777136}, + // (0.1, 0.8], Taylor centre 0.45 + {0.15556349186104046, 0.7015595130271106}, + {0.6363961030678928, 0.4416330793820557}, + {1.1313708498984762, 0.24150063210766093}}; for (const branch_ref& c : cases) { const double y = dir * c.y; double value = 0; double grad = 0; value_and_grad(f, y, &value, &grad); - EXPECT_LT(std::fabs(grad / (dir * c.d1) - 1.0), c.tol) + EXPECT_LT(std::fabs(grad / (dir * c.d1) - 1.0), 1e-12) << "derivative branch drifted at y = " << y; } } @@ -222,7 +255,7 @@ void expect_derivatives_finite(const F& f) { /** * Far-tail value and gradient against 60-digit references, then a sweep * against the five-term Mills asymptotic over the region the implementation - * handles with its own four-term form. + * evaluates with its Cody rational approximation. */ template void expect_far_tail_gradient(const F& f, double dir) { diff --git a/test/unit/math/mix/prob/skew_normal_test.cpp b/test/unit/math/mix/prob/skew_normal_test.cpp new file mode 100644 index 00000000000..205f9e82247 --- /dev/null +++ b/test/unit/math/mix/prob/skew_normal_test.cpp @@ -0,0 +1,22 @@ +#include +#include + +TEST_F(AgradRev, mathMixScalFun_skew_normal_lpdf) { + auto f = [](const auto& p) { + return stan::math::skew_normal_lpdf(p[0], p[1], p[2], p[3]); + }; + + stan::test::expect_ad(f, Eigen::Vector4d(0.3, 0.1, 1.2, 0.5)); + stan::test::expect_ad(f, Eigen::Vector4d(-50, 0, 1, 1)); + stan::test::expect_ad( + [](const auto& mu, const auto& sigma, const auto& alpha) { + return stan::math::skew_normal_lpdf(Eigen::Vector3d(-50, 0, 2), mu, + sigma, alpha); + }, + 0.1, 1.2, 0.5); +} + +TEST_F(AgradRev, mathMixScalFun_skew_normal_lpdf_tail) { + EXPECT_NEAR(-2505.0571524920647, + stan::math::skew_normal_lpdf(-50.0, 0.0, 1.0, 1.0), 4e-12); +} diff --git a/test/unit/math/mix/prob/std_normal_ccdf_log_test.cpp b/test/unit/math/mix/prob/std_normal_ccdf_log_test.cpp index e052be64810..04fdada6562 100644 --- a/test/unit/math/mix/prob/std_normal_ccdf_log_test.cpp +++ b/test/unit/math/mix/prob/std_normal_ccdf_log_test.cpp @@ -49,3 +49,21 @@ TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_branch_accuracy) { normal_lcdf_tail_test::expect_branch_accuracy(std_normal_lccdf_mix_test::fn, std_normal_lccdf_mix_test::dir); } + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_varmat) { + using stan::math::var; + using stan::math::var_value; + Eigen::VectorXd y(3); + y << -1.5, 0.2, 3.0; + Eigen::Matrix y_mat = y; + var lp_mat = stan::math::std_normal_lccdf(y_mat); + lp_mat.grad(); + const Eigen::VectorXd y_adj = y_mat.adj(); + stan::math::set_zero_all_adjoints(); + + var_value y_var(y); + var lp = stan::math::std_normal_lccdf(y_var); + lp.grad(); + EXPECT_DOUBLE_EQ(lp_mat.val(), lp.val()); + EXPECT_MATRIX_NEAR(y_adj, y_var.adj(), 1e-14); +} diff --git a/test/unit/math/mix/prob/std_normal_lcdf_impl_test.cpp b/test/unit/math/mix/prob/std_normal_lcdf_impl_test.cpp new file mode 100644 index 00000000000..10b035c6a00 --- /dev/null +++ b/test/unit/math/mix/prob/std_normal_lcdf_impl_test.cpp @@ -0,0 +1,55 @@ +#include +#include +#include +#include +#include + +using normal_lcdf_tail_test::check_tail_derivatives; + +TEST_F(AgradRev, std_normal_extreme_tail_derivatives) { + using namespace stan::math; + // High-precision references for L', L'', L''', where L(z)=log Phi(z). + // The last row uses L'(-a)~a, L''(-a)~-1, L'''(-a)~2/a^3. Larger |z| + // overflows the squared divisor in nested fvar division. + for (const auto& row : + {std::array{-1e4, 10000.000099999998, -0.9999999900000006, + 1.99999976000003e-12}, + {-1e8, 100000000.00000001, -0.9999999999999999, 1.9999999999999976e-24}, + {-1e10, 1e10, -1, 2e-30}}) { + SCOPED_TRACE(row[0]); + check_tail_derivatives([](const auto& z) { return std_normal_lcdf(z); }, + row[0], {row[1], row[2], row[3]}); + check_tail_derivatives( + [](const auto& z) { return normal_lcdf(z, 0.0, 1.0); }, row[0], + {row[1], row[2], row[3]}); + check_tail_derivatives([](const auto& z) { return std_normal_lccdf(z); }, + -row[0], {-row[1], row[2], -row[3]}); + check_tail_derivatives( + [](const auto& z) { return normal_lccdf(z, 0.0, 1.0); }, -row[0], + {-row[1], row[2], -row[3]}); + check_tail_derivatives( + [](const auto& z) { + return internal::std_normal_lcdf_value_grad(z).first; + }, + row[0], {row[1], row[2], row[3]}); + } +} + +TEST_F(AgradRev, std_normal_extreme_tail_vector_derivatives) { + using namespace stan::math; + const Eigen::Array z(-1e8, -1e10, -1e50, 0, 1e308); + const Eigen::Array expected(2e-24, 2e-30, 2e-150, + 0.21801361414499016, 0); + Eigen::Array>, 5, 1> x; + for (Eigen::Index i = 0; i < x.size(); ++i) { + x(i).val_.val_ = z(i); + x(i).val_.d_ = 1; + x(i).d_.val_ = 1; + } + auto y = std_normal_lcdf(x); + y.d_.d_.grad(); + for (Eigen::Index i = 0; i < x.size(); ++i) { + EXPECT_NEAR(expected(i), x(i).val_.val_.adj(), + 1e-12 * std::abs(expected(i))); + } +} diff --git a/test/unit/math/opencl/rev/lognormal_lcdf_test.cpp b/test/unit/math/opencl/rev/lognormal_lcdf_test.cpp index b3c9163b015..f5286f01b52 100644 --- a/test/unit/math/opencl/rev/lognormal_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/lognormal_lcdf_test.cpp @@ -158,4 +158,13 @@ TEST(ProbDistributionsLognormalLcdf, opencl_matches_cpu_big) { lognormal_lcdf_functor, y.transpose().eval(), mu.transpose().eval(), sigma.transpose().eval()); } + +TEST(ProbDistributionsLognormalLcdf, normal_tail) { + for (double y0 : {std::exp(-50.0), 1.0, std::exp(50.0), stan::math::INFTY}) { + SCOPED_TRACE(y0); + const Eigen::VectorXd y = Eigen::VectorXd::Constant(1, y0); + stan::math::test::compare_cpu_opencl_prim_rev(lognormal_lcdf_functor, y, + 0.0, 1.0); + } +} #endif diff --git a/test/unit/math/opencl/rev/normal_lcdf_test.cpp b/test/unit/math/opencl/rev/normal_lcdf_test.cpp index e264c5e1d9c..6bfab52af83 100644 --- a/test/unit/math/opencl/rev/normal_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/normal_lcdf_test.cpp @@ -79,6 +79,20 @@ TEST(ProbDistributionsNormalLcdf, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsNormalLcdf, opencl_matches_cpu_tail_branches) { + Eigen::VectorXd y(12); + y << -1e100, -40, -6, -4 * stan::math::SQRT_TWO, -1, 0, 0.3, 1, 4, 8, 40, + 1e100; + stan::math::test::compare_cpu_opencl_prim_rev(normal_lcdf_functor, y, 0.0, + 1.0); + // Check each value separately so the extreme tail cannot dominate the sum. + for (Eigen::Index i = 0; i < y.size(); ++i) { + SCOPED_TRACE(y[i]); + stan::math::test::compare_cpu_opencl_prim_rev( + normal_lcdf_functor, y.segment(i, 1).eval(), 0.0, 1.0); + } +} + TEST(ProbDistributionsNormalLcdf, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/opencl/rev/skew_normal_lpdf_test.cpp b/test/unit/math/opencl/rev/skew_normal_lpdf_test.cpp index 4e60315ee71..de457f68407 100644 --- a/test/unit/math/opencl/rev/skew_normal_lpdf_test.cpp +++ b/test/unit/math/opencl/rev/skew_normal_lpdf_test.cpp @@ -230,4 +230,12 @@ TEST(ProbDistributionsSkewNormal, opencl_matches_cpu_big) { sigma.transpose().eval(), alpha.transpose().eval()); } +TEST(ProbDistributionsSkewNormal, normal_tail) { + for (double y0 : {-50.0, -6.0, 0.0, 50.0}) { + SCOPED_TRACE(y0); + const Eigen::VectorXd y = Eigen::VectorXd::Constant(1, y0); + stan::math::test::compare_cpu_opencl_prim_rev(skew_normal_lpdf_functor, y, + 0.0, 1.0, 1.0); + } +} #endif diff --git a/test/unit/math/opencl/rev/std_normal_lcdf_test.cpp b/test/unit/math/opencl/rev/std_normal_lcdf_test.cpp index 52f3562c6fe..c07a6ee16a9 100644 --- a/test/unit/math/opencl/rev/std_normal_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/std_normal_lcdf_test.cpp @@ -46,4 +46,30 @@ TEST(ProbDistributionsStdNormalLcdf, opencl_matches_cpu_big) { stan::math::test::compare_cpu_opencl_prim_rev(std_normal_lcdf_functor, y.transpose().eval()); } + +TEST(ProbDistributionsStdNormalLcdf, opencl_matches_cpu_tail_branches) { + Eigen::VectorXd y(12); + y << -1e100, -40, -6, -4 * stan::math::SQRT_TWO, -1, 0, 0.3, 1, 4, 8, 40, + 1e100; + stan::math::test::compare_cpu_opencl_prim_rev(std_normal_lcdf_functor, y); + // Check each value separately so the extreme tail cannot dominate the sum. + for (Eigen::Index i = 0; i < y.size(); ++i) { + SCOPED_TRACE(y[i]); + stan::math::test::compare_cpu_opencl_prim_rev(std_normal_lcdf_functor, + y.segment(i, 1).eval()); + } +} + +TEST(ProbDistributionsStdNormalLcdf, empty_and_extreme_gradient) { + using namespace stan::math; + matrix_cl empty(Eigen::VectorXd(0)); + EXPECT_EQ(0, std_normal_lcdf(empty)); + EXPECT_EQ(0, std_normal_lccdf(empty)); + nested_rev_autodiff nested; + var_value> y( + to_matrix_cl(Eigen::VectorXd::Constant(1, -1.5e308))); + auto lp = std_normal_lcdf(y); + lp.grad(); + EXPECT_NEAR(1.5e308, from_matrix_cl(y.adj())(0, 0), 1.5e296); +} #endif diff --git a/test/unit/math/prim/fun/Phi_test.cpp b/test/unit/math/prim/fun/Phi_test.cpp index c43e260bf88..b97b498eba3 100644 --- a/test/unit/math/prim/fun/Phi_test.cpp +++ b/test/unit/math/prim/fun/Phi_test.cpp @@ -3,11 +3,10 @@ #include TEST(MathFunctions, Phi) { - EXPECT_EQ(0.5 + 0.5 * stan::math::erf(0.0), stan::math::Phi(0.0)); - EXPECT_FLOAT_EQ(0.5 + 0.5 * stan::math::erf(0.9 / std::sqrt(2.0)), - stan::math::Phi(0.9)); - EXPECT_EQ(0.5 + 0.5 * stan::math::erf(-5.0 / std::sqrt(2.0)), - stan::math::Phi(-5.0)); + EXPECT_NEAR(0.5, stan::math::Phi(0.0), 1e-16); + EXPECT_NEAR(0.81593987465324047, stan::math::Phi(0.9), 1e-15); + EXPECT_NEAR(1, stan::math::Phi(-5.0) / 2.8665157187919391e-07, 1e-14); + EXPECT_NEAR(1, stan::math::Phi(-3.0) / 0.0013498980316300946, 1e-14); } // tests calculating using R 3.0.2 Snow Leopard build (6558) diff --git a/test/unit/math/prim/prob/exp_mod_normal_ccdf_log_test.cpp b/test/unit/math/prim/prob/exp_mod_normal_ccdf_log_test.cpp index 5285b65e2d1..982cd9c1ef3 100644 --- a/test/unit/math/prim/prob/exp_mod_normal_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/exp_mod_normal_ccdf_log_test.cpp @@ -10,8 +10,7 @@ TEST(ProbExpModNormal, ccdf_log_matches_lccdf) { EXPECT_FLOAT_EQ((stan::math::exp_mod_normal_lccdf(y, mu, lambda, sigma)), (stan::math::exp_mod_normal_ccdf_log(y, mu, lambda, sigma))); EXPECT_FLOAT_EQ( - (stan::math::exp_mod_normal_lccdf( - y, mu, lambda, sigma)), + (stan::math::exp_mod_normal_lccdf(y, mu, lambda, sigma)), (stan::math::exp_mod_normal_ccdf_log( y, mu, lambda, sigma))); } diff --git a/test/unit/math/prim/prob/exp_mod_normal_cdf_log_test.cpp b/test/unit/math/prim/prob/exp_mod_normal_cdf_log_test.cpp index b57e3513258..35a191363a5 100644 --- a/test/unit/math/prim/prob/exp_mod_normal_cdf_log_test.cpp +++ b/test/unit/math/prim/prob/exp_mod_normal_cdf_log_test.cpp @@ -10,8 +10,7 @@ TEST(ProbExpModNormal, cdf_log_matches_lcdf) { EXPECT_FLOAT_EQ((stan::math::exp_mod_normal_lcdf(y, mu, lambda, sigma)), (stan::math::exp_mod_normal_cdf_log(y, mu, lambda, sigma))); EXPECT_FLOAT_EQ( - (stan::math::exp_mod_normal_lcdf( - y, mu, lambda, sigma)), + (stan::math::exp_mod_normal_lcdf(y, mu, lambda, sigma)), (stan::math::exp_mod_normal_cdf_log( y, mu, lambda, sigma))); } diff --git a/test/unit/math/prim/prob/lognormal_ccdf_log_test.cpp b/test/unit/math/prim/prob/lognormal_ccdf_log_test.cpp index 1d532ad0e46..d5ec3bb763e 100644 --- a/test/unit/math/prim/prob/lognormal_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/lognormal_ccdf_log_test.cpp @@ -9,6 +9,6 @@ TEST(ProbLognormal, ccdf_log_matches_lccdf) { EXPECT_FLOAT_EQ((stan::math::lognormal_lccdf(y, mu, sigma)), (stan::math::lognormal_ccdf_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::lognormal_lccdf(y, mu, sigma)), + (stan::math::lognormal_lccdf(y, mu, sigma)), (stan::math::lognormal_ccdf_log(y, mu, sigma))); } diff --git a/test/unit/math/prim/prob/lognormal_cdf_log_test.cpp b/test/unit/math/prim/prob/lognormal_cdf_log_test.cpp index aba7f32370f..eccf52e9ff2 100644 --- a/test/unit/math/prim/prob/lognormal_cdf_log_test.cpp +++ b/test/unit/math/prim/prob/lognormal_cdf_log_test.cpp @@ -9,6 +9,6 @@ TEST(ProbLognormal, cdf_log_matches_lcdf) { EXPECT_FLOAT_EQ((stan::math::lognormal_lcdf(y, mu, sigma)), (stan::math::lognormal_cdf_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::lognormal_lcdf(y, mu, sigma)), + (stan::math::lognormal_lcdf(y, mu, sigma)), (stan::math::lognormal_cdf_log(y, mu, sigma))); } diff --git a/test/unit/math/prim/prob/normal_ccdf_log_test.cpp b/test/unit/math/prim/prob/normal_ccdf_log_test.cpp index 5bab9655a32..09d53948938 100644 --- a/test/unit/math/prim/prob/normal_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/normal_ccdf_log_test.cpp @@ -9,7 +9,7 @@ TEST(ProbNormal, ccdf_log_matches_lccdf) { EXPECT_FLOAT_EQ((stan::math::normal_lccdf(y, mu, sigma)), (stan::math::normal_ccdf_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::normal_lccdf(y, mu, sigma)), + (stan::math::normal_lccdf(y, mu, sigma)), (stan::math::normal_ccdf_log(y, mu, sigma))); } diff --git a/test/unit/math/prim/prob/std_normal_ccdf_log_test.cpp b/test/unit/math/prim/prob/std_normal_ccdf_log_test.cpp index a2e39eb25b3..81a2d2a022c 100644 --- a/test/unit/math/prim/prob/std_normal_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/std_normal_ccdf_log_test.cpp @@ -6,7 +6,7 @@ TEST(ProbStdNormal, ccdf_log_matches_lccdf) { EXPECT_FLOAT_EQ((stan::math::std_normal_lccdf(y)), (stan::math::std_normal_ccdf_log(y))); - EXPECT_FLOAT_EQ((stan::math::std_normal_lccdf(y)), + EXPECT_FLOAT_EQ((stan::math::std_normal_lccdf(y)), (stan::math::std_normal_ccdf_log(y))); } diff --git a/test/unit/math/prim/prob/std_normal_lcdf_impl_test.cpp b/test/unit/math/prim/prob/std_normal_lcdf_impl_test.cpp new file mode 100644 index 00000000000..c2069517104 --- /dev/null +++ b/test/unit/math/prim/prob/std_normal_lcdf_impl_test.cpp @@ -0,0 +1,67 @@ +#include +#include +#include + +TEST(ProbStdNormal, scalar_tail_kernel) { + // High-precision references: input, log Phi(input), and its slope. + static constexpr std::array cases[] + = {{-1e150, -4.9999999999999995e299, 1e150}, + {-30, -454.32124395634321, 30.033259667433676}, + {-6, -20.736768949974707, 6.1584826045445986}, + {-1, -1.8410216450092636, 1.5251352761609811}, + {0, -0.69314718055994529, 0.79788456080286541}, + {1, -0.17275377902344988, 0.28759997093917838}, + {5.6568542494923797, -7.7086289798515121e-09, 4.4895039573954219e-08}, + {5.6568542494923806, -7.7086289798514724e-09, 4.4895039573953994e-08}, + {6, -9.865876455243758e-10, 6.0758828558176762e-09}, + {30, -4.9067139271481872e-198, 1.4736461348785476e-196}, + {1e8, 0, 0}, + {1e150, 0, 0}}; + for (const auto& row : cases) { + SCOPED_TRACE(row[0]); + const auto result + = stan::math::internal::std_normal_lcdf_value_grad(row[0]); + EXPECT_NEAR(row[1], result.first, 1e-12 * std::abs(row[1])); + EXPECT_NEAR(row[2], result.second, 1e-12 * std::abs(row[2])); + EXPECT_DOUBLE_EQ( + result.first, + stan::math::internal::std_normal_lcdf_value_grad(row[0]).first); + } +} + +TEST(ProbStdNormal, large_finite_log_density) { + EXPECT_TRUE(std::isfinite(stan::math::std_normal_lcdf(-1.5e154))); +} + +TEST(ProbStdNormal, vectorized_tail_kernel) { + using stan::math::internal::std_normal_lcdf_value_grad; + Eigen::ArrayXd z(25); + z << -1e308, -1e150, -50, -6, -4 * stan::math::SQRT_TWO, -5, -1, + -0.46875 * stan::math::SQRT_TWO, -0.01, 0, 0.01, + 0.46875 * stan::math::SQRT_TWO, 1, 3, 6, 20, 30, 37, 37.1, 38, 38.5, 40, + 50, 1e150, 1e308; + const auto result = std_normal_lcdf_value_grad(z); + const auto values = std_normal_lcdf_value_grad(z); + for (Eigen::Index i = 0; i < z.size(); ++i) { + SCOPED_TRACE(z[i]); + const auto scalar = std_normal_lcdf_value_grad(z[i]); + if (std::isfinite(scalar.first)) { + EXPECT_NEAR(scalar.first, result.first[i], + std::abs(scalar.first) * 1e-12 + 1e-323); + } else { + EXPECT_EQ(scalar.first, result.first[i]); + } + EXPECT_EQ(result.first[i], values.first[i]); + EXPECT_NEAR(scalar.second, result.second[i], + std::abs(scalar.second) * 1e-12 + 1e-323); + } +} + +TEST(ProbStdNormal, integer_vectors) { + using stan::math::std_normal_lcdf; + const Eigen::Vector3i z(-1, 0, 1); + EXPECT_NEAR(std_normal_lcdf(z.cast().eval()), std_normal_lcdf(z), + 1e-14); + EXPECT_NEAR(std_normal_lcdf(z), std_normal_lcdf(std::vector{-1, 0, 1}), + 1e-14); +}