Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simplify the formula for the CDF of the Cauchy distribution. #1234

Merged
merged 1 commit into from
Jan 17, 2025
+37 −22
Merged

Simplify the formula for the CDF of the Cauchy distribution. #1234

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 5 additions & 22 deletions include/boost/math/distributions/cauchy.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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%)";
Expand Down Expand Up @@ -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>
Expand Down
32 changes: 32 additions & 0 deletions test/test_cauchy.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down
Loading