Merge pull request #1234 from WarrenWeckesser/simplify-cauchy-cdf

Simplify the formula for the CDF of the Cauchy distribution.
boostorg · Jan 17, 2025 · 024a4fe · 024a4fe
2 parents 933faf4 + 7d6f320
commit 024a4fe
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 22 deletions.
diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp
@@ -45,23 +45,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
    //
    // This calculates the cdf of the Cauchy distribution and/or its complement.
    //
-   // The usual formula for the Cauchy cdf is:
+   // This implementation uses the formula
    //
-   // cdf = 0.5 + atan(x)/pi
+   //     cdf = atan2(1, -x)/pi
    //
-   // But that suffers from cancellation error as x -> -INF.
-   //
-   // Recall that for x < 0:
-   //
-   // atan(x) = -pi/2 - atan(1/x)
-   //
-   // Substituting into the above we get:
-   //
-   // CDF = -atan(1/x)/pi  ; x < 0
-   //
-   // So the procedure is to calculate the cdf for -fabs(x)
-   // using the above formula, and then subtract from 1 when required
-   // to get the result.
+   // where x is the standardized (i.e. shifted and scaled) domain variable.
    //
    BOOST_MATH_STD_USING // for ADL of std functions
    constexpr auto function = "boost::math::cdf(cauchy<%1%>&, %1%)";
@@ -99,13 +87,8 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
    { // Catches x == NaN
       return result;
    }
-   RealType mx = -fabs((x - location) / scale); // scale is > 0
-   if(mx > -tools::epsilon<RealType>() / 8)
-   {  // special case first: x extremely close to location.
-      return static_cast<RealType>(0.5f);
-   }
-   result = -atan(1 / mx) / constants::pi<RealType>();
-   return (((x > location) != complement) ? 1 - result : result);
+   RealType x_std = static_cast<RealType>((complement) ? 1 : -1)*(x - location) / scale;
+   return atan2(static_cast<RealType>(1), x_std) / constants::pi<RealType>();
 } // cdf
 
 template <class RealType, class Policy>

diff --git a/test/test_cauchy.cpp b/test/test_cauchy.cpp
@@ -124,6 +124,22 @@ void test_spots(RealType T)
          static_cast<RealType>(-10.0)),              // x
          static_cast<RealType>(0.031725517430553569514977118601302L),                // probability.
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         cauchy_distribution<RealType>(),
+         static_cast<RealType>(-15000000.0)),
+         static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      // Test the CDF at -max_value()/4.
+      // For an input x of this magnitude, the reference value is 4/|x|/pi.
+      ::boost::math::cdf(
+         cauchy_distribution<RealType>(),
+         -boost::math::tools::max_value<RealType>()/4),
+         static_cast<RealType>(4)
+                      / boost::math::tools::max_value<RealType>()
+                      / boost::math::constants::pi<RealType>(),
+         tolerance); // %
 
    //
    // Complements:
@@ -188,6 +204,22 @@ void test_spots(RealType T)
          static_cast<RealType>(-10.0))),              // x
          static_cast<RealType>(0.9682744825694464304850228813987L),                // probability.
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         complement(cauchy_distribution<RealType>(),
+         static_cast<RealType>(15000000.0))),
+         static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      // Test the complemented CDF at max_value()/4.
+      // For an input x of this magnitude, the reference value is 4/x/pi.
+      ::boost::math::cdf(
+         complement(cauchy_distribution<RealType>(),
+         boost::math::tools::max_value<RealType>()/4)),
+         static_cast<RealType>(4)
+                      / boost::math::tools::max_value<RealType>()
+                      / boost::math::constants::pi<RealType>(),
+         tolerance); // %
 
    //
    // Quantiles: