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 49e308d7853..c48ce98acea 100644 --- a/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp +++ b/stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp @@ -29,7 +29,7 @@ static constexpr const char* std_normal_lcdf_device_function if (isnan(lcdf_n)) { lcdf_n = 0; } - } else if (scaled_y > -20.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)) { @@ -76,10 +76,15 @@ static constexpr const char* std_normal_lcdf_device_function t = 1.0 / (1.0 + 0.3275911 * scaled_y); t2 = t * t; t4 = pow(t, 4); - dnlcdf = 0.5 * M_2_SQRTPI - / (exp(x2) - 0.254829592 + 0.284496736 * t - - 1.421413741 * t2 + 1.453152027 * t2 * t - - 1.061405429 * t4); + // 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; @@ -116,6 +121,14 @@ static constexpr const char* std_normal_lcdf_device_function 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; @@ -124,10 +137,7 @@ static constexpr const char* std_normal_lcdf_device_function = M_2_SQRTPI / (0.254829592 * t - 0.284496736 * t2 + 1.421413741 * t2 * t - 1.453152027 * t4 + 1.061405429 * t4 * t); - if (scaled_y < -29.0) { - dnlcdf += 0.0015065154280332 * x2 - - 0.3993154819705530 * scaled_y - 4.2919418242931700; - } else if (scaled_y < -17.0) { + if (scaled_y < -17.0) { dnlcdf += 0.0001263257217272 * x2 * scaled_y + 0.0123586859488623 * x2 - 0.0860505264736028 * scaled_y - 1.252783383752970; diff --git a/stan/math/opencl/prim/normal_lccdf.hpp b/stan/math/opencl/prim/normal_lccdf.hpp index 775b9a9238b..25f4f6ed0d8 100644 --- a/stan/math/opencl/prim/normal_lccdf.hpp +++ b/stan/math/opencl/prim/normal_lccdf.hpp @@ -2,16 +2,13 @@ #define STAN_MATH_OPENCL_PRIM_NORMAL_LCCDF_HPP #ifdef STAN_OPENCL -#include -#include -#include -#include -#include -#include -#include +#include 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 @@ -33,77 +30,7 @@ template < require_any_not_stan_scalar_t* = nullptr> inline return_type_t normal_lccdf( const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) { - static constexpr const char* function = "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& 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); - matrix_cl one_m_erf = select( - scaled_diff < -37.5 * INV_SQRT_TWO, 2.0, - select(scaled_diff < -5.0 * INV_SQRT_TWO, 2.0 - erfc(-scaled_diff), - select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0, - 1.0 - erf(scaled_diff)))); - auto lccdf_expr = log(one_m_erf); - auto mu_deriv = select(scaled_diff > 8.25 * INV_SQRT_TWO, INFTY, - SQRT_TWO_OVER_SQRT_PI - * elt_divide(exp(-square(scaled_diff)), - elt_multiply(one_m_erf, sigma_val))); - auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff) * SQRT_TWO; - - matrix_cl lccdf_cl; - matrix_cl y_deriv_cl; - matrix_cl mu_deriv_cl; - matrix_cl sigma_deriv_cl; - - results(check_y_not_nan, check_mu_finite, check_sigma_positive) - = expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_expr); - results(lccdf_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl) - = expressions(lccdf_expr, calc_if>(y_deriv), - calc_if>(mu_deriv), - calc_if>(sigma_deriv)); - - T_partials_return lccdf - = LOG_HALF * lccdf_cl.size() + 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 normal_lcdf(-y, -mu, sigma); } } // namespace math diff --git a/stan/math/opencl/prim/normal_lcdf.hpp b/stan/math/opencl/prim/normal_lcdf.hpp index ccb755917ad..187eda05866 100644 --- a/stan/math/opencl/prim/normal_lcdf.hpp +++ b/stan/math/opencl/prim/normal_lcdf.hpp @@ -13,6 +13,7 @@ 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; @@ -29,7 +30,7 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY( if (isnan(normal_lcdf_n)) { normal_lcdf_n = 0; } - } else if (normal_lcdf_scaled_diff > -20.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)) { @@ -58,7 +59,8 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY( // 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 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) @@ -67,10 +69,15 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY( 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 - / (exp(x2) - 0.254829592 + 0.284496736 * t - 1.421413741 * t2 - + 1.453152027 * t2 * t - 1.061405429 * t4); + = 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 @@ -114,6 +121,17 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY( 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) @@ -128,11 +146,7 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY( - 1.453152027 * t4 + 1.061405429 * t4 * t); // check if we need to add a correction term // (from cubic fit of residuals) - if (normal_lcdf_deriv_scaled_diff < -29.0) { - normal_ldncdf += 0.0015065154280332 * x2 - - 0.3993154819705530 * normal_lcdf_deriv_scaled_diff - - 4.2919418242931700; - } else if (normal_lcdf_deriv_scaled_diff < -17.0) { + if (normal_lcdf_deriv_scaled_diff < -17.0) { normal_ldncdf += 0.0001263257217272 * x2 * normal_lcdf_deriv_scaled_diff + 0.0123586859488623 * x2 - 0.0860505264736028 * normal_lcdf_deriv_scaled_diff @@ -174,13 +188,14 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY( * @return The log of the product of densities. */ template < - typename T_y_cl, typename T_loc_cl, typename T_scale_cl, + const char* func = internal::normal_lcdf_opencl_func, 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) { - static constexpr const char* function = "normal_lcdf(OpenCL)"; + static constexpr const char* function = func; using std::isfinite; using std::isnan; diff --git a/stan/math/opencl/prim/std_normal_lccdf.hpp b/stan/math/opencl/prim/std_normal_lccdf.hpp index 7a3c90c7089..0804d5d9ad6 100644 --- a/stan/math/opencl/prim/std_normal_lccdf.hpp +++ b/stan/math/opencl/prim/std_normal_lccdf.hpp @@ -2,16 +2,13 @@ #define STAN_MATH_OPENCL_PRIM_STD_NORMAL_LCCDF_HPP #ifdef STAN_OPENCL -#include -#include -#include -#include -#include -#include -#include +#include 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 @@ -25,47 +22,7 @@ template * = nullptr, require_any_not_stan_scalar_t* = nullptr> inline return_type_t std_normal_lccdf(const T_y_cl& y) { - static constexpr const char* function = "std_normal_lccdf(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 one_m_erf - = select(y_val < -37.5, 2.0, - select(y_val < -5.0, 2.0 - erfc(-scaled_y), - select(y_val > 8.25, 0.0, 1.0 - erf(scaled_y)))); - auto lccdf_expr = colwise_sum(log(one_m_erf)); - auto y_deriv = -select( - y_val > 8.25, INFTY, - SQRT_TWO_OVER_SQRT_PI * elt_divide(exp(-square(scaled_y)), one_m_erf)); - - matrix_cl lccdf_cl; - matrix_cl y_deriv_cl; - - results(check_y_not_nan, lccdf_cl, y_deriv_cl) = expressions( - y_not_nan_expr, lccdf_expr, calc_if>(y_deriv)); - - T_partials_return lccdf = from_matrix_cl(lccdf_cl).sum() + LOG_HALF * N; - - auto ops_partials = make_partials_propagator(y_col); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = std::move(y_deriv_cl); - } - return ops_partials.build(lccdf); + return std_normal_lcdf(-y); } } // namespace math diff --git a/stan/math/opencl/prim/std_normal_lcdf.hpp b/stan/math/opencl/prim/std_normal_lcdf.hpp index 4a7db54fcd9..463b59b6f49 100644 --- a/stan/math/opencl/prim/std_normal_lcdf.hpp +++ b/stan/math/opencl/prim/std_normal_lcdf.hpp @@ -13,6 +13,9 @@ 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. @@ -21,11 +24,12 @@ namespace math { * @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 = "std_normal_lcdf(OpenCL)"; + static constexpr const char* function = func; using std::isfinite; using std::isnan; diff --git a/stan/math/prim/prob/normal_cdf_log.hpp b/stan/math/prim/prob/normal_cdf_log.hpp index 447a52846e2..8531c1bef76 100644 --- a/stan/math/prim/prob/normal_cdf_log.hpp +++ b/stan/math/prim/prob/normal_cdf_log.hpp @@ -14,7 +14,7 @@ template inline return_type_t normal_cdf_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { - return normal_lcdf(y, mu, sigma); + return normal_lcdf(y, mu, sigma); } } // namespace math diff --git a/stan/math/prim/prob/normal_lccdf.hpp b/stan/math/prim/prob/normal_lccdf.hpp index 123dced2a60..cbdd074418c 100644 --- a/stan/math/prim/prob/normal_lccdf.hpp +++ b/stan/math/prim/prob/normal_lccdf.hpp @@ -1,22 +1,13 @@ #ifndef STAN_MATH_PRIM_PROB_NORMAL_LCCDF_HPP #define STAN_MATH_PRIM_PROB_NORMAL_LCCDF_HPP -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include namespace stan { namespace math { +namespace internal { +constexpr char normal_lccdf_func[] = "normal_lccdf"; +} // namespace internal template normal_lccdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { - using T_partials_return = partials_return_t; - using std::exp; - using std::log; - 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_lccdf"; - 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 ccdf_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); - - T_partials_return one_m_erf; - if (scaled_diff < -37.5 * INV_SQRT_TWO) { - one_m_erf = 2.0; - } else if (scaled_diff < -5.0 * INV_SQRT_TWO) { - one_m_erf = 2.0 - erfc(-scaled_diff); - } else if (scaled_diff > 8.25 * INV_SQRT_TWO) { - one_m_erf = 0.0; - } else { - one_m_erf = 1.0 - erf(scaled_diff); - } - - ccdf_log += LOG_HALF + log(one_m_erf); - - if constexpr (is_any_autodiff_v) { - const T_partials_return rep_deriv_div_sigma - = scaled_diff > 8.25 * INV_SQRT_TWO - ? INFTY - : SQRT_TWO_OVER_SQRT_PI * exp(-scaled_diff * scaled_diff) - / one_m_erf / sigma_dbl; - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] -= rep_deriv_div_sigma; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += rep_deriv_div_sigma; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += rep_deriv_div_sigma * scaled_diff * SQRT_TWO; - } - } - } - return ops_partials.build(ccdf_log); + return normal_lcdf( + -as_array_or_scalar(y), -as_array_or_scalar(mu), sigma); } } // namespace math diff --git a/stan/math/prim/prob/normal_lcdf.hpp b/stan/math/prim/prob/normal_lcdf.hpp index f0a126608d1..b7446c9fd82 100644 --- a/stan/math/prim/prob/normal_lcdf.hpp +++ b/stan/math/prim/prob/normal_lcdf.hpp @@ -7,12 +7,15 @@ #include #include #include +#include #include #include #include +#include #include #include #include +#include #include #include #include @@ -20,23 +23,131 @@ namespace stan { namespace math { +namespace internal { +constexpr char normal_lcdf_func[] = "normal_lcdf"; +} // namespace internal -template 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. + * + * 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. + * @param y (Sequence of) scalar(s). + * @param mu (Sequence of) scalar(s). + * @param sigma (Sequence of) scalar(s). + * @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 std::exp; - using std::fabs; - using std::log; - using std::pow; - using std::sqrt; 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_lcdf"; + 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; @@ -66,7 +177,6 @@ inline return_type_t normal_lcdf(const T_y& y, const T_partials_return scaled_diff = (y_dbl - mu_dbl) / (sigma_dbl * SQRT_TWO); - const T_partials_return sigma_sqrt2 = sigma_dbl * SQRT_TWO; const T_partials_return x2 = square(scaled_diff); // Rigorous numerical approximations are applied here to deal with values @@ -82,8 +192,9 @@ inline return_type_t normal_lcdf(const T_y& y, if (!is_not_nan(cdf_log)) { cdf_log = 0; } - } else if (scaled_diff > -20.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x) + } 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())) { @@ -126,11 +237,16 @@ inline return_type_t normal_lcdf(const T_y& y, 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 - = 1.0 + = exp_m_x2 / (SQRT_PI - * (exp(x2) - 0.254829592 + 0.284496736 * t - 1.421413741 * t2 - + 1.453152027 * t2 * t - 1.061405429 * t4)); + * (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 @@ -174,6 +290,13 @@ inline return_type_t normal_lcdf(const T_y& y, 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 @@ -190,10 +313,7 @@ inline return_type_t normal_lcdf(const T_y& y, - 1.453152027 * t4 + 1.061405429 * t4 * t)); // check if we need to add a correction term // (from cubic fit of residuals) - if (scaled_diff < -29.0) { - dncdf_log += 0.0015065154280332 * x2 - - 0.3993154819705530 * scaled_diff - 4.2919418242931700; - } else if (scaled_diff < -17.0) { + if (scaled_diff < -17.0) { dncdf_log += 0.0001263257217272 * x2 * scaled_diff + 0.0123586859488623 * x2 - 0.0860505264736028 * scaled_diff - 1.252783383752970; @@ -213,7 +333,7 @@ inline return_type_t normal_lcdf(const T_y& y, } else { dncdf_log = stan::math::positive_infinity(); } - + const T_partials_return sigma_sqrt2 = sigma_dbl * SQRT_TWO; if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] += dncdf_log / sigma_sqrt2; } diff --git a/stan/math/prim/prob/std_normal_cdf_log.hpp b/stan/math/prim/prob/std_normal_cdf_log.hpp index 4c0bdaeea5c..14ec45f27dd 100644 --- a/stan/math/prim/prob/std_normal_cdf_log.hpp +++ b/stan/math/prim/prob/std_normal_cdf_log.hpp @@ -12,7 +12,7 @@ namespace math { */ template inline return_type_t std_normal_cdf_log(const T_y& y) { - return std_normal_lcdf(y); + return std_normal_lcdf(y); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_lccdf.hpp b/stan/math/prim/prob/std_normal_lccdf.hpp index 9ee17d204b0..baf52f18156 100644 --- a/stan/math/prim/prob/std_normal_lccdf.hpp +++ b/stan/math/prim/prob/std_normal_lccdf.hpp @@ -1,72 +1,20 @@ #ifndef STAN_MATH_PRIM_PROB_STD_NORMAL_LCCDF_HPP #define STAN_MATH_PRIM_PROB_STD_NORMAL_LCCDF_HPP -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include 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) { - using T_partials_return = partials_return_t; - using std::exp; - using std::log; - using T_y_ref = ref_type_t; - static constexpr const char* function = "std_normal_lccdf"; - T_y_ref y_ref = y; - check_not_nan(function, "Random variable", y_ref); - - if (size_zero(y)) { - return 0; - } - - T_partials_return lccdf(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; - - T_partials_return one_m_erf; - if (y_dbl < -37.5) { - one_m_erf = 2.0; - } else if (y_dbl < -5.0) { - one_m_erf = 2.0 - erfc(-scaled_y); - } else if (y_dbl > 8.25) { - one_m_erf = 0.0; - } else { - one_m_erf = 1.0 - erf(scaled_y); - } - - lccdf += LOG_HALF + log(one_m_erf); - - if constexpr (is_autodiff_v) { - const T_partials_return rep_deriv - = y_dbl > 8.25 - ? INFTY - : SQRT_TWO_OVER_SQRT_PI * exp(-scaled_y * scaled_y) / one_m_erf; - partials<0>(ops_partials)[n] -= rep_deriv; - } - } - - return ops_partials.build(lccdf); + return std_normal_lcdf( + -as_array_or_scalar(y)); } } // namespace math diff --git a/stan/math/prim/prob/std_normal_lcdf.hpp b/stan/math/prim/prob/std_normal_lcdf.hpp index 2946ecd7742..738fe948765 100644 --- a/stan/math/prim/prob/std_normal_lcdf.hpp +++ b/stan/math/prim/prob/std_normal_lcdf.hpp @@ -20,9 +20,38 @@ namespace stan { namespace math { +namespace internal { +constexpr char std_normal_lcdf_func[] = "std_normal_lcdf"; +} // 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. + * + * 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 < - typename T_y, + const char* func = internal::std_normal_lcdf_func, 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; @@ -31,7 +60,7 @@ inline return_type_t std_normal_lcdf(const T_y& y) { using std::log; using std::pow; using T_y_ref = ref_type_t; - static constexpr const char* function = "std_normal_lcdf"; + static constexpr const char* function = func; T_y_ref y_ref = y; check_not_nan(function, "Random variable", y_ref); @@ -63,8 +92,9 @@ inline return_type_t std_normal_lcdf(const T_y& y) { if (!is_not_nan(lcdf)) { lcdf = 0; } - } else if (scaled_y > -20.0) { - // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x) + } 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())) { @@ -107,9 +137,14 @@ inline return_type_t std_normal_lcdf(const T_y& y) { t = 1.0 / (1.0 + 0.3275911 * scaled_y); t2 = square(t); t4 = pow(t, 4); - dnlcdf = INV_SQRT_PI - / (exp(x2) - 0.254829592 + 0.284496736 * t - 1.421413741 * t2 - + 1.453152027 * t2 * t - 1.061405429 * t4); + // 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 @@ -153,6 +188,12 @@ inline return_type_t std_normal_lcdf(const T_y& y) { 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 @@ -167,10 +208,7 @@ inline return_type_t std_normal_lcdf(const T_y& y) { - 1.453152027 * t4 + 1.061405429 * t4 * t); // check if we need to add a correction term // (from cubic fit of residuals) - if (scaled_y < -29.0) { - dnlcdf += 0.0015065154280332 * x2 - 0.3993154819705530 * scaled_y - - 4.2919418242931700; - } else if (scaled_y < -17.0) { + 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) { 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 e573e7da289..f0fce7c3630 100644 --- a/test/unit/math/mix/prob/normal_ccdf_log_test.cpp +++ b/test/unit/math/mix/prob/normal_ccdf_log_test.cpp @@ -1,25 +1,54 @@ #include #include +#include TEST_F(AgradRev, mathMixScalFun_normal_lccdf) { - auto f = [](const double mu, const double sigma) { - return - [=](const auto& y) { return stan::math::normal_lccdf(y, mu, sigma); }; + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::normal_lccdf(y, mu, sigma); }; - stan::test::expect_ad(f(0.0, 1.0), -50.0); - stan::test::expect_ad(f(0.0, 1.0), -20.0 * stan::math::SQRT_TWO); - stan::test::expect_ad(f(0.0, 1.0), -5.5); - stan::test::expect_ad(f(0.0, 1.0), 0.0); - stan::test::expect_ad(f(0.0, 1.0), 0.15); - stan::test::expect_ad(f(0.0, 1.0), 1.14); - stan::test::expect_ad(f(0.0, 1.0), 3.00); - stan::test::expect_ad(f(0.0, 1.0), 10.00); - stan::test::expect_ad(f(-1.0, 2.0), 1.50); - stan::test::expect_ad(f(2.0, 1.0), 0.50); - - // third order autodiff tests can fail at borders of piecewise function - stan::test::ad_tolerances tols; - tols.grad_hessian_grad_hessian_ = 1e1; - stan::test::expect_ad(tols, f(0.0, 1.0), 0.1 * stan::math::SQRT_TWO); + stan::test::expect_ad(f, 50.0, 0.0, 1.0); + stan::test::expect_ad(f, 20.0 * stan::math::SQRT_TWO, 0.0, 1.0); + stan::test::expect_ad(f, 5.5, 0.0, 1.0); + stan::test::expect_ad(f, 0.0, 0.0, 1.0); + stan::test::expect_ad(f, -0.15, 0.0, 1.0); + stan::test::expect_ad(f, -1.14, 0.0, 1.0); + stan::test::expect_ad(f, -3.00, 0.0, 1.0); + stan::test::expect_ad(f, -10.00, 0.0, 1.0); + stan::test::expect_ad(f, -3.50, -1.0, 2.0); + stan::test::expect_ad(f, 3.50, 2.0, 1.0); +} + +namespace normal_lccdf_mix_test { +auto fn = [](const auto& y) { return stan::math::normal_lccdf(y, 0.0, 1.0); }; +constexpr double dir = normal_lcdf_tail_test::orientation::lccdf; +} // namespace normal_lccdf_mix_test + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_defect_inputs) { + normal_lcdf_tail_test::expect_ad_at_defect_inputs(normal_lccdf_mix_test::fn, + normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_branch_cutoffs) { + normal_lcdf_tail_test::expect_ad_across_cutoffs(normal_lccdf_mix_test::fn, + normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_tail_derivatives) { + normal_lcdf_tail_test::expect_tail_derivatives(normal_lccdf_mix_test::fn, + normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_derivatives_are_finite) { + normal_lcdf_tail_test::expect_derivatives_finite(normal_lccdf_mix_test::fn); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lccdf_far_tail_gradient) { + normal_lcdf_tail_test::expect_far_tail_gradient(normal_lccdf_mix_test::fn, + normal_lccdf_mix_test::dir); +} + +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); } 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 7a737c87c77..b17a93685d5 100644 --- a/test/unit/math/mix/prob/normal_cdf_log_test.cpp +++ b/test/unit/math/mix/prob/normal_cdf_log_test.cpp @@ -1,19 +1,54 @@ #include #include +#include TEST_F(AgradRev, mathMixScalFun_normal_lcdf) { - auto f = [](const double mu, const double sigma) { - return [=](const auto& y) { return stan::math::normal_lcdf(y, mu, sigma); }; + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::normal_lcdf(y, mu, sigma); }; - stan::test::expect_ad(f(0.0, 1.0), -50.0); - stan::test::expect_ad(f(0.0, 1.0), -20.0 * stan::math::SQRT_TWO); - stan::test::expect_ad(f(0.0, 1.0), -5.5); - stan::test::expect_ad(f(0.0, 1.0), 0.0); - stan::test::expect_ad(f(0.0, 1.0), 0.15); - stan::test::expect_ad(f(0.0, 1.0), 1.14); - stan::test::expect_ad(f(0.0, 1.0), 3.00); - stan::test::expect_ad(f(0.0, 1.0), 10.00); - stan::test::expect_ad(f(-1.0, 2.0), 1.50); - stan::test::expect_ad(f(2.0, 1.0), 0.50); + stan::test::expect_ad(f, -50.0, 0.0, 1.0); + stan::test::expect_ad(f, -20.0 * stan::math::SQRT_TWO, 0.0, 1.0); + stan::test::expect_ad(f, -5.5, 0.0, 1.0); + stan::test::expect_ad(f, 0.0, 0.0, 1.0); + stan::test::expect_ad(f, 0.15, 0.0, 1.0); + stan::test::expect_ad(f, 1.14, 0.0, 1.0); + stan::test::expect_ad(f, 3.00, 0.0, 1.0); + stan::test::expect_ad(f, 10.00, 0.0, 1.0); + stan::test::expect_ad(f, 1.50, -1.0, 2.0); + stan::test::expect_ad(f, 0.50, 2.0, 1.0); +} + +namespace normal_lcdf_mix_test { +auto fn = [](const auto& y) { return stan::math::normal_lcdf(y, 0.0, 1.0); }; +constexpr double dir = normal_lcdf_tail_test::orientation::lcdf; +} // namespace normal_lcdf_mix_test + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_defect_inputs) { + normal_lcdf_tail_test::expect_ad_at_defect_inputs(normal_lcdf_mix_test::fn, + normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_branch_cutoffs) { + normal_lcdf_tail_test::expect_ad_across_cutoffs(normal_lcdf_mix_test::fn, + normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_tail_derivatives) { + normal_lcdf_tail_test::expect_tail_derivatives(normal_lcdf_mix_test::fn, + normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_derivatives_are_finite) { + normal_lcdf_tail_test::expect_derivatives_finite(normal_lcdf_mix_test::fn); +} + +TEST_F(AgradRev, mathMixScalFun_normal_lcdf_far_tail_gradient) { + normal_lcdf_tail_test::expect_far_tail_gradient(normal_lcdf_mix_test::fn, + normal_lcdf_mix_test::dir); +} + +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); } 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 new file mode 100644 index 00000000000..4cadb205dd5 --- /dev/null +++ b/test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp @@ -0,0 +1,260 @@ +#ifndef TEST_UNIT_MATH_MIX_PROB_NORMAL_LCDF_TAIL_TEST_HELPERS_HPP +#define TEST_UNIT_MATH_MIX_PROB_NORMAL_LCDF_TAIL_TEST_HELPERS_HPP + +#include +#include +#include +#include + +/** + * Shared checks for the tails of normal_lcdf, std_normal_lcdf and their + * reflected counterparts normal_lccdf and std_normal_lccdf. + * + * All four share one piecewise implementation, so they share one set of + * reference values. Those are given here for `normal_lcdf(y, 0, 1)`. A + * reflected function -- `lccdf(y) == lcdf(-y)` -- is exercised at `-y`, and its + * odd-order derivatives flip sign while even-order ones do not. Pass + * `orientation::lcdf` or `orientation::lccdf` to select. + * + * 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. + */ +namespace normal_lcdf_tail_test { + +/** Sign relating the function under test to normal_lcdf. */ +struct orientation { + static constexpr double lcdf = 1.0; + static constexpr double lccdf = -1.0; +}; + +/** d2/dy2 of `f` at `y`, pure forward. */ +template +double d2(const F& f, double y) { + stan::math::fvar> yv; + yv.val_.val_ = y; + yv.val_.d_ = 1.0; + yv.d_.val_ = 1.0; + return f(yv).d_.d_; +} + +/** + * d3/dy3 of `f` at `y`, pure forward. + * + * `fvar>>` is deliberate. test/unit/math/test_ad.hpp defines + * `d_t, v_t, fd_t, ffd_t, fv_t, ffv_t` and stops at `fvar>`, so + * expect_ad never instantiates this type at any input. + */ +template +double d3(const F& f, double y) { + stan::math::fvar>> yv; + yv.val_.val_.val_ = y; + yv.val_.val_.d_ = 1.0; + yv.val_.d_.val_ = 1.0; + yv.d_.val_.val_ = 1.0; + return f(yv).d_.d_.d_; +} + +/** Value and reverse-mode gradient of `f` at `y`. */ +template +void value_and_grad(const F& f, double y, double* value, double* grad) { + stan::math::var yv = y; + stan::math::var lp = f(yv); + lp.grad(); + *value = lp.val(); + *grad = yv.adj(); + stan::math::set_zero_all_adjoints(); +} + +/** + * 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. + */ +inline double mills_dlogphi_dy(double y) { + const double s = y / std::sqrt(2.0); + const double inv = 1.0 / (s * s); + const double series + = 1.0 + inv * (-0.5 + inv * (0.75 + inv * (-1.875 + inv * 6.5625))); + return -2.0 * s / series / std::sqrt(2.0); +} + +/** + * 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. + * + * 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. + */ +template +void expect_ad_across_cutoffs(const F& f, double dir) { + const double cutoffs[] = {2.9, 2.5, 2.1, 1.5, 0.8, 0.1, 0.0, + -2.1, -3.9, -4.0, -7.0, -17.0, -29.0}; + for (double cut : cutoffs) { + stan::test::expect_ad(f, dir * (cut - 0.01) * stan::math::SQRT_TWO); + stan::test::expect_ad(f, dir * (cut + 0.01) * stan::math::SQRT_TWO); + } +} + +/** + * expect_ad at the inputs that exposed each of the three tail defects: the + * `exp(x2)` overflow above `scaled_diff = 26.6`, the `erfc^2` underflow window + * `scaled_diff in (-20, -19.2103)`, and the far-tail residual fit that + * diverges from the linear Mills asymptote. + */ +template +void expect_ad_at_defect_inputs(const F& f, double dir) { + const double inputs[] = {37.6, 40.0, 50.0, -27.5, -28.0, -60.0, -100.0}; + for (double y : inputs) { + stan::test::expect_ad(f, dir * y); + } +} + +/** + * 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. + */ +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}}; + + 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) + << "derivative branch drifted at y = " << y; + } +} + +/** + * Second and third derivatives in the tail against 60-digit references. + * + * Asserted directly rather than through expect_ad because expect_ad cannot see + * either failure. The order-2 rows were silently zero before the `exp(-x2)` + * rearrangement, and a true value of ~1e-158 is invisible to a + * finite-difference comparison whose relative tolerance floors its denominator + * at 1. The order-3 rows use a type expect_ad never instantiates, and fail + * earlier than anything it does build -- silently 0 at 20, NaN from 22, + * against 37.6 for order 2. + */ +template +void expect_tail_derivatives(const F& f, double dir) { + struct tail_ref { + double y; + double expected; + }; + const tail_ref d2_cases[] = {{26.0, -1.67628759256344147e-146}, + {27.0, -5.39430099975435425e-158}, + {30.0, -4.42093840463564223e-195}, + {33.0, -4.42668490225313907e-236}, + {37.0, -7.84402424064104056e-297}}; + for (const tail_ref& c : d2_cases) { + const double y = dir * c.y; + const double got = d2(f, y); + ASSERT_FALSE(std::isnan(got)) << "NaN d2 at y = " << y; + EXPECT_NE(0.0, got) << "d2 collapsed to zero at y = " << y; + EXPECT_LT(std::fabs(got / c.expected - 1.0), 1e-6) << "d2 at y = " << y; + } + + const tail_ref d3_cases[] = {{20.0, 2.20285839650174548e-85}, + {22.0, 1.53317799001460164e-103}, + {25.0, 4.77605215552570066e-134}, + {30.0, 1.32480787525581414e-193}}; + for (const tail_ref& c : d3_cases) { + const double y = dir * c.y; + const double got = d3(f, y); + ASSERT_FALSE(std::isnan(got)) << "NaN d3 at y = " << y; + EXPECT_NE(0.0, got) << "d3 collapsed to zero at y = " << y; + EXPECT_LT(std::fabs(got / (dir * c.expected) - 1.0), 1e-6) + << "d3 at y = " << y; + } +} + +/** + * Sweep for non-finite derivatives across the whole useful range. Cheap at 0.5 + * spacing for two direct fvar evaluations, where expect_ad -- six AD modes plus + * finite differencing per point -- would not be. This single loop would have + * caught every tail defect in this function's history at once. + */ +template +void expect_derivatives_finite(const F& f) { + for (double y = -120.0; y <= 120.0; y += 0.5) { + EXPECT_FALSE(std::isnan(d2(f, y))) << "d2 NaN at y = " << y; + EXPECT_FALSE(std::isnan(d3(f, y))) << "d3 NaN at y = " << y; + } +} + +/** + * 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. + */ +template +void expect_far_tail_gradient(const F& f, double dir) { + struct far_ref { + double y; + double value; + double grad; + }; + const far_ref cases[] + = {{-60.0, -1.80501356068056725e+03, 6.00166574202411240e+01}, + {-100.0, -5.00552420869420530e+03, 1.00009998000999261e+02}, + {-300.0, -4.50066227321186598e+04, 3.00003333259263400e+02}}; + for (const far_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(value / c.value - 1.0), 1e-12) << "value at y = " << y; + EXPECT_LT(std::fabs(grad / (dir * c.grad) - 1.0), 1e-10) + << "gradient at y = " << y; + } + + for (double ly = -300.0; ly <= -42.0; ly += 0.25) { + const double y = dir * ly; + double value = 0; + double grad = 0; + value_and_grad(f, y, &value, &grad); + EXPECT_LT(std::fabs(grad / (dir * mills_dlogphi_dy(ly)) - 1.0), 1e-9) + << "gradient drifts from the asymptotic Mills ratio at y = " << y; + } +} + +} // namespace normal_lcdf_tail_test + +#endif 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 3c2b01d395a..e052be64810 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 @@ -1,15 +1,51 @@ #include #include +#include TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf) { auto f = [](const auto& y) { return stan::math::std_normal_lccdf(y); }; - stan::test::expect_ad(f, -50.0); - stan::test::expect_ad(f, -20.0 * stan::math::SQRT_TWO); - stan::test::expect_ad(f, -5.5); + stan::test::expect_ad(f, 50.0); + stan::test::expect_ad(f, 20.0 * stan::math::SQRT_TWO); + stan::test::expect_ad(f, 5.5); stan::test::expect_ad(f, 0.0); - stan::test::expect_ad(f, 0.15); - stan::test::expect_ad(f, 1.14); - stan::test::expect_ad(f, 3.00); - stan::test::expect_ad(f, 10.00); + stan::test::expect_ad(f, -0.15); + stan::test::expect_ad(f, -1.14); + stan::test::expect_ad(f, -3.00); + stan::test::expect_ad(f, -10.00); +} + +namespace std_normal_lccdf_mix_test { +auto fn = [](const auto& y) { return stan::math::std_normal_lccdf(y); }; +constexpr double dir = normal_lcdf_tail_test::orientation::lccdf; +} // namespace std_normal_lccdf_mix_test + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_defect_inputs) { + normal_lcdf_tail_test::expect_ad_at_defect_inputs( + std_normal_lccdf_mix_test::fn, std_normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_branch_cutoffs) { + normal_lcdf_tail_test::expect_ad_across_cutoffs( + std_normal_lccdf_mix_test::fn, std_normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_tail_derivatives) { + normal_lcdf_tail_test::expect_tail_derivatives( + std_normal_lccdf_mix_test::fn, std_normal_lccdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_derivatives_are_finite) { + normal_lcdf_tail_test::expect_derivatives_finite( + std_normal_lccdf_mix_test::fn); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lccdf_far_tail_gradient) { + normal_lcdf_tail_test::expect_far_tail_gradient( + std_normal_lccdf_mix_test::fn, std_normal_lccdf_mix_test::dir); +} + +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); } diff --git a/test/unit/math/mix/prob/std_normal_cdf_log_test.cpp b/test/unit/math/mix/prob/std_normal_cdf_log_test.cpp index 7f566f3452f..50eb84baec1 100644 --- a/test/unit/math/mix/prob/std_normal_cdf_log_test.cpp +++ b/test/unit/math/mix/prob/std_normal_cdf_log_test.cpp @@ -1,5 +1,6 @@ #include #include +#include TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf) { auto f = [](const auto& y) { return stan::math::std_normal_lcdf(y); }; @@ -13,3 +14,38 @@ TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf) { stan::test::expect_ad(f, 3.00); stan::test::expect_ad(f, 10.00); } + +namespace std_normal_lcdf_mix_test { +auto fn = [](const auto& y) { return stan::math::std_normal_lcdf(y); }; +constexpr double dir = normal_lcdf_tail_test::orientation::lcdf; +} // namespace std_normal_lcdf_mix_test + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_defect_inputs) { + normal_lcdf_tail_test::expect_ad_at_defect_inputs( + std_normal_lcdf_mix_test::fn, std_normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_branch_cutoffs) { + normal_lcdf_tail_test::expect_ad_across_cutoffs( + std_normal_lcdf_mix_test::fn, std_normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_tail_derivatives) { + normal_lcdf_tail_test::expect_tail_derivatives(std_normal_lcdf_mix_test::fn, + std_normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_derivatives_are_finite) { + normal_lcdf_tail_test::expect_derivatives_finite( + std_normal_lcdf_mix_test::fn); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_far_tail_gradient) { + normal_lcdf_tail_test::expect_far_tail_gradient( + std_normal_lcdf_mix_test::fn, std_normal_lcdf_mix_test::dir); +} + +TEST_F(AgradRev, mathMixScalFun_std_normal_lcdf_branch_accuracy) { + normal_lcdf_tail_test::expect_branch_accuracy(std_normal_lcdf_mix_test::fn, + std_normal_lcdf_mix_test::dir); +} 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 cc17ed90eab..5bab9655a32 100644 --- a/test/unit/math/prim/prob/normal_ccdf_log_test.cpp +++ b/test/unit/math/prim/prob/normal_ccdf_log_test.cpp @@ -16,34 +16,325 @@ TEST(ProbNormal, ccdf_log_matches_lccdf) { TEST(ProbNormal, lccdf_tail) { using stan::math::normal_lccdf; - EXPECT_FLOAT_EQ(-6.661338147750941214694e-16, normal_lccdf(-8.0, 0, 1)); - EXPECT_FLOAT_EQ(-3.186340080674249758114e-14, normal_lccdf(-7.5, 0, 1)); - EXPECT_FLOAT_EQ(-1.279865102788699562477e-12, normal_lccdf(-7.0, 0, 1)); - EXPECT_FLOAT_EQ(-4.015998644826973564545e-11, normal_lccdf(-6.5, 0, 1)); - EXPECT_FLOAT_EQ(-9.865877009111571184118e-10, normal_lccdf(-6.0, 0, 1)); - EXPECT_FLOAT_EQ(-1.898956265833514866414e-08, normal_lccdf(-5.5, 0, 1)); - EXPECT_FLOAT_EQ(-2.866516130081049047962e-07, normal_lccdf(-5.0, 0, 1)); - EXPECT_FLOAT_EQ(-3.397678896843115195074e-06, normal_lccdf(-4.5, 0, 1)); - EXPECT_FLOAT_EQ(-3.167174337748932124543e-05, normal_lccdf(-4.0, 0, 1)); - EXPECT_FLOAT_EQ(-0.0002326561413768195969113, normal_lccdf(-3.5, 0, 1)); - EXPECT_FLOAT_EQ(-0.001350809964748202673598, normal_lccdf(-3.0, 0, 1)); - EXPECT_FLOAT_EQ(-0.0062290254858600267035, normal_lccdf(-2.5, 0, 1)); - EXPECT_FLOAT_EQ(-0.02301290932896348992442, normal_lccdf(-2.0, 0, 1)); - EXPECT_FLOAT_EQ(-0.06914345561223400604689, normal_lccdf(-1.5, 0, 1)); - EXPECT_FLOAT_EQ(-0.1727537790234499048836, normal_lccdf(-1.0, 0, 1)); - EXPECT_FLOAT_EQ(-0.3689464152886565151412, normal_lccdf(-0.5, 0, 1)); - EXPECT_FLOAT_EQ(-0.6931471805599452862268, normal_lccdf(0, 0, 1)); - EXPECT_FLOAT_EQ(-1.175911761593618320987, normal_lccdf(0.5, 0, 1)); - EXPECT_FLOAT_EQ(-1.841021645009263352222, normal_lccdf(1.0, 0, 1)); - EXPECT_FLOAT_EQ(-2.705944400823889317564, normal_lccdf(1.5, 0, 1)); - EXPECT_FLOAT_EQ(-3.78318433368203210776, normal_lccdf(2.0, 0, 1)); - EXPECT_FLOAT_EQ(-5.081648277278686620662, normal_lccdf(2.5, 0, 1)); - EXPECT_FLOAT_EQ(-6.607726221510342945464, normal_lccdf(3.0, 0, 1)); - EXPECT_FLOAT_EQ(-8.366065308344028395027, normal_lccdf(3.5, 0, 1)); - EXPECT_FLOAT_EQ(-10.36010148652728979357, normal_lccdf(4.0, 0, 1)); - EXPECT_FLOAT_EQ(-12.59241973571053385683, normal_lccdf(4.5, 0, 1)); - EXPECT_FLOAT_EQ(-15.06499839383403838156, normal_lccdf(5.0, 0, 1)); - EXPECT_FLOAT_EQ(-17.77937635198566113104, normal_lccdf(5.5, 0, 1)); - EXPECT_FLOAT_EQ(-20.73676889383495947072, normal_lccdf(6.0, 0, 1)); - EXPECT_FLOAT_EQ(-23.93814997800869548428, normal_lccdf(6.5, 0, 1)); + // The test values come from R 4.6.1 and cover the expected useful range of + // the function. When z <= -38.5, even the log of the CCDF is + // indistinguishable from 0.0 in double precision. + // + // q <- + // c( + // seq(-38, -11, by = 1.0), + // seq(-10, 10, by = 0.5), + // seq(11, 50, by = 1.0), + // 10^seq(2, 8, by = 1) + // ) + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, normal_lccdf(%.17g, 0, 1));\n", + // pnorm(q[i], lower.tail = FALSE, log.p = TRUE), + // q[i] + // ) + // ) + // } + + EXPECT_FLOAT_EQ(-2.8854283510039645e-316, normal_lccdf(-38, 0, 1)); + EXPECT_FLOAT_EQ(-5.7255712225245771e-300, normal_lccdf(-37, 0, 1)); + EXPECT_FLOAT_EQ(-4.1826240657972830e-284, normal_lccdf(-36, 0, 1)); + EXPECT_FLOAT_EQ(-1.1249107064724062e-268, normal_lccdf(-35, 0, 1)); + EXPECT_FLOAT_EQ(-1.1138987855743795e-253, normal_lccdf(-34, 0, 1)); + EXPECT_FLOAT_EQ(-4.0611856209158557e-239, normal_lccdf(-33, 0, 1)); + EXPECT_FLOAT_EQ(-5.4520806035123956e-225, normal_lccdf(-32, 0, 1)); + EXPECT_FLOAT_EQ(-2.6952500812005002e-211, normal_lccdf(-31, 0, 1)); + EXPECT_FLOAT_EQ(-4.9067139271481872e-198, normal_lccdf(-30, 0, 1)); + EXPECT_FLOAT_EQ(-3.2897852667043802e-185, normal_lccdf(-29, 0, 1)); + EXPECT_FLOAT_EQ(-8.1238694696594273e-173, normal_lccdf(-28, 0, 1)); + EXPECT_FLOAT_EQ(-7.3894810068850200e-161, normal_lccdf(-27, 0, 1)); + EXPECT_FLOAT_EQ(-2.4760633155033892e-149, normal_lccdf(-26, 0, 1)); + EXPECT_FLOAT_EQ(-3.0566967063825616e-138, normal_lccdf(-25, 0, 1)); + EXPECT_FLOAT_EQ(-1.3903921185497032e-127, normal_lccdf(-24, 0, 1)); + EXPECT_FLOAT_EQ(-2.3306370062206492e-117, normal_lccdf(-23, 0, 1)); + EXPECT_FLOAT_EQ(-1.4398924351450790e-107, normal_lccdf(-22, 0, 1)); + EXPECT_FLOAT_EQ(-3.2792780189790367e-98, normal_lccdf(-21, 0, 1)); + EXPECT_FLOAT_EQ(-2.7536241186062337e-89, normal_lccdf(-20, 0, 1)); + EXPECT_FLOAT_EQ(-8.5272239526309772e-81, normal_lccdf(-19, 0, 1)); + EXPECT_FLOAT_EQ(-9.7409489189371508e-73, normal_lccdf(-18, 0, 1)); + EXPECT_FLOAT_EQ(-4.1059962020989074e-65, normal_lccdf(-17, 0, 1)); + EXPECT_FLOAT_EQ(-6.3887544005380882e-58, normal_lccdf(-16, 0, 1)); + EXPECT_FLOAT_EQ(-3.6709661993127514e-51, normal_lccdf(-15, 0, 1)); + EXPECT_FLOAT_EQ(-7.7935368191928000e-45, normal_lccdf(-14, 0, 1)); + EXPECT_FLOAT_EQ(-6.1171643995498803e-39, normal_lccdf(-13, 0, 1)); + EXPECT_FLOAT_EQ(-1.7764821120776790e-33, normal_lccdf(-12, 0, 1)); + EXPECT_FLOAT_EQ(-1.9106595744986757e-28, normal_lccdf(-11, 0, 1)); + EXPECT_FLOAT_EQ(-7.6198530241605269e-24, normal_lccdf(-10, 0, 1)); + EXPECT_FLOAT_EQ(-1.0494515075362608e-21, normal_lccdf(-9.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.1285884059538408e-19, normal_lccdf(-9, 0, 1)); + EXPECT_FLOAT_EQ(-9.4795348222033192e-18, normal_lccdf(-8.5, 0, 1)); + EXPECT_FLOAT_EQ(-6.2209605742717868e-16, normal_lccdf(-8, 0, 1)); + EXPECT_FLOAT_EQ(-3.1908916729109475e-14, normal_lccdf(-7.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.2798125438866541e-12, normal_lccdf(-7, 0, 1)); + EXPECT_FLOAT_EQ(-4.0160005839397589e-11, normal_lccdf(-6.5, 0, 1)); + EXPECT_FLOAT_EQ(-9.8658764552437559e-10, normal_lccdf(-6, 0, 1)); + EXPECT_FLOAT_EQ(-1.8989562646189464e-08, normal_lccdf(-5.5, 0, 1)); + EXPECT_FLOAT_EQ(-2.8665161296376358e-07, normal_lccdf(-5, 0, 1)); + EXPECT_FLOAT_EQ(-3.3976788968344657e-06, normal_lccdf(-4.5, 0, 1)); + EXPECT_FLOAT_EQ(-3.1671743377489267e-05, normal_lccdf(-4, 0, 1)); + EXPECT_FLOAT_EQ(-0.00023265614137680445, normal_lccdf(-3.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.0013508099647481938, normal_lccdf(-3, 0, 1)); + EXPECT_FLOAT_EQ(-0.0062290254858600024, normal_lccdf(-2.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.023012909328963493, normal_lccdf(-2, 0, 1)); + EXPECT_FLOAT_EQ(-0.069143455612233992, normal_lccdf(-1.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.17275377902344988, normal_lccdf(-1, 0, 1)); + EXPECT_FLOAT_EQ(-0.36894641528865652, normal_lccdf(-0.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.69314718055994529, normal_lccdf(0, 0, 1)); + EXPECT_FLOAT_EQ(-1.1759117615936185, normal_lccdf(0.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.8410216450092636, normal_lccdf(1, 0, 1)); + EXPECT_FLOAT_EQ(-2.7059444008238898, normal_lccdf(1.5, 0, 1)); + EXPECT_FLOAT_EQ(-3.7831843336820317, normal_lccdf(2, 0, 1)); + EXPECT_FLOAT_EQ(-5.0816482772786902, normal_lccdf(2.5, 0, 1)); + EXPECT_FLOAT_EQ(-6.6077262215103492, normal_lccdf(3, 0, 1)); + EXPECT_FLOAT_EQ(-8.3660653083440941, normal_lccdf(3.5, 0, 1)); + EXPECT_FLOAT_EQ(-10.360101486527292, normal_lccdf(4, 0, 1)); + EXPECT_FLOAT_EQ(-12.592419735713079, normal_lccdf(4.5, 0, 1)); + EXPECT_FLOAT_EQ(-15.064998393988725, normal_lccdf(5, 0, 1)); + EXPECT_FLOAT_EQ(-17.779376352625260, normal_lccdf(5.5, 0, 1)); + EXPECT_FLOAT_EQ(-20.736768949974707, normal_lccdf(6, 0, 1)); + EXPECT_FLOAT_EQ(-23.938149495161838, normal_lccdf(6.5, 0, 1)); + EXPECT_FLOAT_EQ(-27.384307498811076, normal_lccdf(7, 0, 1)); + EXPECT_FLOAT_EQ(-31.075890902890002, normal_lccdf(7.5, 0, 1)); + EXPECT_FLOAT_EQ(-35.013437159914552, normal_lccdf(8, 0, 1)); + EXPECT_FLOAT_EQ(-39.197396428217672, normal_lccdf(8.5, 0, 1)); + EXPECT_FLOAT_EQ(-43.628149113332114, normal_lccdf(9, 0, 1)); + EXPECT_FLOAT_EQ(-48.306019298965232, normal_lccdf(9.5, 0, 1)); + EXPECT_FLOAT_EQ(-53.231285150512470, normal_lccdf(10, 0, 1)); + EXPECT_FLOAT_EQ(-63.824934094423718, normal_lccdf(11, 0, 1)); + EXPECT_FLOAT_EQ(-75.410673001568796, normal_lccdf(12, 0, 1)); + EXPECT_FLOAT_EQ(-87.989719971022524, normal_lccdf(13, 0, 1)); + EXPECT_FLOAT_EQ(-101.56303440744996, normal_lccdf(14, 0, 1)); + EXPECT_FLOAT_EQ(-116.13138484571169, normal_lccdf(15, 0, 1)); + EXPECT_FLOAT_EQ(-131.69539607375970, normal_lccdf(16, 0, 1)); + EXPECT_FLOAT_EQ(-148.25558265098039, normal_lccdf(17, 0, 1)); + EXPECT_FLOAT_EQ(-165.81237325071419, normal_lccdf(18, 0, 1)); + EXPECT_FLOAT_EQ(-184.36612866916096, normal_lccdf(19, 0, 1)); + EXPECT_FLOAT_EQ(-203.91715537109727, normal_lccdf(20, 0, 1)); + EXPECT_FLOAT_EQ(-224.46571583141449, normal_lccdf(21, 0, 1)); + EXPECT_FLOAT_EQ(-246.01203653738091, normal_lccdf(22, 0, 1)); + EXPECT_FLOAT_EQ(-268.55631425686312, normal_lccdf(23, 0, 1)); + EXPECT_FLOAT_EQ(-292.09872100320780, normal_lccdf(24, 0, 1)); + EXPECT_FLOAT_EQ(-316.63940800802027, normal_lccdf(25, 0, 1)); + EXPECT_FLOAT_EQ(-342.17850892992783, normal_lccdf(26, 0, 1)); + EXPECT_FLOAT_EQ(-368.71614246865636, normal_lccdf(27, 0, 1)); + EXPECT_FLOAT_EQ(-396.25241451163106, normal_lccdf(28, 0, 1)); + EXPECT_FLOAT_EQ(-424.78741990973015, normal_lccdf(29, 0, 1)); + EXPECT_FLOAT_EQ(-454.32124395634321, normal_lccdf(30, 0, 1)); + EXPECT_FLOAT_EQ(-484.85396362717927, normal_lccdf(31, 0, 1)); + EXPECT_FLOAT_EQ(-516.38564862572537, normal_lccdf(32, 0, 1)); + EXPECT_FLOAT_EQ(-548.91636226973810, normal_lccdf(33, 0, 1)); + EXPECT_FLOAT_EQ(-582.44616224687172, normal_lccdf(34, 0, 1)); + EXPECT_FLOAT_EQ(-616.97510126192253, normal_lccdf(35, 0, 1)); + EXPECT_FLOAT_EQ(-652.50322759379844, normal_lccdf(36, 0, 1)); + EXPECT_FLOAT_EQ(-689.03058557689064, normal_lccdf(37, 0, 1)); + EXPECT_FLOAT_EQ(-726.55721601882010, normal_lccdf(38, 0, 1)); + EXPECT_FLOAT_EQ(-765.08315656437753, normal_lccdf(39, 0, 1)); + EXPECT_FLOAT_EQ(-804.60844201375380, normal_lccdf(40, 0, 1)); + EXPECT_FLOAT_EQ(-845.13310460177456, normal_lccdf(41, 0, 1)); + EXPECT_FLOAT_EQ(-886.65717424372951, normal_lccdf(42, 0, 1)); + EXPECT_FLOAT_EQ(-929.18067875247391, normal_lccdf(43, 0, 1)); + EXPECT_FLOAT_EQ(-972.70364403073665, normal_lccdf(44, 0, 1)); + EXPECT_FLOAT_EQ(-1017.2260942419524, normal_lccdf(45, 0, 1)); + EXPECT_FLOAT_EQ(-1062.7480519624305, normal_lccdf(46, 0, 1)); + EXPECT_FLOAT_EQ(-1109.2695383172531, normal_lccdf(47, 0, 1)); + EXPECT_FLOAT_EQ(-1156.7905731019453, normal_lccdf(48, 0, 1)); + EXPECT_FLOAT_EQ(-1205.3111748916654, normal_lccdf(49, 0, 1)); + EXPECT_FLOAT_EQ(-1254.8313611394199, normal_lccdf(50, 0, 1)); + EXPECT_FLOAT_EQ(-5005.5242086942053, normal_lccdf(100, 0, 1)); + EXPECT_FLOAT_EQ(-500007.82669481216, normal_lccdf(1000, 0, 1)); + EXPECT_FLOAT_EQ(-50000010.129278913, normal_lccdf(10000, 0, 1)); + EXPECT_FLOAT_EQ(-5000000012.4318638, normal_lccdf(100000, 0, 1)); + EXPECT_FLOAT_EQ(-500000000014.73444, normal_lccdf(1000000, 0, 1)); + EXPECT_FLOAT_EQ(-50000000000017.039, normal_lccdf(10000000, 0, 1)); + EXPECT_FLOAT_EQ(-5000000000000019.0, normal_lccdf(100000000, 0, 1)); +} + +TEST(ProbNormal, lccdf_with_mu_and_sigma) { + using stan::math::normal_lccdf; + + // The test values come from R 4.6.1 and check behaviour around the branch + // points in `normal_lcdf` when mu != 0 and sigma != 1. + // + // parameter_sets <- list( + // c(mu = 2, sigma = 3), + // c(mu = -50, sigma = 20), + // c(mu = 0.5, sigma = 0.01) + // ) + // branch_points <- c( + // 2.9, + // 2.5, + // 2.1, + // 1.5, + // 0.8, + // 0.1, + // 0.0, + // -2.1, + // -3.9, + // -4.0, + // -7.0, + // -17.0, + // -29.0 + // ) + // test_vals <- sort(-c(branch_points - 0.01, branch_points + 0.01)) + // for (pars in parameter_sets) { + // mu <- pars[["mu"]] + // sigma <- pars[["sigma"]] + // q <- mu + sigma * sqrt(2) * test_vals + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, normal_lccdf(%.17g, %.17g, %.17g));\n", + // pnorm(q[i], mean = mu, sd = sigma, lower.tail = FALSE, log.p = + // TRUE), q[i], mu, sigma + // ) + // ) + // } + // } + + EXPECT_FLOAT_EQ(-1.9328842921129592e-05, + normal_lccdf(-10.346084399517119, 2, 3)); + EXPECT_FLOAT_EQ(-2.1842328461060815e-05, + normal_lccdf(-10.261231585774736, 2, 3)); + EXPECT_FLOAT_EQ(-0.00019287133934318566, + normal_lccdf(-8.6490281246694067, 2, 3)); + EXPECT_FLOAT_EQ(-0.00021466697317301965, + normal_lccdf(-8.5641753109270216, 2, 3)); + EXPECT_FLOAT_EQ(-0.0014235903084217206, + normal_lccdf(-6.9519718498216925, 2, 3)); + EXPECT_FLOAT_EQ(-0.0015609874450218889, + normal_lccdf(-6.8671190360793091, 2, 3)); + EXPECT_FLOAT_EQ(-0.016496955508933941, + normal_lccdf(-4.4063874375501211, 2, 3)); + EXPECT_FLOAT_EQ(-0.017706913768612061, + normal_lccdf(-4.3215346238077359, 2, 3)); + EXPECT_FLOAT_EQ(-0.13467302655617183, + normal_lccdf(-1.4365389565666216, 2, 3)); + EXPECT_FLOAT_EQ(-0.14150397288270808, + normal_lccdf(-1.3516861428242359, 2, 3)); + EXPECT_FLOAT_EQ(-0.57658898698359096, normal_lccdf(1.5333095244168786, 2, 3)); + EXPECT_FLOAT_EQ(-0.59667350159732879, normal_lccdf(1.6181623381592642, 2, 3)); + EXPECT_FLOAT_EQ(-0.68192694790236175, normal_lccdf(1.9575735931288072, 2, 3)); + EXPECT_FLOAT_EQ(-0.70449473678943686, normal_lccdf(2.042426406871193, 2, 3)); + EXPECT_FLOAT_EQ(-6.4632170725824736, normal_lccdf(10.867119036079309, 2, 3)); + EXPECT_FLOAT_EQ(-6.5552849228594905, normal_lccdf(10.951971849821692, 2, 3)); + EXPECT_FLOAT_EQ(-17.786697578589479, normal_lccdf(18.503872272894021, 2, 3)); + EXPECT_FLOAT_EQ(-17.947533669545574, normal_lccdf(18.588725086636405, 2, 3)); + EXPECT_FLOAT_EQ(-18.598659064395832, normal_lccdf(18.928136341605949, 2, 3)); + EXPECT_FLOAT_EQ(-18.763386609836139, normal_lccdf(19.012989155348336, 2, 3)); + EXPECT_FLOAT_EQ(-52.080076508792267, normal_lccdf(31.656058402963808, 2, 3)); + EXPECT_FLOAT_EQ(-52.362878114253100, normal_lccdf(31.740911216706191, 2, 3)); + EXPECT_FLOAT_EQ(-292.75996176272099, normal_lccdf(74.082465274156661, 2, 3)); + EXPECT_FLOAT_EQ(-293.44113419738113, normal_lccdf(74.167318087899062, 2, 3)); + EXPECT_FLOAT_EQ(-845.05315712467279, normal_lccdf(124.99415351958808, 2, 3)); + EXPECT_FLOAT_EQ(-846.21384596225801, normal_lccdf(125.07900633333048, 2, 3)); + EXPECT_FLOAT_EQ(-1.9328842921129592e-05, + normal_lccdf(-132.30722933011413, -50, 20)); + EXPECT_FLOAT_EQ(-2.1842328461060734e-05, + normal_lccdf(-131.74154390516492, -50, 20)); + EXPECT_FLOAT_EQ(-0.00019287133934318566, + normal_lccdf(-120.99352083112937, -50, 20)); + EXPECT_FLOAT_EQ(-0.00021466697317302003, + normal_lccdf(-120.42783540618014, -50, 20)); + EXPECT_FLOAT_EQ(-0.0014235903084217224, + normal_lccdf(-109.67981233214461, -50, 20)); + EXPECT_FLOAT_EQ(-0.0015609874450218915, + normal_lccdf(-109.11412690719538, -50, 20)); + EXPECT_FLOAT_EQ(-0.016496955508933965, + normal_lccdf(-92.709249583667471, -50, 20)); + EXPECT_FLOAT_EQ(-0.017706913768612061, + normal_lccdf(-92.143564158718235, -50, 20)); + EXPECT_FLOAT_EQ(-0.13467302655617194, + normal_lccdf(-72.910259710444137, -50, 20)); + EXPECT_FLOAT_EQ(-0.14150397288270816, + normal_lccdf(-72.3445742854949, -50, 20)); + EXPECT_FLOAT_EQ(-0.57658898698359096, + normal_lccdf(-53.11126983722081, -50, 20)); + EXPECT_FLOAT_EQ(-0.59667350159732879, + normal_lccdf(-52.545584412271573, -50, 20)); + EXPECT_FLOAT_EQ(-0.68192694790236175, + normal_lccdf(-50.282842712474618, -50, 20)); + EXPECT_FLOAT_EQ(-0.70449473678943675, + normal_lccdf(-49.717157287525382, -50, 20)); + EXPECT_FLOAT_EQ(-6.4632170725824709, + normal_lccdf(9.1141269071953843, -50, 20)); + EXPECT_FLOAT_EQ(-6.5552849228594887, + normal_lccdf(9.679812332144607, -50, 20)); + EXPECT_FLOAT_EQ(-17.786697578589472, + normal_lccdf(60.025815152626805, -50, 20)); + EXPECT_FLOAT_EQ(-17.947533669545571, + normal_lccdf(60.591500577576028, -50, 20)); + EXPECT_FLOAT_EQ(-18.598659064395825, + normal_lccdf(62.85424227737299, -50, 20)); + EXPECT_FLOAT_EQ(-18.763386609836136, + normal_lccdf(63.419927702322227, -50, 20)); + EXPECT_FLOAT_EQ(-52.080076508792246, + normal_lccdf(147.70705601975871, -50, 20)); + EXPECT_FLOAT_EQ(-52.362878114253085, + normal_lccdf(148.27274144470792, -50, 20)); + EXPECT_FLOAT_EQ(-292.75996176272093, + normal_lccdf(430.54976849437764, -50, 20)); + EXPECT_FLOAT_EQ(-293.44113419738102, + normal_lccdf(431.11545391932702, -50, 20)); + EXPECT_FLOAT_EQ(-845.05315712467279, + normal_lccdf(769.96102346392047, -50, 20)); + EXPECT_FLOAT_EQ(-846.21384596225778, + normal_lccdf(770.52670888886985, -50, 20)); + EXPECT_FLOAT_EQ(-1.9328842921129514e-05, + normal_lccdf(0.45884638533494293, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-2.1842328461060476e-05, + normal_lccdf(0.45912922804741751, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.00019287133934318404, + normal_lccdf(0.46450323958443529, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.00021466697317302003, + normal_lccdf(0.46478608229690993, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.0014235903084217287, + normal_lccdf(0.47016009383392771, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.0015609874450218850, + normal_lccdf(0.47044293654640229, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.016496955508933885, + normal_lccdf(0.47864537520816625, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.017706913768612078, + normal_lccdf(0.47892821792064089, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.13467302655617194, + normal_lccdf(0.48854487014477793, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.14150397288270869, + normal_lccdf(0.48882771285725257, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.57658898698359251, + normal_lccdf(0.49844436508138962, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.59667350159732790, + normal_lccdf(0.4987272077938642, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.68192694790236308, + normal_lccdf(0.49985857864376271, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.70449473678943542, + normal_lccdf(0.50014142135623729, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-6.4632170725824585, + normal_lccdf(0.52955706345359765, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-6.5552849228595038, + normal_lccdf(0.52983990616607235, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-17.786697578589461, + normal_lccdf(0.55501290757631339, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-17.947533669545543, + normal_lccdf(0.55529575028878797, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-18.598659064395846, + normal_lccdf(0.55642712113868653, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-18.763386609836136, + normal_lccdf(0.55670996385116112, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-52.080076508792246, + normal_lccdf(0.59885352800987934, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-52.362878114253050, + normal_lccdf(0.59913637072235393, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-292.75996176272110, + normal_lccdf(0.7402748842471889, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-293.44113419738096, + normal_lccdf(0.74055772695966349, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-845.05315712467279, + normal_lccdf(0.90998051173196026, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-846.21384596225801, + normal_lccdf(0.91026335444443496, 0.5, 0.01)); } diff --git a/test/unit/math/prim/prob/normal_cdf_log_test.cpp b/test/unit/math/prim/prob/normal_cdf_log_test.cpp index 2cc71187fb1..fe8ab52090b 100644 --- a/test/unit/math/prim/prob/normal_cdf_log_test.cpp +++ b/test/unit/math/prim/prob/normal_cdf_log_test.cpp @@ -10,108 +10,397 @@ TEST(ProbNormal, cdf_log_matches_lcdf) { EXPECT_FLOAT_EQ((stan::math::normal_lcdf(y, mu, sigma)), (stan::math::normal_cdf_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::normal_lcdf(y, mu, sigma)), + (stan::math::normal_lcdf(y, mu, sigma)), (stan::math::normal_cdf_log(y, mu, sigma))); } TEST(ProbNormal, lcdf_tails) { using stan::math::normal_lcdf; - using std::exp; - EXPECT_FLOAT_EQ(4.60535300958196e-308, exp(normal_lcdf(-37.5, 0, 1))); - EXPECT_FLOAT_EQ(5.72557122252458e-300, exp(normal_lcdf(-37, 0, 1))); - EXPECT_FLOAT_EQ(5.54472571307484e-292, exp(normal_lcdf(-36.5, 0, 1))); - EXPECT_FLOAT_EQ(4.18262406579728e-284, exp(normal_lcdf(-36, 0, 1))); - EXPECT_FLOAT_EQ(2.45769154066194e-276, exp(normal_lcdf(-35.5, 0, 1))); - EXPECT_FLOAT_EQ(1.12491070647241e-268, exp(normal_lcdf(-35, 0, 1))); - EXPECT_FLOAT_EQ(4.01072896657726e-261, exp(normal_lcdf(-34.5, 0, 1))); - EXPECT_FLOAT_EQ(1.11389878557438e-253, exp(normal_lcdf(-34, 0, 1))); - EXPECT_FLOAT_EQ(2.40983869512039e-246, exp(normal_lcdf(-33.5, 0, 1))); - EXPECT_FLOAT_EQ(4.06118562091586e-239, exp(normal_lcdf(-33, 0, 1))); - EXPECT_FLOAT_EQ(5.33142435967881e-232, exp(normal_lcdf(-32.5, 0, 1))); - EXPECT_FLOAT_EQ(5.4520806035124e-225, exp(normal_lcdf(-32, 0, 1))); - EXPECT_FLOAT_EQ(4.34323260103177e-218, exp(normal_lcdf(-31.5, 0, 1))); - EXPECT_FLOAT_EQ(2.6952500812005e-211, exp(normal_lcdf(-31, 0, 1))); - EXPECT_FLOAT_EQ(1.30293791317808e-204, exp(normal_lcdf(-30.5, 0, 1))); - EXPECT_FLOAT_EQ(4.90671392714819e-198, exp(normal_lcdf(-30, 0, 1))); - EXPECT_FLOAT_EQ(1.43947455222918e-191, exp(normal_lcdf(-29.5, 0, 1))); - EXPECT_FLOAT_EQ(3.28978526670438e-185, exp(normal_lcdf(-29, 0, 1))); - EXPECT_FLOAT_EQ(5.85714125380634e-179, exp(normal_lcdf(-28.5, 0, 1))); - EXPECT_FLOAT_EQ(8.12386946965943e-173, exp(normal_lcdf(-28, 0, 1))); - EXPECT_FLOAT_EQ(8.77817055687808e-167, exp(normal_lcdf(-27.5, 0, 1))); - EXPECT_FLOAT_EQ(7.38948100688502e-161, exp(normal_lcdf(-27, 0, 1))); - EXPECT_FLOAT_EQ(4.84616266030332e-155, exp(normal_lcdf(-26.5, 0, 1))); - EXPECT_FLOAT_EQ(2.47606331550339e-149, exp(normal_lcdf(-26, 0, 1))); - EXPECT_FLOAT_EQ(9.85623651896393e-144, exp(normal_lcdf(-25.5, 0, 1))); - EXPECT_FLOAT_EQ(3.05669670638256e-138, exp(normal_lcdf(-25, 0, 1))); - EXPECT_FLOAT_EQ(7.38570686148941e-133, exp(normal_lcdf(-24.5, 0, 1))); - EXPECT_FLOAT_EQ(1.3903921185497e-127, exp(normal_lcdf(-24, 0, 1))); - EXPECT_FLOAT_EQ(2.03936756324998e-122, exp(normal_lcdf(-23.5, 0, 1))); - EXPECT_FLOAT_EQ(2.33063700622065e-117, exp(normal_lcdf(-23, 0, 1))); - EXPECT_FLOAT_EQ(2.07531079906636e-112, exp(normal_lcdf(-22.5, 0, 1))); - EXPECT_FLOAT_EQ(1.43989243514508e-107, exp(normal_lcdf(-22, 0, 1))); - EXPECT_FLOAT_EQ(7.78439707718263e-103, exp(normal_lcdf(-21.5, 0, 1))); - EXPECT_FLOAT_EQ(3.27927801897904e-98, exp(normal_lcdf(-21, 0, 1))); - EXPECT_FLOAT_EQ(1.0764673258791e-93, exp(normal_lcdf(-20.5, 0, 1))); - EXPECT_FLOAT_EQ(2.75362411860623e-89, exp(normal_lcdf(-20, 0, 1))); - EXPECT_FLOAT_EQ(5.48911547566041e-85, exp(normal_lcdf(-19.5, 0, 1))); - EXPECT_FLOAT_EQ(8.52722395263098e-81, exp(normal_lcdf(-19, 0, 1))); - EXPECT_FLOAT_EQ(1.03236986895633e-76, exp(normal_lcdf(-18.5, 0, 1))); - EXPECT_FLOAT_EQ(9.74094891893715e-73, exp(normal_lcdf(-18, 0, 1))); - EXPECT_FLOAT_EQ(7.16345876623504e-69, exp(normal_lcdf(-17.5, 0, 1))); - EXPECT_FLOAT_EQ(4.10599620209891e-65, exp(normal_lcdf(-17, 0, 1))); - EXPECT_FLOAT_EQ(1.83446300316473e-61, exp(normal_lcdf(-16.5, 0, 1))); - EXPECT_FLOAT_EQ(6.38875440053809e-58, exp(normal_lcdf(-16, 0, 1))); - EXPECT_FLOAT_EQ(1.73446079179387e-54, exp(normal_lcdf(-15.5, 0, 1))); - EXPECT_FLOAT_EQ(3.67096619931275e-51, exp(normal_lcdf(-15, 0, 1))); - EXPECT_FLOAT_EQ(6.05749476441522e-48, exp(normal_lcdf(-14.5, 0, 1))); - EXPECT_FLOAT_EQ(7.7935368191928e-45, exp(normal_lcdf(-14, 0, 1))); - EXPECT_FLOAT_EQ(7.81880730565789e-42, exp(normal_lcdf(-13.5, 0, 1))); - EXPECT_FLOAT_EQ(6.11716439954988e-39, exp(normal_lcdf(-13, 0, 1))); - EXPECT_FLOAT_EQ(3.73256429887771e-36, exp(normal_lcdf(-12.5, 0, 1))); - EXPECT_FLOAT_EQ(1.77648211207768e-33, exp(normal_lcdf(-12, 0, 1))); - EXPECT_FLOAT_EQ(6.59577144611367e-31, exp(normal_lcdf(-11.5, 0, 1))); - EXPECT_FLOAT_EQ(1.91065957449868e-28, exp(normal_lcdf(-11, 0, 1))); - EXPECT_FLOAT_EQ(4.31900631780923e-26, exp(normal_lcdf(-10.5, 0, 1))); - EXPECT_FLOAT_EQ(7.61985302416053e-24, exp(normal_lcdf(-10, 0, 1))); - EXPECT_FLOAT_EQ(1.04945150753626e-21, exp(normal_lcdf(-9.5, 0, 1))); - EXPECT_FLOAT_EQ(1.12858840595384e-19, exp(normal_lcdf(-9, 0, 1))); - EXPECT_FLOAT_EQ(9.47953482220332e-18, exp(normal_lcdf(-8.5, 0, 1))); - EXPECT_FLOAT_EQ(6.22096057427178e-16, exp(normal_lcdf(-8, 0, 1))); - EXPECT_FLOAT_EQ(3.1908916729109e-14, exp(normal_lcdf(-7.5, 0, 1))); - EXPECT_FLOAT_EQ(1.27981254388584e-12, exp(normal_lcdf(-7, 0, 1))); - EXPECT_FLOAT_EQ(4.01600058385912e-11, exp(normal_lcdf(-6.5, 0, 1))); - EXPECT_FLOAT_EQ(9.86587645037698e-10, exp(normal_lcdf(-6, 0, 1))); - EXPECT_FLOAT_EQ(1.89895624658877e-08, exp(normal_lcdf(-5.5, 0, 1))); - EXPECT_FLOAT_EQ(2.86651571879194e-07, exp(normal_lcdf(-5, 0, 1))); - EXPECT_FLOAT_EQ(3.39767312473006e-06, exp(normal_lcdf(-4.5, 0, 1))); - EXPECT_FLOAT_EQ(3.16712418331199e-05, exp(normal_lcdf(-4, 0, 1))); - EXPECT_FLOAT_EQ(0.000232629079035525, exp(normal_lcdf(-3.5, 0, 1))); - EXPECT_FLOAT_EQ(0.00134989803163009, exp(normal_lcdf(-3, 0, 1))); - EXPECT_FLOAT_EQ(0.00620966532577613, exp(normal_lcdf(-2.5, 0, 1))); - EXPECT_FLOAT_EQ(0.0227501319481792, exp(normal_lcdf(-2, 0, 1))); - EXPECT_FLOAT_EQ(0.0668072012688581, exp(normal_lcdf(-1.5, 0, 1))); - EXPECT_FLOAT_EQ(0.158655253931457, exp(normal_lcdf(-1, 0, 1))); - EXPECT_FLOAT_EQ(0.308537538725987, exp(normal_lcdf(-0.5, 0, 1))); - EXPECT_FLOAT_EQ(0.5, exp(normal_lcdf(0, 0, 1))); - EXPECT_FLOAT_EQ(0.691462461274013, exp(normal_lcdf(0.5, 0, 1))); - EXPECT_FLOAT_EQ(0.841344746068543, exp(normal_lcdf(1, 0, 1))); - EXPECT_FLOAT_EQ(0.933192798731142, exp(normal_lcdf(1.5, 0, 1))); - EXPECT_FLOAT_EQ(0.977249868051821, exp(normal_lcdf(2, 0, 1))); - EXPECT_FLOAT_EQ(0.993790334674224, exp(normal_lcdf(2.5, 0, 1))); - EXPECT_FLOAT_EQ(0.99865010196837, exp(normal_lcdf(3, 0, 1))); - EXPECT_FLOAT_EQ(0.999767370920964, exp(normal_lcdf(3.5, 0, 1))); - EXPECT_FLOAT_EQ(0.999968328758167, exp(normal_lcdf(4, 0, 1))); - EXPECT_FLOAT_EQ(0.999996602326875, exp(normal_lcdf(4.5, 0, 1))); - EXPECT_FLOAT_EQ(0.999999713348428, exp(normal_lcdf(5, 0, 1))); - EXPECT_FLOAT_EQ(0.999999981010438, exp(normal_lcdf(5.5, 0, 1))); - EXPECT_FLOAT_EQ(0.999999999013412, exp(normal_lcdf(6, 0, 1))); - EXPECT_FLOAT_EQ(0.99999999995984, exp(normal_lcdf(6.5, 0, 1))); - EXPECT_FLOAT_EQ(0.99999999999872, exp(normal_lcdf(7, 0, 1))); - EXPECT_FLOAT_EQ(0.999999999999968, exp(normal_lcdf(7.5, 0, 1))); - EXPECT_FLOAT_EQ(0.999999999999999, exp(normal_lcdf(8, 0, 1))); - EXPECT_FLOAT_EQ(1, exp(normal_lcdf(8.5, 0, 1))); - EXPECT_FLOAT_EQ(1, exp(normal_lcdf(9, 0, 1))); - EXPECT_FLOAT_EQ(1, exp(normal_lcdf(9.5, 0, 1))); - EXPECT_FLOAT_EQ(1, exp(normal_lcdf(10, 0, 1))); + // The test values come from R 4.6.1 and cover the expected useful range of + // the function. When z >= 38.5, even the log of the CDF is + // indistinguishable from 0.0 in double precision. + // + // q <- + // c( + // -10^seq(8, 2, by = -1), + // seq(-50, -11, by = 1.0), + // seq(-10, 10, by = 0.5), + // seq(11, 38, by = 1.0) + // ) + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, normal_lcdf(%.17g, 0, 1));\n", + // pnorm(q[i], lower.tail = TRUE, log.p = TRUE), + // q[i] + // ) + // ) + // } + + EXPECT_FLOAT_EQ(-5000000000000019.0, normal_lcdf(-100000000, 0, 1)); + EXPECT_FLOAT_EQ(-50000000000017.039, normal_lcdf(-10000000, 0, 1)); + EXPECT_FLOAT_EQ(-500000000014.73444, normal_lcdf(-1000000, 0, 1)); + EXPECT_FLOAT_EQ(-5000000012.4318638, normal_lcdf(-100000, 0, 1)); + EXPECT_FLOAT_EQ(-50000010.129278913, normal_lcdf(-10000, 0, 1)); + EXPECT_FLOAT_EQ(-500007.82669481216, normal_lcdf(-1000, 0, 1)); + EXPECT_FLOAT_EQ(-5005.5242086942053, normal_lcdf(-100, 0, 1)); + EXPECT_FLOAT_EQ(-1254.8313611394199, normal_lcdf(-50, 0, 1)); + EXPECT_FLOAT_EQ(-1205.3111748916654, normal_lcdf(-49, 0, 1)); + EXPECT_FLOAT_EQ(-1156.7905731019453, normal_lcdf(-48, 0, 1)); + EXPECT_FLOAT_EQ(-1109.2695383172531, normal_lcdf(-47, 0, 1)); + EXPECT_FLOAT_EQ(-1062.7480519624305, normal_lcdf(-46, 0, 1)); + EXPECT_FLOAT_EQ(-1017.2260942419524, normal_lcdf(-45, 0, 1)); + EXPECT_FLOAT_EQ(-972.70364403073665, normal_lcdf(-44, 0, 1)); + EXPECT_FLOAT_EQ(-929.18067875247391, normal_lcdf(-43, 0, 1)); + EXPECT_FLOAT_EQ(-886.65717424372951, normal_lcdf(-42, 0, 1)); + EXPECT_FLOAT_EQ(-845.13310460177456, normal_lcdf(-41, 0, 1)); + EXPECT_FLOAT_EQ(-804.60844201375380, normal_lcdf(-40, 0, 1)); + EXPECT_FLOAT_EQ(-765.08315656437753, normal_lcdf(-39, 0, 1)); + EXPECT_FLOAT_EQ(-726.55721601882010, normal_lcdf(-38, 0, 1)); + EXPECT_FLOAT_EQ(-689.03058557689064, normal_lcdf(-37, 0, 1)); + EXPECT_FLOAT_EQ(-652.50322759379844, normal_lcdf(-36, 0, 1)); + EXPECT_FLOAT_EQ(-616.97510126192253, normal_lcdf(-35, 0, 1)); + EXPECT_FLOAT_EQ(-582.44616224687172, normal_lcdf(-34, 0, 1)); + EXPECT_FLOAT_EQ(-548.91636226973810, normal_lcdf(-33, 0, 1)); + EXPECT_FLOAT_EQ(-516.38564862572537, normal_lcdf(-32, 0, 1)); + EXPECT_FLOAT_EQ(-484.85396362717927, normal_lcdf(-31, 0, 1)); + EXPECT_FLOAT_EQ(-454.32124395634321, normal_lcdf(-30, 0, 1)); + EXPECT_FLOAT_EQ(-424.78741990973015, normal_lcdf(-29, 0, 1)); + EXPECT_FLOAT_EQ(-396.25241451163106, normal_lcdf(-28, 0, 1)); + EXPECT_FLOAT_EQ(-368.71614246865636, normal_lcdf(-27, 0, 1)); + EXPECT_FLOAT_EQ(-342.17850892992783, normal_lcdf(-26, 0, 1)); + EXPECT_FLOAT_EQ(-316.63940800802027, normal_lcdf(-25, 0, 1)); + EXPECT_FLOAT_EQ(-292.09872100320780, normal_lcdf(-24, 0, 1)); + EXPECT_FLOAT_EQ(-268.55631425686312, normal_lcdf(-23, 0, 1)); + EXPECT_FLOAT_EQ(-246.01203653738091, normal_lcdf(-22, 0, 1)); + EXPECT_FLOAT_EQ(-224.46571583141449, normal_lcdf(-21, 0, 1)); + EXPECT_FLOAT_EQ(-203.91715537109727, normal_lcdf(-20, 0, 1)); + EXPECT_FLOAT_EQ(-184.36612866916096, normal_lcdf(-19, 0, 1)); + EXPECT_FLOAT_EQ(-165.81237325071419, normal_lcdf(-18, 0, 1)); + EXPECT_FLOAT_EQ(-148.25558265098039, normal_lcdf(-17, 0, 1)); + EXPECT_FLOAT_EQ(-131.69539607375970, normal_lcdf(-16, 0, 1)); + EXPECT_FLOAT_EQ(-116.13138484571169, normal_lcdf(-15, 0, 1)); + EXPECT_FLOAT_EQ(-101.56303440744996, normal_lcdf(-14, 0, 1)); + EXPECT_FLOAT_EQ(-87.989719971022524, normal_lcdf(-13, 0, 1)); + EXPECT_FLOAT_EQ(-75.410673001568796, normal_lcdf(-12, 0, 1)); + EXPECT_FLOAT_EQ(-63.824934094423718, normal_lcdf(-11, 0, 1)); + EXPECT_FLOAT_EQ(-53.231285150512470, normal_lcdf(-10, 0, 1)); + EXPECT_FLOAT_EQ(-48.306019298965232, normal_lcdf(-9.5, 0, 1)); + EXPECT_FLOAT_EQ(-43.628149113332114, normal_lcdf(-9, 0, 1)); + EXPECT_FLOAT_EQ(-39.197396428217672, normal_lcdf(-8.5, 0, 1)); + EXPECT_FLOAT_EQ(-35.013437159914552, normal_lcdf(-8, 0, 1)); + EXPECT_FLOAT_EQ(-31.075890902890002, normal_lcdf(-7.5, 0, 1)); + EXPECT_FLOAT_EQ(-27.384307498811076, normal_lcdf(-7, 0, 1)); + EXPECT_FLOAT_EQ(-23.938149495161838, normal_lcdf(-6.5, 0, 1)); + EXPECT_FLOAT_EQ(-20.736768949974707, normal_lcdf(-6, 0, 1)); + EXPECT_FLOAT_EQ(-17.779376352625260, normal_lcdf(-5.5, 0, 1)); + EXPECT_FLOAT_EQ(-15.064998393988725, normal_lcdf(-5, 0, 1)); + EXPECT_FLOAT_EQ(-12.592419735713079, normal_lcdf(-4.5, 0, 1)); + EXPECT_FLOAT_EQ(-10.360101486527292, normal_lcdf(-4, 0, 1)); + EXPECT_FLOAT_EQ(-8.3660653083440941, normal_lcdf(-3.5, 0, 1)); + EXPECT_FLOAT_EQ(-6.6077262215103492, normal_lcdf(-3, 0, 1)); + EXPECT_FLOAT_EQ(-5.0816482772786902, normal_lcdf(-2.5, 0, 1)); + EXPECT_FLOAT_EQ(-3.7831843336820317, normal_lcdf(-2, 0, 1)); + EXPECT_FLOAT_EQ(-2.7059444008238898, normal_lcdf(-1.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.8410216450092636, normal_lcdf(-1, 0, 1)); + EXPECT_FLOAT_EQ(-1.1759117615936185, normal_lcdf(-0.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.69314718055994529, normal_lcdf(0, 0, 1)); + EXPECT_FLOAT_EQ(-0.36894641528865652, normal_lcdf(0.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.17275377902344988, normal_lcdf(1, 0, 1)); + EXPECT_FLOAT_EQ(-0.069143455612233992, normal_lcdf(1.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.023012909328963493, normal_lcdf(2, 0, 1)); + EXPECT_FLOAT_EQ(-0.0062290254858600024, normal_lcdf(2.5, 0, 1)); + EXPECT_FLOAT_EQ(-0.0013508099647481938, normal_lcdf(3, 0, 1)); + EXPECT_FLOAT_EQ(-0.00023265614137680445, normal_lcdf(3.5, 0, 1)); + EXPECT_FLOAT_EQ(-3.1671743377489267e-05, normal_lcdf(4, 0, 1)); + EXPECT_FLOAT_EQ(-3.3976788968344657e-06, normal_lcdf(4.5, 0, 1)); + EXPECT_FLOAT_EQ(-2.8665161296376358e-07, normal_lcdf(5, 0, 1)); + EXPECT_FLOAT_EQ(-1.8989562646189464e-08, normal_lcdf(5.5, 0, 1)); + EXPECT_FLOAT_EQ(-9.8658764552437559e-10, normal_lcdf(6, 0, 1)); + EXPECT_FLOAT_EQ(-4.0160005839397589e-11, normal_lcdf(6.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.2798125438866541e-12, normal_lcdf(7, 0, 1)); + EXPECT_FLOAT_EQ(-3.1908916729109475e-14, normal_lcdf(7.5, 0, 1)); + EXPECT_FLOAT_EQ(-6.2209605742717868e-16, normal_lcdf(8, 0, 1)); + EXPECT_FLOAT_EQ(-9.4795348222033192e-18, normal_lcdf(8.5, 0, 1)); + EXPECT_FLOAT_EQ(-1.1285884059538408e-19, normal_lcdf(9, 0, 1)); + EXPECT_FLOAT_EQ(-1.0494515075362608e-21, normal_lcdf(9.5, 0, 1)); + EXPECT_FLOAT_EQ(-7.6198530241605269e-24, normal_lcdf(10, 0, 1)); + EXPECT_FLOAT_EQ(-1.9106595744986757e-28, normal_lcdf(11, 0, 1)); + EXPECT_FLOAT_EQ(-1.7764821120776790e-33, normal_lcdf(12, 0, 1)); + EXPECT_FLOAT_EQ(-6.1171643995498803e-39, normal_lcdf(13, 0, 1)); + EXPECT_FLOAT_EQ(-7.7935368191928000e-45, normal_lcdf(14, 0, 1)); + EXPECT_FLOAT_EQ(-3.6709661993127514e-51, normal_lcdf(15, 0, 1)); + EXPECT_FLOAT_EQ(-6.3887544005380882e-58, normal_lcdf(16, 0, 1)); + EXPECT_FLOAT_EQ(-4.1059962020989074e-65, normal_lcdf(17, 0, 1)); + EXPECT_FLOAT_EQ(-9.7409489189371508e-73, normal_lcdf(18, 0, 1)); + EXPECT_FLOAT_EQ(-8.5272239526309772e-81, normal_lcdf(19, 0, 1)); + EXPECT_FLOAT_EQ(-2.7536241186062337e-89, normal_lcdf(20, 0, 1)); + EXPECT_FLOAT_EQ(-3.2792780189790367e-98, normal_lcdf(21, 0, 1)); + EXPECT_FLOAT_EQ(-1.4398924351450790e-107, normal_lcdf(22, 0, 1)); + EXPECT_FLOAT_EQ(-2.3306370062206492e-117, normal_lcdf(23, 0, 1)); + EXPECT_FLOAT_EQ(-1.3903921185497032e-127, normal_lcdf(24, 0, 1)); + EXPECT_FLOAT_EQ(-3.0566967063825616e-138, normal_lcdf(25, 0, 1)); + EXPECT_FLOAT_EQ(-2.4760633155033892e-149, normal_lcdf(26, 0, 1)); + EXPECT_FLOAT_EQ(-7.3894810068850200e-161, normal_lcdf(27, 0, 1)); + EXPECT_FLOAT_EQ(-8.1238694696594273e-173, normal_lcdf(28, 0, 1)); + EXPECT_FLOAT_EQ(-3.2897852667043802e-185, normal_lcdf(29, 0, 1)); + EXPECT_FLOAT_EQ(-4.9067139271481872e-198, normal_lcdf(30, 0, 1)); + EXPECT_FLOAT_EQ(-2.6952500812005002e-211, normal_lcdf(31, 0, 1)); + EXPECT_FLOAT_EQ(-5.4520806035123956e-225, normal_lcdf(32, 0, 1)); + EXPECT_FLOAT_EQ(-4.0611856209158557e-239, normal_lcdf(33, 0, 1)); + EXPECT_FLOAT_EQ(-1.1138987855743795e-253, normal_lcdf(34, 0, 1)); + EXPECT_FLOAT_EQ(-1.1249107064724062e-268, normal_lcdf(35, 0, 1)); + EXPECT_FLOAT_EQ(-4.1826240657972830e-284, normal_lcdf(36, 0, 1)); + EXPECT_FLOAT_EQ(-5.7255712225245771e-300, normal_lcdf(37, 0, 1)); + EXPECT_FLOAT_EQ(-2.8854283510039645e-316, normal_lcdf(38, 0, 1)); +} + +/** + * Value of log Phi against high-precision references. + * + * Each row carries the reference to 40 significant digits in a comment and + * the correctly-rounded double beside it, so a reviewer can check the + * rounding here and the reference itself in any arbitrary-precision tool. + * They were produced at 60 significant digits as `log(erfc(-s)/2)` with + * `s = (y - mu) / (sigma * sqrt(2))`, then rounded to double. + * + * The rows span all three branches of the value: `log1p(-erfc(s)/2)` for + * `s > 0`, `log(erfc(-s)) + LOG_HALF` for `-4 < s <= 0`, and the Cody + * rational approximation below that. + * + * Every row holds to 1e-14 except the last. There `s = 11.31` puts erfc into + * its far tail, where this platform's libm loses relative accuracy: + * `std::erfc(11.313708498984761)` returns 1.27775088010759500e-57 against a + * true 1.27774061192683864e-57, a relative error of 8.0e-06 that passes + * straight through. That is a libm limit rather than this function's, so the + * tolerance on that row is loose enough to stay portable. + */ +TEST(ProbNormal, lcdf_matches_high_precision_reference) { + struct value_ref { + double y; + double mu; + double sigma; + double expected; + double tol; + }; + const value_ref cases[] + = {// s=2.12132 (log1p); -0.001350809964748193798841110490517380092729 + {3.0, 0.0, 1.0, -0.0013508099647481938, 1e-14}, + // s=0.707107 (log1p); -0.1727537790234498895264831735208007300094 + {1.0, 0.0, 1.0, -0.17275377902344988, 1e-14}, + // s=0.353553 (log1p); -0.3689464152886563930656156390431773405547 + {0.5, 0.0, 1.0, -0.3689464152886564, 1e-14}, + // s=0 (erfc); -0.6931471805599453094172321214581765680755 + {0.0, 0.0, 1.0, -0.6931471805599453, 1e-14}, + // s=-0.707107 (erfc); -1.841021645009263505770783073232529021548 + {-1.0, 0.0, 1.0, -1.8410216450092636, 1e-14}, + // s=-2.12132 (erfc); -6.607726221510349543276077083251514438781 + {-3.0, 0.0, 1.0, -6.607726221510349, 1e-14}, + // s=-3.53553 (erfc); -15.06499839398872573608370479189672560507 + {-5.0, 0.0, 1.0, -15.064998393988725, 1e-14}, + // s=-4.24264 (Cody); -20.7367689499747056549688537180999188204 + {-6.0, 0.0, 1.0, -20.736768949974707, 1e-14}, + // s=-7.07107 (Cody); -53.23128515051247057834702735413120987892 + {-10.0, 0.0, 1.0, -53.23128515051247, 1e-14}, + // s=-14.1421 (Cody); -203.9171553710972639368044586545269000525 + {-20.0, 0.0, 1.0, -203.91715537109727, 1e-14}, + // s=-26.5165 (Cody); -707.6689893175071910661131734572604413308 + {-37.5, 0.0, 1.0, -707.6689893175072, 1e-14}, + // s=-42.4264 (Cody); -1805.013560680567138700666808059022928381 + {-60.0, 0.0, 1.0, -1805.0135606805673, 1e-14}, + // s=-70.7107 (Cody); -5005.52420869420508862630245733002553134 + {-100.0, 0.0, 1.0, -5005.524208694205, 1e-14}, + // s=-212.132 (Cody); -45006.62273211866335985382181649804744182 + {-300.0, 0.0, 1.0, -45006.62273211866, 1e-14}, + // s=-3.29983 (erfc); -13.38983327471644676126972305600605540117 + {-12.0, 2.0, 3.0, -13.389833274716446, 1e-14}, + // s=11.3137 (log1p); -6.388703059634194605096774759212020048496e-58 + {7.0, -1.0, 0.5, -6.388703059634195e-58, 1e-4}}; + + for (const value_ref& c : cases) { + const double got = stan::math::normal_lcdf(c.y, c.mu, c.sigma); + EXPECT_LT(std::fabs(got / c.expected - 1.0), c.tol) + << "log Phi at y = " << c.y << " got " << got << " want " << c.expected; + } +} + +TEST(ProbNormal, lcdf_with_mu_and_sigma) { + using stan::math::normal_lcdf; + + // The test values come from R 4.6.1 and check behaviour around the branch + // points in `normal_lcdf` when mu != 0 and sigma != 1. + // + // parameter_sets <- list( + // c(mu = 2, sigma = 3), + // c(mu = -50, sigma = 20), + // c(mu = 0.5, sigma = 0.01) + // ) + // branch_points <- c( + // 2.9, + // 2.5, + // 2.1, + // 1.5, + // 0.8, + // 0.1, + // 0.0, + // -2.1, + // -3.9, + // -4.0, + // -7.0, + // -17.0, + // -29.0 + // ) + // test_vals <- sort(c(branch_points - 0.01, branch_points + 0.01)) + // for (pars in parameter_sets) { + // mu <- pars[["mu"]] + // sigma <- pars[["sigma"]] + // q <- mu + sigma * sqrt(2) * test_vals + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, normal_lcdf(%.17g, %.17g, %.17g));\n", + // pnorm(q[i], mean = mu, sd = sigma, lower.tail = TRUE, log.p = + // TRUE), q[i], mu, sigma + // ) + // ) + // } + // } + + EXPECT_FLOAT_EQ(-846.21384596225801, normal_lcdf(-121.07900633333048, 2, 3)); + EXPECT_FLOAT_EQ(-845.05315712467279, normal_lcdf(-120.99415351958808, 2, 3)); + EXPECT_FLOAT_EQ(-293.44113419738113, normal_lcdf(-70.167318087899062, 2, 3)); + EXPECT_FLOAT_EQ(-292.75996176272099, normal_lcdf(-70.082465274156661, 2, 3)); + EXPECT_FLOAT_EQ(-52.362878114253100, normal_lcdf(-27.740911216706191, 2, 3)); + EXPECT_FLOAT_EQ(-52.080076508792267, normal_lcdf(-27.656058402963808, 2, 3)); + EXPECT_FLOAT_EQ(-18.763386609836139, normal_lcdf(-15.012989155348336, 2, 3)); + EXPECT_FLOAT_EQ(-18.598659064395832, normal_lcdf(-14.928136341605949, 2, 3)); + EXPECT_FLOAT_EQ(-17.947533669545574, normal_lcdf(-14.588725086636405, 2, 3)); + EXPECT_FLOAT_EQ(-17.786697578589479, normal_lcdf(-14.503872272894021, 2, 3)); + EXPECT_FLOAT_EQ(-6.5552849228594905, normal_lcdf(-6.9519718498216925, 2, 3)); + EXPECT_FLOAT_EQ(-6.4632170725824736, normal_lcdf(-6.8671190360793091, 2, 3)); + EXPECT_FLOAT_EQ(-0.70449473678943686, normal_lcdf(1.9575735931288072, 2, 3)); + EXPECT_FLOAT_EQ(-0.68192694790236175, normal_lcdf(2.042426406871193, 2, 3)); + EXPECT_FLOAT_EQ(-0.59667350159732879, normal_lcdf(2.381837661840736, 2, 3)); + EXPECT_FLOAT_EQ(-0.57658898698359096, normal_lcdf(2.4666904755831216, 2, 3)); + EXPECT_FLOAT_EQ(-0.14150397288270808, normal_lcdf(5.3516861428242359, 2, 3)); + EXPECT_FLOAT_EQ(-0.13467302655617183, normal_lcdf(5.436538956566622, 2, 3)); + EXPECT_FLOAT_EQ(-0.017706913768612061, normal_lcdf(8.3215346238077359, 2, 3)); + EXPECT_FLOAT_EQ(-0.016496955508933941, normal_lcdf(8.4063874375501211, 2, 3)); + EXPECT_FLOAT_EQ(-0.0015609874450218889, + normal_lcdf(10.867119036079309, 2, 3)); + EXPECT_FLOAT_EQ(-0.0014235903084217206, + normal_lcdf(10.951971849821692, 2, 3)); + EXPECT_FLOAT_EQ(-0.00021466697317301965, + normal_lcdf(12.564175310927022, 2, 3)); + EXPECT_FLOAT_EQ(-0.00019287133934318566, + normal_lcdf(12.649028124669407, 2, 3)); + EXPECT_FLOAT_EQ(-2.1842328461060815e-05, + normal_lcdf(14.261231585774736, 2, 3)); + EXPECT_FLOAT_EQ(-1.9328842921129592e-05, + normal_lcdf(14.346084399517119, 2, 3)); + EXPECT_FLOAT_EQ(-846.21384596225778, + normal_lcdf(-870.52670888886985, -50, 20)); + EXPECT_FLOAT_EQ(-845.05315712467279, + normal_lcdf(-869.96102346392047, -50, 20)); + EXPECT_FLOAT_EQ(-293.44113419738102, + normal_lcdf(-531.11545391932702, -50, 20)); + EXPECT_FLOAT_EQ(-292.75996176272093, + normal_lcdf(-530.54976849437764, -50, 20)); + EXPECT_FLOAT_EQ(-52.362878114253085, + normal_lcdf(-248.27274144470792, -50, 20)); + EXPECT_FLOAT_EQ(-52.080076508792246, + normal_lcdf(-247.70705601975871, -50, 20)); + EXPECT_FLOAT_EQ(-18.763386609836139, + normal_lcdf(-163.41992770232224, -50, 20)); + EXPECT_FLOAT_EQ(-18.598659064395825, + normal_lcdf(-162.85424227737298, -50, 20)); + EXPECT_FLOAT_EQ(-17.947533669545571, + normal_lcdf(-160.59150057757603, -50, 20)); + EXPECT_FLOAT_EQ(-17.786697578589479, + normal_lcdf(-160.02581515262682, -50, 20)); + EXPECT_FLOAT_EQ(-6.5552849228594887, + normal_lcdf(-109.67981233214461, -50, 20)); + EXPECT_FLOAT_EQ(-6.4632170725824709, + normal_lcdf(-109.11412690719538, -50, 20)); + EXPECT_FLOAT_EQ(-0.70449473678943675, + normal_lcdf(-50.282842712474618, -50, 20)); + EXPECT_FLOAT_EQ(-0.68192694790236175, + normal_lcdf(-49.717157287525382, -50, 20)); + EXPECT_FLOAT_EQ(-0.59667350159732879, + normal_lcdf(-47.454415587728427, -50, 20)); + EXPECT_FLOAT_EQ(-0.57658898698359096, + normal_lcdf(-46.88873016277919, -50, 20)); + EXPECT_FLOAT_EQ(-0.14150397288270808, + normal_lcdf(-27.655425714505096, -50, 20)); + EXPECT_FLOAT_EQ(-0.13467302655617192, + normal_lcdf(-27.089740289555859, -50, 20)); + EXPECT_FLOAT_EQ(-0.017706913768612061, + normal_lcdf(-7.8564358412817654, -50, 20)); + EXPECT_FLOAT_EQ(-0.016496955508933965, + normal_lcdf(-7.2907504163325285, -50, 20)); + EXPECT_FLOAT_EQ(-0.0015609874450218915, + normal_lcdf(9.1141269071953843, -50, 20)); + EXPECT_FLOAT_EQ(-0.0014235903084217224, + normal_lcdf(9.679812332144607, -50, 20)); + EXPECT_FLOAT_EQ(-0.00021466697317302003, + normal_lcdf(20.427835406180137, -50, 20)); + EXPECT_FLOAT_EQ(-0.00019287133934318566, + normal_lcdf(20.993520831129374, -50, 20)); + EXPECT_FLOAT_EQ(-2.1842328461060815e-05, + normal_lcdf(31.741543905164903, -50, 20)); + EXPECT_FLOAT_EQ(-1.9328842921129592e-05, + normal_lcdf(32.307229330114126, -50, 20)); + EXPECT_FLOAT_EQ(-846.21384596225801, + normal_lcdf(0.089736645555565042, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-845.05315712467279, + normal_lcdf(0.090019488268039738, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-293.44113419738096, + normal_lcdf(0.25944227304033651, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-292.75996176272099, + normal_lcdf(0.25972511575281115, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-52.362878114253100, + normal_lcdf(0.40086362927764602, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-52.080076508792246, + normal_lcdf(0.40114647199012066, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-18.763386609836136, + normal_lcdf(0.44329003614883888, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-18.598659064395846, + normal_lcdf(0.44357287886131347, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-17.947533669545582, + normal_lcdf(0.44470424971121197, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-17.786697578589461, + normal_lcdf(0.44498709242368661, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-6.5552849228594852, + normal_lcdf(0.47016009383392771, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-6.4632170725824754, + normal_lcdf(0.47044293654640229, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.70449473678943542, + normal_lcdf(0.49985857864376271, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.68192694790236308, + normal_lcdf(0.50014142135623729, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.59667350159733179, + normal_lcdf(0.50127279220613574, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.57658898698358851, + normal_lcdf(0.50155563491861044, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.14150397288270869, + normal_lcdf(0.51117228714274743, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.13467302655617061, + normal_lcdf(0.51145512985522212, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.017706913768612078, + normal_lcdf(0.52107178207935911, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.016496955508934111, + normal_lcdf(0.5213546247918337, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.0015609874450219117, + normal_lcdf(0.52955706345359765, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.0014235903084217016, + normal_lcdf(0.52983990616607235, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.00021466697317301539, + normal_lcdf(0.53521391770309013, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-0.00019287133934318404, + normal_lcdf(0.53549676041556471, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-2.1842328461060476e-05, + normal_lcdf(0.54087077195258249, 0.5, 0.01)); + EXPECT_FLOAT_EQ(-1.9328842921129514e-05, + normal_lcdf(0.54115361466505707, 0.5, 0.01)); } 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 33aafcd71bd..a2e39eb25b3 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 @@ -13,34 +13,141 @@ TEST(ProbStdNormal, ccdf_log_matches_lccdf) { TEST(ProbStdNormal, lccdf_tail) { using stan::math::std_normal_lccdf; - EXPECT_FLOAT_EQ(-6.661338147750941214694e-16, std_normal_lccdf(-8.0)); - EXPECT_FLOAT_EQ(-3.186340080674249758114e-14, std_normal_lccdf(-7.5)); - EXPECT_FLOAT_EQ(-1.279865102788699562477e-12, std_normal_lccdf(-7.0)); - EXPECT_FLOAT_EQ(-4.015998644826973564545e-11, std_normal_lccdf(-6.5)); - EXPECT_FLOAT_EQ(-9.865877009111571184118e-10, std_normal_lccdf(-6.0)); - EXPECT_FLOAT_EQ(-1.898956265833514866414e-08, std_normal_lccdf(-5.5)); - EXPECT_FLOAT_EQ(-2.866516130081049047962e-07, std_normal_lccdf(-5.0)); - EXPECT_FLOAT_EQ(-3.397678896843115195074e-06, std_normal_lccdf(-4.5)); - EXPECT_FLOAT_EQ(-3.167174337748932124543e-05, std_normal_lccdf(-4.0)); - EXPECT_FLOAT_EQ(-0.0002326561413768195969113, std_normal_lccdf(-3.5)); - EXPECT_FLOAT_EQ(-0.001350809964748202673598, std_normal_lccdf(-3.0)); - EXPECT_FLOAT_EQ(-0.0062290254858600267035, std_normal_lccdf(-2.5)); - EXPECT_FLOAT_EQ(-0.02301290932896348992442, std_normal_lccdf(-2.0)); - EXPECT_FLOAT_EQ(-0.06914345561223400604689, std_normal_lccdf(-1.5)); - EXPECT_FLOAT_EQ(-0.1727537790234499048836, std_normal_lccdf(-1.0)); - EXPECT_FLOAT_EQ(-0.3689464152886565151412, std_normal_lccdf(-0.5)); - EXPECT_FLOAT_EQ(-0.6931471805599452862268, std_normal_lccdf(0)); - EXPECT_FLOAT_EQ(-1.175911761593618320987, std_normal_lccdf(0.5)); - EXPECT_FLOAT_EQ(-1.841021645009263352222, std_normal_lccdf(1.0)); - EXPECT_FLOAT_EQ(-2.705944400823889317564, std_normal_lccdf(1.5)); - EXPECT_FLOAT_EQ(-3.78318433368203210776, std_normal_lccdf(2.0)); - EXPECT_FLOAT_EQ(-5.081648277278686620662, std_normal_lccdf(2.5)); - EXPECT_FLOAT_EQ(-6.607726221510342945464, std_normal_lccdf(3.0)); - EXPECT_FLOAT_EQ(-8.366065308344028395027, std_normal_lccdf(3.5)); - EXPECT_FLOAT_EQ(-10.36010148652728979357, std_normal_lccdf(4.0)); - EXPECT_FLOAT_EQ(-12.59241973571053385683, std_normal_lccdf(4.5)); - EXPECT_FLOAT_EQ(-15.06499839383403838156, std_normal_lccdf(5.0)); - EXPECT_FLOAT_EQ(-17.77937635198566113104, std_normal_lccdf(5.5)); - EXPECT_FLOAT_EQ(-20.73676889383495947072, std_normal_lccdf(6.0)); - EXPECT_FLOAT_EQ(-23.93814997800869548428, std_normal_lccdf(6.5)); + // The test values come from R 4.6.1 and cover the expected useful range of + // the function. When z <= -38.5, even the log of the CCDF is + // indistinguishable from 0.0 in double precision. + // + // q <- + // c( + // seq(-38, -11, by = 1.0), + // seq(-10, 10, by = 0.5), + // seq(11, 50, by = 1.0), + // 10^seq(2, 8, by = 1) + // ) + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, std_normal_lccdf(%.17g));\n", + // pnorm(q[i], lower.tail = FALSE, log.p = TRUE), + // q[i] + // ) + // ) + // } + + EXPECT_FLOAT_EQ(-2.8854283510039645e-316, std_normal_lccdf(-38)); + EXPECT_FLOAT_EQ(-5.7255712225245771e-300, std_normal_lccdf(-37)); + EXPECT_FLOAT_EQ(-4.1826240657972830e-284, std_normal_lccdf(-36)); + EXPECT_FLOAT_EQ(-1.1249107064724062e-268, std_normal_lccdf(-35)); + EXPECT_FLOAT_EQ(-1.1138987855743795e-253, std_normal_lccdf(-34)); + EXPECT_FLOAT_EQ(-4.0611856209158557e-239, std_normal_lccdf(-33)); + EXPECT_FLOAT_EQ(-5.4520806035123956e-225, std_normal_lccdf(-32)); + EXPECT_FLOAT_EQ(-2.6952500812005002e-211, std_normal_lccdf(-31)); + EXPECT_FLOAT_EQ(-4.9067139271481872e-198, std_normal_lccdf(-30)); + EXPECT_FLOAT_EQ(-3.2897852667043802e-185, std_normal_lccdf(-29)); + EXPECT_FLOAT_EQ(-8.1238694696594273e-173, std_normal_lccdf(-28)); + EXPECT_FLOAT_EQ(-7.3894810068850200e-161, std_normal_lccdf(-27)); + EXPECT_FLOAT_EQ(-2.4760633155033892e-149, std_normal_lccdf(-26)); + EXPECT_FLOAT_EQ(-3.0566967063825616e-138, std_normal_lccdf(-25)); + EXPECT_FLOAT_EQ(-1.3903921185497032e-127, std_normal_lccdf(-24)); + EXPECT_FLOAT_EQ(-2.3306370062206492e-117, std_normal_lccdf(-23)); + EXPECT_FLOAT_EQ(-1.4398924351450790e-107, std_normal_lccdf(-22)); + EXPECT_FLOAT_EQ(-3.2792780189790367e-98, std_normal_lccdf(-21)); + EXPECT_FLOAT_EQ(-2.7536241186062337e-89, std_normal_lccdf(-20)); + EXPECT_FLOAT_EQ(-8.5272239526309772e-81, std_normal_lccdf(-19)); + EXPECT_FLOAT_EQ(-9.7409489189371508e-73, std_normal_lccdf(-18)); + EXPECT_FLOAT_EQ(-4.1059962020989074e-65, std_normal_lccdf(-17)); + EXPECT_FLOAT_EQ(-6.3887544005380882e-58, std_normal_lccdf(-16)); + EXPECT_FLOAT_EQ(-3.6709661993127514e-51, std_normal_lccdf(-15)); + EXPECT_FLOAT_EQ(-7.7935368191928000e-45, std_normal_lccdf(-14)); + EXPECT_FLOAT_EQ(-6.1171643995498803e-39, std_normal_lccdf(-13)); + EXPECT_FLOAT_EQ(-1.7764821120776790e-33, std_normal_lccdf(-12)); + EXPECT_FLOAT_EQ(-1.9106595744986757e-28, std_normal_lccdf(-11)); + EXPECT_FLOAT_EQ(-7.6198530241605269e-24, std_normal_lccdf(-10)); + EXPECT_FLOAT_EQ(-1.0494515075362608e-21, std_normal_lccdf(-9.5)); + EXPECT_FLOAT_EQ(-1.1285884059538408e-19, std_normal_lccdf(-9)); + EXPECT_FLOAT_EQ(-9.4795348222033192e-18, std_normal_lccdf(-8.5)); + EXPECT_FLOAT_EQ(-6.2209605742717868e-16, std_normal_lccdf(-8)); + EXPECT_FLOAT_EQ(-3.1908916729109475e-14, std_normal_lccdf(-7.5)); + EXPECT_FLOAT_EQ(-1.2798125438866541e-12, std_normal_lccdf(-7)); + EXPECT_FLOAT_EQ(-4.0160005839397589e-11, std_normal_lccdf(-6.5)); + EXPECT_FLOAT_EQ(-9.8658764552437559e-10, std_normal_lccdf(-6)); + EXPECT_FLOAT_EQ(-1.8989562646189464e-08, std_normal_lccdf(-5.5)); + EXPECT_FLOAT_EQ(-2.8665161296376358e-07, std_normal_lccdf(-5)); + EXPECT_FLOAT_EQ(-3.3976788968344657e-06, std_normal_lccdf(-4.5)); + EXPECT_FLOAT_EQ(-3.1671743377489267e-05, std_normal_lccdf(-4)); + EXPECT_FLOAT_EQ(-0.00023265614137680445, std_normal_lccdf(-3.5)); + EXPECT_FLOAT_EQ(-0.0013508099647481938, std_normal_lccdf(-3)); + EXPECT_FLOAT_EQ(-0.0062290254858600024, std_normal_lccdf(-2.5)); + EXPECT_FLOAT_EQ(-0.023012909328963493, std_normal_lccdf(-2)); + EXPECT_FLOAT_EQ(-0.069143455612233992, std_normal_lccdf(-1.5)); + EXPECT_FLOAT_EQ(-0.17275377902344988, std_normal_lccdf(-1)); + EXPECT_FLOAT_EQ(-0.36894641528865652, std_normal_lccdf(-0.5)); + EXPECT_FLOAT_EQ(-0.69314718055994529, std_normal_lccdf(0)); + EXPECT_FLOAT_EQ(-1.1759117615936185, std_normal_lccdf(0.5)); + EXPECT_FLOAT_EQ(-1.8410216450092636, std_normal_lccdf(1)); + EXPECT_FLOAT_EQ(-2.7059444008238898, std_normal_lccdf(1.5)); + EXPECT_FLOAT_EQ(-3.7831843336820317, std_normal_lccdf(2)); + EXPECT_FLOAT_EQ(-5.0816482772786902, std_normal_lccdf(2.5)); + EXPECT_FLOAT_EQ(-6.6077262215103492, std_normal_lccdf(3)); + EXPECT_FLOAT_EQ(-8.3660653083440941, std_normal_lccdf(3.5)); + EXPECT_FLOAT_EQ(-10.360101486527292, std_normal_lccdf(4)); + EXPECT_FLOAT_EQ(-12.592419735713079, std_normal_lccdf(4.5)); + EXPECT_FLOAT_EQ(-15.064998393988725, std_normal_lccdf(5)); + EXPECT_FLOAT_EQ(-17.779376352625260, std_normal_lccdf(5.5)); + EXPECT_FLOAT_EQ(-20.736768949974707, std_normal_lccdf(6)); + EXPECT_FLOAT_EQ(-23.938149495161838, std_normal_lccdf(6.5)); + EXPECT_FLOAT_EQ(-27.384307498811076, std_normal_lccdf(7)); + EXPECT_FLOAT_EQ(-31.075890902890002, std_normal_lccdf(7.5)); + EXPECT_FLOAT_EQ(-35.013437159914552, std_normal_lccdf(8)); + EXPECT_FLOAT_EQ(-39.197396428217672, std_normal_lccdf(8.5)); + EXPECT_FLOAT_EQ(-43.628149113332114, std_normal_lccdf(9)); + EXPECT_FLOAT_EQ(-48.306019298965232, std_normal_lccdf(9.5)); + EXPECT_FLOAT_EQ(-53.231285150512470, std_normal_lccdf(10)); + EXPECT_FLOAT_EQ(-63.824934094423718, std_normal_lccdf(11)); + EXPECT_FLOAT_EQ(-75.410673001568796, std_normal_lccdf(12)); + EXPECT_FLOAT_EQ(-87.989719971022524, std_normal_lccdf(13)); + EXPECT_FLOAT_EQ(-101.56303440744996, std_normal_lccdf(14)); + EXPECT_FLOAT_EQ(-116.13138484571169, std_normal_lccdf(15)); + EXPECT_FLOAT_EQ(-131.69539607375970, std_normal_lccdf(16)); + EXPECT_FLOAT_EQ(-148.25558265098039, std_normal_lccdf(17)); + EXPECT_FLOAT_EQ(-165.81237325071419, std_normal_lccdf(18)); + EXPECT_FLOAT_EQ(-184.36612866916096, std_normal_lccdf(19)); + EXPECT_FLOAT_EQ(-203.91715537109727, std_normal_lccdf(20)); + EXPECT_FLOAT_EQ(-224.46571583141449, std_normal_lccdf(21)); + EXPECT_FLOAT_EQ(-246.01203653738091, std_normal_lccdf(22)); + EXPECT_FLOAT_EQ(-268.55631425686312, std_normal_lccdf(23)); + EXPECT_FLOAT_EQ(-292.09872100320780, std_normal_lccdf(24)); + EXPECT_FLOAT_EQ(-316.63940800802027, std_normal_lccdf(25)); + EXPECT_FLOAT_EQ(-342.17850892992783, std_normal_lccdf(26)); + EXPECT_FLOAT_EQ(-368.71614246865636, std_normal_lccdf(27)); + EXPECT_FLOAT_EQ(-396.25241451163106, std_normal_lccdf(28)); + EXPECT_FLOAT_EQ(-424.78741990973015, std_normal_lccdf(29)); + EXPECT_FLOAT_EQ(-454.32124395634321, std_normal_lccdf(30)); + EXPECT_FLOAT_EQ(-484.85396362717927, std_normal_lccdf(31)); + EXPECT_FLOAT_EQ(-516.38564862572537, std_normal_lccdf(32)); + EXPECT_FLOAT_EQ(-548.91636226973810, std_normal_lccdf(33)); + EXPECT_FLOAT_EQ(-582.44616224687172, std_normal_lccdf(34)); + EXPECT_FLOAT_EQ(-616.97510126192253, std_normal_lccdf(35)); + EXPECT_FLOAT_EQ(-652.50322759379844, std_normal_lccdf(36)); + EXPECT_FLOAT_EQ(-689.03058557689064, std_normal_lccdf(37)); + EXPECT_FLOAT_EQ(-726.55721601882010, std_normal_lccdf(38)); + EXPECT_FLOAT_EQ(-765.08315656437753, std_normal_lccdf(39)); + EXPECT_FLOAT_EQ(-804.60844201375380, std_normal_lccdf(40)); + EXPECT_FLOAT_EQ(-845.13310460177456, std_normal_lccdf(41)); + EXPECT_FLOAT_EQ(-886.65717424372951, std_normal_lccdf(42)); + EXPECT_FLOAT_EQ(-929.18067875247391, std_normal_lccdf(43)); + EXPECT_FLOAT_EQ(-972.70364403073665, std_normal_lccdf(44)); + EXPECT_FLOAT_EQ(-1017.2260942419524, std_normal_lccdf(45)); + EXPECT_FLOAT_EQ(-1062.7480519624305, std_normal_lccdf(46)); + EXPECT_FLOAT_EQ(-1109.2695383172531, std_normal_lccdf(47)); + EXPECT_FLOAT_EQ(-1156.7905731019453, std_normal_lccdf(48)); + EXPECT_FLOAT_EQ(-1205.3111748916654, std_normal_lccdf(49)); + EXPECT_FLOAT_EQ(-1254.8313611394199, std_normal_lccdf(50)); + EXPECT_FLOAT_EQ(-5005.5242086942053, std_normal_lccdf(100)); + EXPECT_FLOAT_EQ(-500007.82669481216, std_normal_lccdf(1000)); + EXPECT_FLOAT_EQ(-50000010.129278913, std_normal_lccdf(10000)); + EXPECT_FLOAT_EQ(-5000000012.4318638, std_normal_lccdf(100000)); + EXPECT_FLOAT_EQ(-500000000014.73444, std_normal_lccdf(1000000)); + EXPECT_FLOAT_EQ(-50000000000017.039, std_normal_lccdf(10000000)); + EXPECT_FLOAT_EQ(-5000000000000019.0, std_normal_lccdf(100000000)); } diff --git a/test/unit/math/prim/prob/std_normal_cdf_log_test.cpp b/test/unit/math/prim/prob/std_normal_cdf_log_test.cpp index 6cb75bb93f8..8181bd29c7f 100644 --- a/test/unit/math/prim/prob/std_normal_cdf_log_test.cpp +++ b/test/unit/math/prim/prob/std_normal_cdf_log_test.cpp @@ -1,113 +1,219 @@ #include #include +#include TEST(ProbStdNormal, cdf_log_matches_lcdf) { double y = 0.8; EXPECT_FLOAT_EQ((stan::math::std_normal_lcdf(y)), (stan::math::std_normal_cdf_log(y))); - EXPECT_FLOAT_EQ((stan::math::std_normal_lcdf(y)), + EXPECT_FLOAT_EQ((stan::math::std_normal_lcdf(y)), (stan::math::std_normal_cdf_log(y))); } TEST(ProbStdNormal, lcdf_tails) { using stan::math::std_normal_lcdf; - using std::exp; - EXPECT_FLOAT_EQ(4.60535300958196e-308, exp(std_normal_lcdf(-37.5))); - EXPECT_FLOAT_EQ(5.72557122252458e-300, exp(std_normal_lcdf(-37))); - EXPECT_FLOAT_EQ(5.54472571307484e-292, exp(std_normal_lcdf(-36.5))); - EXPECT_FLOAT_EQ(4.18262406579728e-284, exp(std_normal_lcdf(-36))); - EXPECT_FLOAT_EQ(2.45769154066194e-276, exp(std_normal_lcdf(-35.5))); - EXPECT_FLOAT_EQ(1.12491070647241e-268, exp(std_normal_lcdf(-35))); - EXPECT_FLOAT_EQ(4.01072896657726e-261, exp(std_normal_lcdf(-34.5))); - EXPECT_FLOAT_EQ(1.11389878557438e-253, exp(std_normal_lcdf(-34))); - EXPECT_FLOAT_EQ(2.40983869512039e-246, exp(std_normal_lcdf(-33.5))); - EXPECT_FLOAT_EQ(4.06118562091586e-239, exp(std_normal_lcdf(-33))); - EXPECT_FLOAT_EQ(5.33142435967881e-232, exp(std_normal_lcdf(-32.5))); - EXPECT_FLOAT_EQ(5.4520806035124e-225, exp(std_normal_lcdf(-32))); - EXPECT_FLOAT_EQ(4.34323260103177e-218, exp(std_normal_lcdf(-31.5))); - EXPECT_FLOAT_EQ(2.6952500812005e-211, exp(std_normal_lcdf(-31))); - EXPECT_FLOAT_EQ(1.30293791317808e-204, exp(std_normal_lcdf(-30.5))); - EXPECT_FLOAT_EQ(4.90671392714819e-198, exp(std_normal_lcdf(-30))); - EXPECT_FLOAT_EQ(1.43947455222918e-191, exp(std_normal_lcdf(-29.5))); - EXPECT_FLOAT_EQ(3.28978526670438e-185, exp(std_normal_lcdf(-29))); - EXPECT_FLOAT_EQ(5.85714125380634e-179, exp(std_normal_lcdf(-28.5))); - EXPECT_FLOAT_EQ(8.12386946965943e-173, exp(std_normal_lcdf(-28))); - EXPECT_FLOAT_EQ(8.77817055687808e-167, exp(std_normal_lcdf(-27.5))); - EXPECT_FLOAT_EQ(7.38948100688502e-161, exp(std_normal_lcdf(-27))); - EXPECT_FLOAT_EQ(4.84616266030332e-155, exp(std_normal_lcdf(-26.5))); - EXPECT_FLOAT_EQ(2.47606331550339e-149, exp(std_normal_lcdf(-26))); - EXPECT_FLOAT_EQ(9.85623651896393e-144, exp(std_normal_lcdf(-25.5))); - EXPECT_FLOAT_EQ(3.05669670638256e-138, exp(std_normal_lcdf(-25))); - EXPECT_FLOAT_EQ(7.38570686148941e-133, exp(std_normal_lcdf(-24.5))); - EXPECT_FLOAT_EQ(1.3903921185497e-127, exp(std_normal_lcdf(-24))); - EXPECT_FLOAT_EQ(2.03936756324998e-122, exp(std_normal_lcdf(-23.5))); - EXPECT_FLOAT_EQ(2.33063700622065e-117, exp(std_normal_lcdf(-23))); - EXPECT_FLOAT_EQ(2.07531079906636e-112, exp(std_normal_lcdf(-22.5))); - EXPECT_FLOAT_EQ(1.43989243514508e-107, exp(std_normal_lcdf(-22))); - EXPECT_FLOAT_EQ(7.78439707718263e-103, exp(std_normal_lcdf(-21.5))); - EXPECT_FLOAT_EQ(3.27927801897904e-98, exp(std_normal_lcdf(-21))); - EXPECT_FLOAT_EQ(1.0764673258791e-93, exp(std_normal_lcdf(-20.5))); - EXPECT_FLOAT_EQ(2.75362411860623e-89, exp(std_normal_lcdf(-20))); - EXPECT_FLOAT_EQ(5.48911547566041e-85, exp(std_normal_lcdf(-19.5))); - EXPECT_FLOAT_EQ(8.52722395263098e-81, exp(std_normal_lcdf(-19))); - EXPECT_FLOAT_EQ(1.03236986895633e-76, exp(std_normal_lcdf(-18.5))); - EXPECT_FLOAT_EQ(9.74094891893715e-73, exp(std_normal_lcdf(-18))); - EXPECT_FLOAT_EQ(7.16345876623504e-69, exp(std_normal_lcdf(-17.5))); - EXPECT_FLOAT_EQ(4.10599620209891e-65, exp(std_normal_lcdf(-17))); - EXPECT_FLOAT_EQ(1.83446300316473e-61, exp(std_normal_lcdf(-16.5))); - EXPECT_FLOAT_EQ(6.38875440053809e-58, exp(std_normal_lcdf(-16))); - EXPECT_FLOAT_EQ(1.73446079179387e-54, exp(std_normal_lcdf(-15.5))); - EXPECT_FLOAT_EQ(3.67096619931275e-51, exp(std_normal_lcdf(-15))); - EXPECT_FLOAT_EQ(6.05749476441522e-48, exp(std_normal_lcdf(-14.5))); - EXPECT_FLOAT_EQ(7.7935368191928e-45, exp(std_normal_lcdf(-14))); - EXPECT_FLOAT_EQ(7.81880730565789e-42, exp(std_normal_lcdf(-13.5))); - EXPECT_FLOAT_EQ(6.11716439954988e-39, exp(std_normal_lcdf(-13))); - EXPECT_FLOAT_EQ(3.73256429887771e-36, exp(std_normal_lcdf(-12.5))); - EXPECT_FLOAT_EQ(1.77648211207768e-33, exp(std_normal_lcdf(-12))); - EXPECT_FLOAT_EQ(6.59577144611367e-31, exp(std_normal_lcdf(-11.5))); - EXPECT_FLOAT_EQ(1.91065957449868e-28, exp(std_normal_lcdf(-11))); - EXPECT_FLOAT_EQ(4.31900631780923e-26, exp(std_normal_lcdf(-10.5))); - EXPECT_FLOAT_EQ(7.61985302416053e-24, exp(std_normal_lcdf(-10))); - EXPECT_FLOAT_EQ(1.04945150753626e-21, exp(std_normal_lcdf(-9.5))); - EXPECT_FLOAT_EQ(1.12858840595384e-19, exp(std_normal_lcdf(-9))); - EXPECT_FLOAT_EQ(9.47953482220332e-18, exp(std_normal_lcdf(-8.5))); - EXPECT_FLOAT_EQ(6.22096057427178e-16, exp(std_normal_lcdf(-8))); - EXPECT_FLOAT_EQ(3.1908916729109e-14, exp(std_normal_lcdf(-7.5))); - EXPECT_FLOAT_EQ(1.27981254388584e-12, exp(std_normal_lcdf(-7))); - EXPECT_FLOAT_EQ(4.01600058385912e-11, exp(std_normal_lcdf(-6.5))); - EXPECT_FLOAT_EQ(9.86587645037698e-10, exp(std_normal_lcdf(-6))); - EXPECT_FLOAT_EQ(1.89895624658877e-08, exp(std_normal_lcdf(-5.5))); - EXPECT_FLOAT_EQ(2.86651571879194e-07, exp(std_normal_lcdf(-5))); - EXPECT_FLOAT_EQ(3.39767312473006e-06, exp(std_normal_lcdf(-4.5))); - EXPECT_FLOAT_EQ(3.16712418331199e-05, exp(std_normal_lcdf(-4))); - EXPECT_FLOAT_EQ(0.000232629079035525, exp(std_normal_lcdf(-3.5))); - EXPECT_FLOAT_EQ(0.00134989803163009, exp(std_normal_lcdf(-3))); - EXPECT_FLOAT_EQ(0.00620966532577613, exp(std_normal_lcdf(-2.5))); - EXPECT_FLOAT_EQ(0.0227501319481792, exp(std_normal_lcdf(-2))); - EXPECT_FLOAT_EQ(0.0668072012688581, exp(std_normal_lcdf(-1.5))); - EXPECT_FLOAT_EQ(0.158655253931457, exp(std_normal_lcdf(-1))); - EXPECT_FLOAT_EQ(0.308537538725987, exp(std_normal_lcdf(-0.5))); - EXPECT_FLOAT_EQ(0.5, exp(std_normal_lcdf(0))); - EXPECT_FLOAT_EQ(0.691462461274013, exp(std_normal_lcdf(0.5))); - EXPECT_FLOAT_EQ(0.841344746068543, exp(std_normal_lcdf(1))); - EXPECT_FLOAT_EQ(0.933192798731142, exp(std_normal_lcdf(1.5))); - EXPECT_FLOAT_EQ(0.977249868051821, exp(std_normal_lcdf(2))); - EXPECT_FLOAT_EQ(0.993790334674224, exp(std_normal_lcdf(2.5))); - EXPECT_FLOAT_EQ(0.99865010196837, exp(std_normal_lcdf(3))); - EXPECT_FLOAT_EQ(0.999767370920964, exp(std_normal_lcdf(3.5))); - EXPECT_FLOAT_EQ(0.999968328758167, exp(std_normal_lcdf(4))); - EXPECT_FLOAT_EQ(0.999996602326875, exp(std_normal_lcdf(4.5))); - EXPECT_FLOAT_EQ(0.999999713348428, exp(std_normal_lcdf(5))); - EXPECT_FLOAT_EQ(0.999999981010438, exp(std_normal_lcdf(5.5))); - EXPECT_FLOAT_EQ(0.999999999013412, exp(std_normal_lcdf(6))); - EXPECT_FLOAT_EQ(0.99999999995984, exp(std_normal_lcdf(6.5))); - EXPECT_FLOAT_EQ(0.99999999999872, exp(std_normal_lcdf(7))); - EXPECT_FLOAT_EQ(0.999999999999968, exp(std_normal_lcdf(7.5))); - EXPECT_FLOAT_EQ(0.999999999999999, exp(std_normal_lcdf(8))); - EXPECT_FLOAT_EQ(1, exp(std_normal_lcdf(8.5))); - EXPECT_FLOAT_EQ(1, exp(std_normal_lcdf(9))); - EXPECT_FLOAT_EQ(1, exp(std_normal_lcdf(9.5))); - EXPECT_FLOAT_EQ(1, exp(std_normal_lcdf(10))); + // The test values come from R 4.6.1 and cover the expected useful range of + // the function. When z >= 38.5, even the log of the CDF is + // indistinguishable from 0.0 in double precision. + // + // q <- + // c( + // -10^seq(8, 2, by = -1), + // seq(-50, -11, by = 1.0), + // seq(-10, 10, by = 0.5), + // seq(11, 38, by = 1.0) + // ) + // for (i in 1:length(q)) { + // cat( + // sprintf( + // "EXPECT_FLOAT_EQ(%#.17g, std_normal_lcdf(%.17g));\n", + // pnorm(q[i], lower.tail = TRUE, log.p = TRUE), + // q[i] + // ) + // ) + // } + + EXPECT_FLOAT_EQ(-5000000000000019.0, std_normal_lcdf(-100000000)); + EXPECT_FLOAT_EQ(-50000000000017.039, std_normal_lcdf(-10000000)); + EXPECT_FLOAT_EQ(-500000000014.73444, std_normal_lcdf(-1000000)); + EXPECT_FLOAT_EQ(-5000000012.4318638, std_normal_lcdf(-100000)); + EXPECT_FLOAT_EQ(-50000010.129278913, std_normal_lcdf(-10000)); + EXPECT_FLOAT_EQ(-500007.82669481216, std_normal_lcdf(-1000)); + EXPECT_FLOAT_EQ(-5005.5242086942053, std_normal_lcdf(-100)); + EXPECT_FLOAT_EQ(-1254.8313611394199, std_normal_lcdf(-50)); + EXPECT_FLOAT_EQ(-1205.3111748916654, std_normal_lcdf(-49)); + EXPECT_FLOAT_EQ(-1156.7905731019453, std_normal_lcdf(-48)); + EXPECT_FLOAT_EQ(-1109.2695383172531, std_normal_lcdf(-47)); + EXPECT_FLOAT_EQ(-1062.7480519624305, std_normal_lcdf(-46)); + EXPECT_FLOAT_EQ(-1017.2260942419524, std_normal_lcdf(-45)); + EXPECT_FLOAT_EQ(-972.70364403073665, std_normal_lcdf(-44)); + EXPECT_FLOAT_EQ(-929.18067875247391, std_normal_lcdf(-43)); + EXPECT_FLOAT_EQ(-886.65717424372951, std_normal_lcdf(-42)); + EXPECT_FLOAT_EQ(-845.13310460177456, std_normal_lcdf(-41)); + EXPECT_FLOAT_EQ(-804.60844201375380, std_normal_lcdf(-40)); + EXPECT_FLOAT_EQ(-765.08315656437753, std_normal_lcdf(-39)); + EXPECT_FLOAT_EQ(-726.55721601882010, std_normal_lcdf(-38)); + EXPECT_FLOAT_EQ(-689.03058557689064, std_normal_lcdf(-37)); + EXPECT_FLOAT_EQ(-652.50322759379844, std_normal_lcdf(-36)); + EXPECT_FLOAT_EQ(-616.97510126192253, std_normal_lcdf(-35)); + EXPECT_FLOAT_EQ(-582.44616224687172, std_normal_lcdf(-34)); + EXPECT_FLOAT_EQ(-548.91636226973810, std_normal_lcdf(-33)); + EXPECT_FLOAT_EQ(-516.38564862572537, std_normal_lcdf(-32)); + EXPECT_FLOAT_EQ(-484.85396362717927, std_normal_lcdf(-31)); + EXPECT_FLOAT_EQ(-454.32124395634321, std_normal_lcdf(-30)); + EXPECT_FLOAT_EQ(-424.78741990973015, std_normal_lcdf(-29)); + EXPECT_FLOAT_EQ(-396.25241451163106, std_normal_lcdf(-28)); + EXPECT_FLOAT_EQ(-368.71614246865636, std_normal_lcdf(-27)); + EXPECT_FLOAT_EQ(-342.17850892992783, std_normal_lcdf(-26)); + EXPECT_FLOAT_EQ(-316.63940800802027, std_normal_lcdf(-25)); + EXPECT_FLOAT_EQ(-292.09872100320780, std_normal_lcdf(-24)); + EXPECT_FLOAT_EQ(-268.55631425686312, std_normal_lcdf(-23)); + EXPECT_FLOAT_EQ(-246.01203653738091, std_normal_lcdf(-22)); + EXPECT_FLOAT_EQ(-224.46571583141449, std_normal_lcdf(-21)); + EXPECT_FLOAT_EQ(-203.91715537109727, std_normal_lcdf(-20)); + EXPECT_FLOAT_EQ(-184.36612866916096, std_normal_lcdf(-19)); + EXPECT_FLOAT_EQ(-165.81237325071419, std_normal_lcdf(-18)); + EXPECT_FLOAT_EQ(-148.25558265098039, std_normal_lcdf(-17)); + EXPECT_FLOAT_EQ(-131.69539607375970, std_normal_lcdf(-16)); + EXPECT_FLOAT_EQ(-116.13138484571169, std_normal_lcdf(-15)); + EXPECT_FLOAT_EQ(-101.56303440744996, std_normal_lcdf(-14)); + EXPECT_FLOAT_EQ(-87.989719971022524, std_normal_lcdf(-13)); + EXPECT_FLOAT_EQ(-75.410673001568796, std_normal_lcdf(-12)); + EXPECT_FLOAT_EQ(-63.824934094423718, std_normal_lcdf(-11)); + EXPECT_FLOAT_EQ(-53.231285150512470, std_normal_lcdf(-10)); + EXPECT_FLOAT_EQ(-48.306019298965232, std_normal_lcdf(-9.5)); + EXPECT_FLOAT_EQ(-43.628149113332114, std_normal_lcdf(-9)); + EXPECT_FLOAT_EQ(-39.197396428217672, std_normal_lcdf(-8.5)); + EXPECT_FLOAT_EQ(-35.013437159914552, std_normal_lcdf(-8)); + EXPECT_FLOAT_EQ(-31.075890902890002, std_normal_lcdf(-7.5)); + EXPECT_FLOAT_EQ(-27.384307498811076, std_normal_lcdf(-7)); + EXPECT_FLOAT_EQ(-23.938149495161838, std_normal_lcdf(-6.5)); + EXPECT_FLOAT_EQ(-20.736768949974707, std_normal_lcdf(-6)); + EXPECT_FLOAT_EQ(-17.779376352625260, std_normal_lcdf(-5.5)); + EXPECT_FLOAT_EQ(-15.064998393988725, std_normal_lcdf(-5)); + EXPECT_FLOAT_EQ(-12.592419735713079, std_normal_lcdf(-4.5)); + EXPECT_FLOAT_EQ(-10.360101486527292, std_normal_lcdf(-4)); + EXPECT_FLOAT_EQ(-8.3660653083440941, std_normal_lcdf(-3.5)); + EXPECT_FLOAT_EQ(-6.6077262215103492, std_normal_lcdf(-3)); + EXPECT_FLOAT_EQ(-5.0816482772786902, std_normal_lcdf(-2.5)); + EXPECT_FLOAT_EQ(-3.7831843336820317, std_normal_lcdf(-2)); + EXPECT_FLOAT_EQ(-2.7059444008238898, std_normal_lcdf(-1.5)); + EXPECT_FLOAT_EQ(-1.8410216450092636, std_normal_lcdf(-1)); + EXPECT_FLOAT_EQ(-1.1759117615936185, std_normal_lcdf(-0.5)); + EXPECT_FLOAT_EQ(-0.69314718055994529, std_normal_lcdf(0)); + EXPECT_FLOAT_EQ(-0.36894641528865652, std_normal_lcdf(0.5)); + EXPECT_FLOAT_EQ(-0.17275377902344988, std_normal_lcdf(1)); + EXPECT_FLOAT_EQ(-0.069143455612233992, std_normal_lcdf(1.5)); + EXPECT_FLOAT_EQ(-0.023012909328963493, std_normal_lcdf(2)); + EXPECT_FLOAT_EQ(-0.0062290254858600024, std_normal_lcdf(2.5)); + EXPECT_FLOAT_EQ(-0.0013508099647481938, std_normal_lcdf(3)); + EXPECT_FLOAT_EQ(-0.00023265614137680445, std_normal_lcdf(3.5)); + EXPECT_FLOAT_EQ(-3.1671743377489267e-05, std_normal_lcdf(4)); + EXPECT_FLOAT_EQ(-3.3976788968344657e-06, std_normal_lcdf(4.5)); + EXPECT_FLOAT_EQ(-2.8665161296376358e-07, std_normal_lcdf(5)); + EXPECT_FLOAT_EQ(-1.8989562646189464e-08, std_normal_lcdf(5.5)); + EXPECT_FLOAT_EQ(-9.8658764552437559e-10, std_normal_lcdf(6)); + EXPECT_FLOAT_EQ(-4.0160005839397589e-11, std_normal_lcdf(6.5)); + EXPECT_FLOAT_EQ(-1.2798125438866541e-12, std_normal_lcdf(7)); + EXPECT_FLOAT_EQ(-3.1908916729109475e-14, std_normal_lcdf(7.5)); + EXPECT_FLOAT_EQ(-6.2209605742717868e-16, std_normal_lcdf(8)); + EXPECT_FLOAT_EQ(-9.4795348222033192e-18, std_normal_lcdf(8.5)); + EXPECT_FLOAT_EQ(-1.1285884059538408e-19, std_normal_lcdf(9)); + EXPECT_FLOAT_EQ(-1.0494515075362608e-21, std_normal_lcdf(9.5)); + EXPECT_FLOAT_EQ(-7.6198530241605269e-24, std_normal_lcdf(10)); + EXPECT_FLOAT_EQ(-1.9106595744986757e-28, std_normal_lcdf(11)); + EXPECT_FLOAT_EQ(-1.7764821120776790e-33, std_normal_lcdf(12)); + EXPECT_FLOAT_EQ(-6.1171643995498803e-39, std_normal_lcdf(13)); + EXPECT_FLOAT_EQ(-7.7935368191928000e-45, std_normal_lcdf(14)); + EXPECT_FLOAT_EQ(-3.6709661993127514e-51, std_normal_lcdf(15)); + EXPECT_FLOAT_EQ(-6.3887544005380882e-58, std_normal_lcdf(16)); + EXPECT_FLOAT_EQ(-4.1059962020989074e-65, std_normal_lcdf(17)); + EXPECT_FLOAT_EQ(-9.7409489189371508e-73, std_normal_lcdf(18)); + EXPECT_FLOAT_EQ(-8.5272239526309772e-81, std_normal_lcdf(19)); + EXPECT_FLOAT_EQ(-2.7536241186062337e-89, std_normal_lcdf(20)); + EXPECT_FLOAT_EQ(-3.2792780189790367e-98, std_normal_lcdf(21)); + EXPECT_FLOAT_EQ(-1.4398924351450790e-107, std_normal_lcdf(22)); + EXPECT_FLOAT_EQ(-2.3306370062206492e-117, std_normal_lcdf(23)); + EXPECT_FLOAT_EQ(-1.3903921185497032e-127, std_normal_lcdf(24)); + EXPECT_FLOAT_EQ(-3.0566967063825616e-138, std_normal_lcdf(25)); + EXPECT_FLOAT_EQ(-2.4760633155033892e-149, std_normal_lcdf(26)); + EXPECT_FLOAT_EQ(-7.3894810068850200e-161, std_normal_lcdf(27)); + EXPECT_FLOAT_EQ(-8.1238694696594273e-173, std_normal_lcdf(28)); + EXPECT_FLOAT_EQ(-3.2897852667043802e-185, std_normal_lcdf(29)); + EXPECT_FLOAT_EQ(-4.9067139271481872e-198, std_normal_lcdf(30)); + EXPECT_FLOAT_EQ(-2.6952500812005002e-211, std_normal_lcdf(31)); + EXPECT_FLOAT_EQ(-5.4520806035123956e-225, std_normal_lcdf(32)); + EXPECT_FLOAT_EQ(-4.0611856209158557e-239, std_normal_lcdf(33)); + EXPECT_FLOAT_EQ(-1.1138987855743795e-253, std_normal_lcdf(34)); + EXPECT_FLOAT_EQ(-1.1249107064724062e-268, std_normal_lcdf(35)); + EXPECT_FLOAT_EQ(-4.1826240657972830e-284, std_normal_lcdf(36)); + EXPECT_FLOAT_EQ(-5.7255712225245771e-300, std_normal_lcdf(37)); + EXPECT_FLOAT_EQ(-2.8854283510039645e-316, std_normal_lcdf(38)); +} + +/** + * Value of log Phi against high-precision references. + * + * Each row carries the reference to 40 significant digits in a comment and + * the correctly-rounded double beside it, so a reviewer can check the + * rounding here and the reference itself in any arbitrary-precision tool. + * They were produced at 60 significant digits as `log(erfc(-s)/2)` with + * `s = y * INV_SQRT_TWO`, then rounded to double. + * + * The rows span all three branches of the value: `log1p(-erfc(s)/2)` for + * `s > 0`, `log(erfc(-s)) + LOG_HALF` for `-4 < s <= 0`, and the Cody + * rational approximation below that. + * + * Every row holds to 1e-14 except the last. There `s = 11.31` puts erfc into + * its far tail, where this platform's libm loses relative accuracy: + * `std::erfc(11.313708498984761)` returns 1.27775088010759500e-57 against a + * true 1.27774061192683864e-57, a relative error of 8.0e-06 that passes + * straight through. That is a libm limit rather than this function's, so the + * tolerance on that row is loose enough to stay portable. + */ +TEST(ProbStdNormal, lcdf_matches_high_precision_reference) { + struct value_ref { + double y; + double expected; + double tol; + }; + const value_ref cases[] + = {// s=2.12132 (log1p); -0.001350809964748193798841110490517380092729 + {3.0, -0.0013508099647481938, 1e-14}, + // s=0.707107 (log1p); -0.1727537790234498895264831735208007300094 + {1.0, -0.17275377902344988, 1e-14}, + // s=0.353553 (log1p); -0.3689464152886563930656156390431773405547 + {0.5, -0.3689464152886564, 1e-14}, + // s=0 (erfc); -0.6931471805599453094172321214581765680755 + {0.0, -0.6931471805599453, 1e-14}, + // s=-0.707107 (erfc); -1.841021645009263505770783073232529021548 + {-1.0, -1.8410216450092636, 1e-14}, + // s=-2.12132 (erfc); -6.607726221510349543276077083251514438781 + {-3.0, -6.607726221510349, 1e-14}, + // s=-3.53553 (erfc); -15.06499839398872573608370479189672560507 + {-5.0, -15.064998393988725, 1e-14}, + // s=-4.24264 (Cody); -20.7367689499747056549688537180999188204 + {-6.0, -20.736768949974707, 1e-14}, + // s=-7.07107 (Cody); -53.23128515051247057834702735413120987892 + {-10.0, -53.23128515051247, 1e-14}, + // s=-14.1421 (Cody); -203.9171553710972639368044586545269000525 + {-20.0, -203.91715537109727, 1e-14}, + // s=-26.5165 (Cody); -707.6689893175071910661131734572604413308 + {-37.5, -707.6689893175072, 1e-14}, + // s=-42.4264 (Cody); -1805.013560680567138700666808059022928381 + {-60.0, -1805.0135606805673, 1e-14}, + // s=-70.7107 (Cody); -5005.52420869420508862630245733002553134 + {-100.0, -5005.524208694205, 1e-14}, + // s=-212.132 (Cody); -45006.62273211866335985382181649804744182 + {-300.0, -45006.62273211866, 1e-14}, + // s=11.3137 (log1p); -6.388703059634194605096774759212020048496e-58 + {16.0, -6.388703059634195e-58, 1e-4}}; + + for (const value_ref& c : cases) { + const double got = stan::math::std_normal_lcdf(c.y); + EXPECT_LT(std::fabs(got / c.expected - 1.0), c.tol) + << "log Phi at y = " << c.y << " got " << got << " want " << c.expected; + } }