cauchy.hpp 12 KB

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