ceph/src/boost/libs/math/include/boost/math/distributions/cauchy.hpp

   1 // Copyright John Maddock 2006, 2007.
   2 // Copyright Paul A. Bristow 2007.
   3
   4 //  Use, modification and distribution are subject to the
   5 //  Boost Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_STATS_CAUCHY_HPP
   9 #define BOOST_STATS_CAUCHY_HPP
  10
  11 #ifdef _MSC_VER
  12 #pragma warning(push)
  13 #pragma warning(disable : 4127) // conditional expression is constant
  14 #endif
  15
  16 #include <boost/math/distributions/fwd.hpp>
  17 #include <boost/math/constants/constants.hpp>
  18 #include <boost/math/distributions/complement.hpp>
  19 #include <boost/math/distributions/detail/common_error_handling.hpp>
  20 #include <boost/config/no_tr1/cmath.hpp>
  21
  22 #include <utility>
  23
  24 namespace boost{ namespace math
  25 {
  26
  27 template <class RealType, class Policy>
  28 class cauchy_distribution;
  29
  30 namespace detail
  31 {
  32
  33 template <class RealType, class Policy>
  34 RealType cdf_imp(const cauchy_distribution<RealType, Policy>& dist, const RealType& x, bool complement)
  35 {
  36    //
  37    // This calculates the cdf of the Cauchy distribution and/or its complement.
  38    //
  39    // The usual formula for the Cauchy cdf is:
  40    //
  41    // cdf = 0.5 + atan(x)/pi
  42    //
  43    // But that suffers from cancellation error as x -> -INF.
  44    //
  45    // Recall that for x < 0:
  46    //
  47    // atan(x) = -pi/2 - atan(1/x)
  48    //
  49    // Substituting into the above we get:
  50    //
  51    // CDF = -atan(1/x)  ; x < 0
  52    //
  53    // So the proceedure is to calculate the cdf for -fabs(x)
  54    // using the above formula, and then subtract from 1 when required
  55    // to get the result.
  56    //
  57    BOOST_MATH_STD_USING // for ADL of std functions
  58    static const char* function = "boost::math::cdf(cauchy<%1%>&, %1%)";
  59    RealType result = 0;
  60    RealType location = dist.location();
  61    RealType scale = dist.scale();
  62    if(false == detail::check_location(function, location, &result, Policy()))
  63    {
  64      return result;
  65    }
  66    if(false == detail::check_scale(function, scale, &result, Policy()))
  67    {
  68       return result;
  69    }
  70    if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
  71    { // cdf +infinity is unity.
  72      return static_cast<RealType>((complement) ? 0 : 1);
  73    }
  74    if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
  75    { // cdf -infinity is zero.
  76      return static_cast<RealType>((complement) ? 1 : 0);
  77    }
  78    if(false == detail::check_x(function, x, &result, Policy()))
  79    { // Catches x == NaN
  80       return result;
  81    }
  82    RealType mx = -fabs((x - location) / scale); // scale is > 0
  83    if(mx > -tools::epsilon<RealType>() / 8)
  84    {  // special case first: x extremely close to location.
  85       return 0.5;
  86    }
  87    result = -atan(1 / mx) / constants::pi<RealType>();
  88    return (((x > location) != complement) ? 1 - result : result);
  89 } // cdf
  90
  91 template <class RealType, class Policy>
  92 RealType quantile_imp(
  93       const cauchy_distribution<RealType, Policy>& dist,
  94       const RealType& p,
  95       bool complement)
  96 {
  97    // This routine implements the quantile for the Cauchy distribution,
  98    // the value p may be the probability, or its complement if complement=true.
  99    //
 100    // The procedure first performs argument reduction on p to avoid error
 101    // when calculating the tangent, then calulates the distance from the
 102    // mid-point of the distribution.  This is either added or subtracted
 103    // from the location parameter depending on whether `complement` is true.
 104    //
 105    static const char* function = "boost::math::quantile(cauchy<%1%>&, %1%)";
 106    BOOST_MATH_STD_USING // for ADL of std functions
 107
 108    RealType result = 0;
 109    RealType location = dist.location();
 110    RealType scale = dist.scale();
 111    if(false == detail::check_location(function, location, &result, Policy()))
 112    {
 113      return result;
 114    }
 115    if(false == detail::check_scale(function, scale, &result, Policy()))
 116    {
 117       return result;
 118    }
 119    if(false == detail::check_probability(function, p, &result, Policy()))
 120    {
 121       return result;
 122    }
 123    // Special cases:
 124    if(p == 1)
 125    {
 126       return (complement ? -1 : 1) * policies::raise_overflow_error<RealType>(function, 0, Policy());
 127    }
 128    if(p == 0)
 129    {
 130       return (complement ? 1 : -1) * policies::raise_overflow_error<RealType>(function, 0, Policy());
 131    }
 132
 133    RealType P = p - floor(p);   // argument reduction of p:
 134    if(P > 0.5)
 135    {
 136       P = P - 1;
 137    }
 138    if(P == 0.5)   // special case:
 139    {
 140       return location;
 141    }
 142    result = -scale / tan(constants::pi<RealType>() * P);
 143    return complement ? RealType(location - result) : RealType(location + result);
 144 } // quantile
 145
 146 } // namespace detail
 147
 148 template <class RealType = double, class Policy = policies::policy<> >
 149 class cauchy_distribution
 150 {
 151 public:
 152    typedef RealType value_type;
 153    typedef Policy policy_type;
 154
 155    cauchy_distribution(RealType l_location = 0, RealType l_scale = 1)
 156       : m_a(l_location), m_hg(l_scale)
 157    {
 158     static const char* function = "boost::math::cauchy_distribution<%1%>::cauchy_distribution";
 159      RealType result;
 160      detail::check_location(function, l_location, &result, Policy());
 161      detail::check_scale(function, l_scale, &result, Policy());
 162    } // cauchy_distribution
 163
 164    RealType location()const
 165    {
 166       return m_a;
 167    }
 168    RealType scale()const
 169    {
 170       return m_hg;
 171    }
 172
 173 private:
 174    RealType m_a;    // The location, this is the median of the distribution.
 175    RealType m_hg;   // The scale )or shape), this is the half width at half height.
 176 };
 177
 178 typedef cauchy_distribution<double> cauchy;
 179
 180 template <class RealType, class Policy>
 181 inline const std::pair<RealType, RealType> range(const cauchy_distribution<RealType, Policy>&)
 182 { // Range of permissible values for random variable x.
 183   if (std::numeric_limits<RealType>::has_infinity)
 184   {
 185      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
 186   }
 187   else
 188   { // Can only use max_value.
 189    using boost::math::tools::max_value;
 190    return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max.
 191   }
 192 }
 193
 194 template <class RealType, class Policy>
 195 inline const std::pair<RealType, RealType> support(const cauchy_distribution<RealType, Policy>& )
 196 { // Range of supported values for random variable x.
 197    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 198   if (std::numeric_limits<RealType>::has_infinity)
 199   {
 200      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
 201   }
 202   else
 203   { // Can only use max_value.
 204      using boost::math::tools::max_value;
 205      return std::pair<RealType, RealType>(-tools::max_value<RealType>(), max_value<RealType>()); // - to + max.
 206   }
 207 }
 208
 209 template <class RealType, class Policy>
 210 inline RealType pdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x)
 211 {
 212    BOOST_MATH_STD_USING  // for ADL of std functions
 213
 214    static const char* function = "boost::math::pdf(cauchy<%1%>&, %1%)";
 215    RealType result = 0;
 216    RealType location = dist.location();
 217    RealType scale = dist.scale();
 218    if(false == detail::check_scale("boost::math::pdf(cauchy<%1%>&, %1%)", scale, &result, Policy()))
 219    {
 220       return result;
 221    }
 222    if(false == detail::check_location("boost::math::pdf(cauchy<%1%>&, %1%)", location, &result, Policy()))
 223    {
 224       return result;
 225    }
 226    if((boost::math::isinf)(x))
 227    {
 228      return 0; // pdf + and - infinity is zero.
 229    }
 230    // These produce MSVC 4127 warnings, so the above used instead.
 231    //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())
 232    //{ // pdf + and - infinity is zero.
 233    //  return 0;
 234    //}
 235
 236    if(false == detail::check_x(function, x, &result, Policy()))
 237    { // Catches x = NaN
 238       return result;
 239    }
 240
 241    RealType xs = (x - location) / scale;
 242    result = 1 / (constants::pi<RealType>() * scale * (1 + xs * xs));
 243    return result;
 244 } // pdf
 245
 246 template <class RealType, class Policy>
 247 inline RealType cdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x)
 248 {
 249    return detail::cdf_imp(dist, x, false);
 250 } // cdf
 251
 252 template <class RealType, class Policy>
 253 inline RealType quantile(const cauchy_distribution<RealType, Policy>& dist, const RealType& p)
 254 {
 255    return detail::quantile_imp(dist, p, false);
 256 } // quantile
 257
 258 template <class RealType, class Policy>
 259 inline RealType cdf(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c)
 260 {
 261    return detail::cdf_imp(c.dist, c.param, true);
 262 } //  cdf complement
 263
 264 template <class RealType, class Policy>
 265 inline RealType quantile(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c)
 266 {
 267    return detail::quantile_imp(c.dist, c.param, true);
 268 } // quantile complement
 269
 270 template <class RealType, class Policy>
 271 inline RealType mean(const cauchy_distribution<RealType, Policy>&)
 272 {  // There is no mean:
 273    typedef typename Policy::assert_undefined_type assert_type;
 274    BOOST_STATIC_ASSERT(assert_type::value == 0);
 275
 276    return policies::raise_domain_error<RealType>(
 277       "boost::math::mean(cauchy<%1%>&)",
 278       "The Cauchy distribution does not have a mean: "
 279       "the only possible return value is %1%.",
 280       std::numeric_limits<RealType>::quiet_NaN(), Policy());
 281 }
 282
 283 template <class RealType, class Policy>
 284 inline RealType variance(const cauchy_distribution<RealType, Policy>& /*dist*/)
 285 {
 286    // There is no variance:
 287    typedef typename Policy::assert_undefined_type assert_type;
 288    BOOST_STATIC_ASSERT(assert_type::value == 0);
 289
 290    return policies::raise_domain_error<RealType>(
 291       "boost::math::variance(cauchy<%1%>&)",
 292       "The Cauchy distribution does not have a variance: "
 293       "the only possible return value is %1%.",
 294       std::numeric_limits<RealType>::quiet_NaN(), Policy());
 295 }
 296
 297 template <class RealType, class Policy>
 298 inline RealType mode(const cauchy_distribution<RealType, Policy>& dist)
 299 {
 300    return dist.location();
 301 }
 302
 303 template <class RealType, class Policy>
 304 inline RealType median(const cauchy_distribution<RealType, Policy>& dist)
 305 {
 306    return dist.location();
 307 }
 308 template <class RealType, class Policy>
 309 inline RealType skewness(const cauchy_distribution<RealType, Policy>& /*dist*/)
 310 {
 311    // There is no skewness:
 312    typedef typename Policy::assert_undefined_type assert_type;
 313    BOOST_STATIC_ASSERT(assert_type::value == 0);
 314
 315    return policies::raise_domain_error<RealType>(
 316       "boost::math::skewness(cauchy<%1%>&)",
 317       "The Cauchy distribution does not have a skewness: "
 318       "the only possible return value is %1%.",
 319       std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
 320 }
 321
 322 template <class RealType, class Policy>
 323 inline RealType kurtosis(const cauchy_distribution<RealType, Policy>& /*dist*/)
 324 {
 325    // There is no kurtosis:
 326    typedef typename Policy::assert_undefined_type assert_type;
 327    BOOST_STATIC_ASSERT(assert_type::value == 0);
 328
 329    return policies::raise_domain_error<RealType>(
 330       "boost::math::kurtosis(cauchy<%1%>&)",
 331       "The Cauchy distribution does not have a kurtosis: "
 332       "the only possible return value is %1%.",
 333       std::numeric_limits<RealType>::quiet_NaN(), Policy());
 334 }
 335
 336 template <class RealType, class Policy>
 337 inline RealType kurtosis_excess(const cauchy_distribution<RealType, Policy>& /*dist*/)
 338 {
 339    // There is no kurtosis excess:
 340    typedef typename Policy::assert_undefined_type assert_type;
 341    BOOST_STATIC_ASSERT(assert_type::value == 0);
 342
 343    return policies::raise_domain_error<RealType>(
 344       "boost::math::kurtosis_excess(cauchy<%1%>&)",
 345       "The Cauchy distribution does not have a kurtosis: "
 346       "the only possible return value is %1%.",
 347       std::numeric_limits<RealType>::quiet_NaN(), Policy());
 348 }
 349
 350 } // namespace math
 351 } // namespace boost
 352
 353 #ifdef _MSC_VER
 354 #pragma warning(pop)
 355 #endif
 356
 357 // This include must be at the end, *after* the accessors
 358 // for this distribution have been defined, in order to
 359 // keep compilers that support two-phase lookup happy.
 360 #include <boost/math/distributions/detail/derived_accessors.hpp>
 361
 362 #endif // BOOST_STATS_CAUCHY_HPP