cauchy.hpp 13 KB

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