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

source code of boost/libs/math/include/boost/math/distributions/cauchy.hpp