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.

Statistics Download as Zip
| Branch: | Tag: | Revision:

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