To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.
The primary repository for this project is hosted at https://github.com/sonic-visualiser/sv-dependency-builds .
This repository is a read-only copy which is updated automatically every hour.
root / any / include / boost / math / distributions / cauchy.hpp @ 160:cff480c41f97
History | View | Annotate | Download (12.2 KB)
| 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 |