boost_mirror/math
1
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2025-05-11 21:33:52 +00:00
Code Issues Packages Projects Releases Wiki Activity

Simplify the formula for the CDF of the Cauchy distribution.

Browse Source 
The Cauchy CDF function may be expressed as atan2(1, -x)/pi.
This commit is contained in:
Warren Weckesser 2025-01-16 21:15:40 -05:00
parent 933faf4d2c
commit 7d6f32045a
2 changed files with 37 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
include/boost/math/distributions/cauchy.hpp
View File

@ -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>

32
test/test_cauchy.cpp
View File

@ -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: