Merge pull request #3363 from jaburgoyne/fix/normal_lccdf

The normal_lccdf() and std_normal_lccdf() are updated to use normal_lcdf() and std_normal_lcdf() implementations the test suites for all four functions are improved.

Author: jaburgoyne

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;

stan/math/opencl/prim/normal_lccdf.hpp

@@ -2,16 +2,13 @@
 #define STAN_MATH_OPENCL_PRIM_NORMAL_LCCDF_HPP
 #ifdef STAN_OPENCL
 
-#include <stan/math/prim/meta.hpp>
-#include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/elt_divide.hpp>
-#include <stan/math/prim/fun/elt_multiply.hpp>
-#include <stan/math/opencl/kernel_generator.hpp>
-#include <stan/math/prim/functor/partials_propagator.hpp>
+#include <stan/math/opencl/prim/normal_lcdf.hpp>
 
 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<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
 inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> 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<T_y_cl, T_loc_cl, T_scale_cl>;
-  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<double> 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<double> lccdf_cl;
-  matrix_cl<double> y_deriv_cl;
-  matrix_cl<double> mu_deriv_cl;
-  matrix_cl<double> 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<is_autodiff_v<T_y_cl>>(y_deriv),
-                    calc_if<is_autodiff_v<T_loc_cl>>(mu_deriv),
-                    calc_if<is_autodiff_v<T_scale_cl>>(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<T_y_cl>) {
-    partials<0>(ops_partials) = std::move(y_deriv_cl);
-  }
-  if constexpr (is_autodiff_v<T_loc_cl>) {
-    partials<1>(ops_partials) = std::move(mu_deriv_cl);
-  }
-  if constexpr (is_autodiff_v<T_scale_cl>) {
-    partials<2>(ops_partials) = std::move(sigma_deriv_cl);
-  }
-  return ops_partials.build(lccdf);
+  return normal_lcdf<internal::normal_lccdf_opencl_func>(-y, -mu, sigma);
 }
 
 }  // namespace math

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<T_y_cl, T_loc_cl,
                                                 T_scale_cl>* = nullptr,
     require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
 inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> 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;
 

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 <stan/math/prim/meta.hpp>
-#include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/elt_divide.hpp>
-#include <stan/math/prim/fun/elt_multiply.hpp>
-#include <stan/math/opencl/kernel_generator.hpp>
-#include <stan/math/prim/functor/partials_propagator.hpp>
+#include <stan/math/opencl/prim/std_normal_lcdf.hpp>
 
 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 <typename T_y_cl,
           require_all_prim_or_rev_kernel_expression_t<T_y_cl>* = nullptr,
           require_any_not_stan_scalar_t<T_y_cl>* = nullptr>
 inline return_type_t<T_y_cl> std_normal_lccdf(const T_y_cl& y) {
-  static constexpr const char* function = "std_normal_lccdf(OpenCL)";
-  using T_partials_return = partials_return_t<T_y_cl>;
-  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<double> lccdf_cl;
-  matrix_cl<double> y_deriv_cl;
-
-  results(check_y_not_nan, lccdf_cl, y_deriv_cl) = expressions(
-      y_not_nan_expr, lccdf_expr, calc_if<is_autodiff_v<T_y_cl>>(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<T_y_cl>) {
-    partials<0>(ops_partials) = std::move(y_deriv_cl);
-  }
-  return ops_partials.build(lccdf);
+  return std_normal_lcdf<internal::std_normal_lccdf_opencl_func>(-y);
 }
 
 }  // namespace math

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 <typename T_y_cl,
+template <const char* func = internal::std_normal_lcdf_opencl_func,
+          typename T_y_cl,
           require_all_prim_or_rev_kernel_expression_t<T_y_cl>* = nullptr,
           require_any_not_stan_scalar_t<T_y_cl>* = nullptr>
 inline return_type_t<T_y_cl> 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;
 

stan/math/prim/prob/normal_cdf_log.hpp

@@ -14,7 +14,7 @@ template <typename T_y, typename T_loc, typename T_scale>
 inline return_type_t<T_y, T_loc, T_scale> normal_cdf_log(const T_y& y,
                                                          const T_loc& mu,
                                                          const T_scale& sigma) {
-  return normal_lcdf<T_y, T_loc, T_scale>(y, mu, sigma);
+  return normal_lcdf(y, mu, sigma);
 }
 
 }  // namespace math