1 // Copyright John Maddock 2006, 2007.
2 // Copyright Paul A. Bristow 2007.
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)
8 #ifndef BOOST_STATS_CAUCHY_HPP
9 #define BOOST_STATS_CAUCHY_HPP
13 #pragma warning(disable : 4127) // conditional expression is constant
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>
24 namespace boost{ namespace math
27 template <class RealType, class Policy>
28 class cauchy_distribution;
33 template <class RealType, class Policy>
34 RealType cdf_imp(const cauchy_distribution<RealType, Policy>& dist, const RealType& x, bool complement)
37 // This calculates the cdf of the Cauchy distribution and/or its complement.
39 // The usual formula for the Cauchy cdf is:
41 // cdf = 0.5 + atan(x)/pi
43 // But that suffers from cancellation error as x -> -INF.
45 // Recall that for x < 0:
47 // atan(x) = -pi/2 - atan(1/x)
49 // Substituting into the above we get:
51 // CDF = -atan(1/x) ; x < 0
53 // So the proceedure is to calculate the cdf for -fabs(x)
54 // using the above formula, and then subtract from 1 when required
57 BOOST_MATH_STD_USING // for ADL of std functions
58 static const char* function = "boost::math::cdf(cauchy<%1%>&, %1%)";
60 RealType location = dist.location();
61 RealType scale = dist.scale();
62 if(false == detail::check_location(function, location, &result, Policy()))
66 if(false == detail::check_scale(function, scale, &result, Policy()))
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);
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);
78 if(false == detail::check_x(function, x, &result, Policy()))
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.
87 result = -atan(1 / mx) / constants::pi<RealType>();
88 return (((x > location) != complement) ? 1 - result : result);
91 template <class RealType, class Policy>
92 RealType quantile_imp(
93 const cauchy_distribution<RealType, Policy>& dist,
97 // This routine implements the quantile for the Cauchy distribution,
98 // the value p may be the probability, or its complement if complement=true.
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.
105 static const char* function = "boost::math::quantile(cauchy<%1%>&, %1%)";
106 BOOST_MATH_STD_USING // for ADL of std functions
109 RealType location = dist.location();
110 RealType scale = dist.scale();
111 if(false == detail::check_location(function, location, &result, Policy()))
115 if(false == detail::check_scale(function, scale, &result, Policy()))
119 if(false == detail::check_probability(function, p, &result, Policy()))
126 return (complement ? -1 : 1) * policies::raise_overflow_error<RealType>(function, 0, Policy());
130 return (complement ? 1 : -1) * policies::raise_overflow_error<RealType>(function, 0, Policy());
133 RealType P = p - floor(p); // argument reduction of p:
138 if(P == 0.5) // special case:
142 result = -scale / tan(constants::pi<RealType>() * P);
143 return complement ? RealType(location - result) : RealType(location + result);
146 } // namespace detail
148 template <class RealType = double, class Policy = policies::policy<> >
149 class cauchy_distribution
152 typedef RealType value_type;
153 typedef Policy policy_type;
155 cauchy_distribution(RealType l_location = 0, RealType l_scale = 1)
156 : m_a(l_location), m_hg(l_scale)
158 static const char* function = "boost::math::cauchy_distribution<%1%>::cauchy_distribution";
160 detail::check_location(function, l_location, &result, Policy());
161 detail::check_scale(function, l_scale, &result, Policy());
162 } // cauchy_distribution
164 RealType location()const
168 RealType scale()const
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.
178 typedef cauchy_distribution<double> cauchy;
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)
185 return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
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.
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)
200 return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
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.
209 template <class RealType, class Policy>
210 inline RealType pdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x)
212 BOOST_MATH_STD_USING // for ADL of std functions
214 static const char* function = "boost::math::pdf(cauchy<%1%>&, %1%)";
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()))
222 if(false == detail::check_location("boost::math::pdf(cauchy<%1%>&, %1%)", location, &result, Policy()))
226 if((boost::math::isinf)(x))
228 return 0; // pdf + and - infinity is zero.
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.
236 if(false == detail::check_x(function, x, &result, Policy()))
241 RealType xs = (x - location) / scale;
242 result = 1 / (constants::pi<RealType>() * scale * (1 + xs * xs));
246 template <class RealType, class Policy>
247 inline RealType cdf(const cauchy_distribution<RealType, Policy>& dist, const RealType& x)
249 return detail::cdf_imp(dist, x, false);
252 template <class RealType, class Policy>
253 inline RealType quantile(const cauchy_distribution<RealType, Policy>& dist, const RealType& p)
255 return detail::quantile_imp(dist, p, false);
258 template <class RealType, class Policy>
259 inline RealType cdf(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c)
261 return detail::cdf_imp(c.dist, c.param, true);
264 template <class RealType, class Policy>
265 inline RealType quantile(const complemented2_type<cauchy_distribution<RealType, Policy>, RealType>& c)
267 return detail::quantile_imp(c.dist, c.param, true);
268 } // quantile complement
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);
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());
283 template <class RealType, class Policy>
284 inline RealType variance(const cauchy_distribution<RealType, Policy>& /*dist*/)
286 // There is no variance:
287 typedef typename Policy::assert_undefined_type assert_type;
288 BOOST_STATIC_ASSERT(assert_type::value == 0);
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());
297 template <class RealType, class Policy>
298 inline RealType mode(const cauchy_distribution<RealType, Policy>& dist)
300 return dist.location();
303 template <class RealType, class Policy>
304 inline RealType median(const cauchy_distribution<RealType, Policy>& dist)
306 return dist.location();
308 template <class RealType, class Policy>
309 inline RealType skewness(const cauchy_distribution<RealType, Policy>& /*dist*/)
311 // There is no skewness:
312 typedef typename Policy::assert_undefined_type assert_type;
313 BOOST_STATIC_ASSERT(assert_type::value == 0);
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?
322 template <class RealType, class Policy>
323 inline RealType kurtosis(const cauchy_distribution<RealType, Policy>& /*dist*/)
325 // There is no kurtosis:
326 typedef typename Policy::assert_undefined_type assert_type;
327 BOOST_STATIC_ASSERT(assert_type::value == 0);
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());
336 template <class RealType, class Policy>
337 inline RealType kurtosis_excess(const cauchy_distribution<RealType, Policy>& /*dist*/)
339 // There is no kurtosis excess:
340 typedef typename Policy::assert_undefined_type assert_type;
341 BOOST_STATIC_ASSERT(assert_type::value == 0);
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());
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>
362 #endif // BOOST_STATS_CAUCHY_HPP