rayleigh.hpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382
  1. // Copyright Paul A. Bristow 2007.
  2. // Copyright Matt Borland 2023.
  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_rayleigh_HPP
  7. #define BOOST_STATS_rayleigh_HPP
  8. #include <boost/math/tools/config.hpp>
  9. #include <boost/math/tools/cstdint.hpp>
  10. #include <boost/math/tools/tuple.hpp>
  11. #include <boost/math/tools/numeric_limits.hpp>
  12. #include <boost/math/distributions/fwd.hpp>
  13. #include <boost/math/constants/constants.hpp>
  14. #include <boost/math/special_functions/log1p.hpp>
  15. #include <boost/math/special_functions/expm1.hpp>
  16. #include <boost/math/distributions/complement.hpp>
  17. #include <boost/math/distributions/detail/common_error_handling.hpp>
  18. #ifdef _MSC_VER
  19. # pragma warning(push)
  20. # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  21. #endif
  22. namespace boost{ namespace math{
  23. namespace detail
  24. { // Error checks:
  25. template <class RealType, class Policy>
  26. BOOST_MATH_GPU_ENABLED inline bool verify_sigma(const char* function, RealType sigma, RealType* presult, const Policy& pol)
  27. {
  28. if((sigma <= 0) || (!(boost::math::isfinite)(sigma)))
  29. {
  30. *presult = policies::raise_domain_error<RealType>(
  31. function,
  32. "The scale parameter \"sigma\" must be > 0 and finite, but was: %1%.", sigma, pol);
  33. return false;
  34. }
  35. return true;
  36. } // bool verify_sigma
  37. template <class RealType, class Policy>
  38. BOOST_MATH_GPU_ENABLED inline bool verify_rayleigh_x(const char* function, RealType x, RealType* presult, const Policy& pol)
  39. {
  40. if((x < 0) || (boost::math::isnan)(x))
  41. {
  42. *presult = policies::raise_domain_error<RealType>(
  43. function,
  44. "The random variable must be >= 0, but was: %1%.", x, pol);
  45. return false;
  46. }
  47. return true;
  48. } // bool verify_rayleigh_x
  49. } // namespace detail
  50. template <class RealType = double, class Policy = policies::policy<> >
  51. class rayleigh_distribution
  52. {
  53. public:
  54. using value_type = RealType;
  55. using policy_type = Policy;
  56. BOOST_MATH_GPU_ENABLED explicit rayleigh_distribution(RealType l_sigma = 1)
  57. : m_sigma(l_sigma)
  58. {
  59. RealType err;
  60. detail::verify_sigma("boost::math::rayleigh_distribution<%1%>::rayleigh_distribution", l_sigma, &err, Policy());
  61. } // rayleigh_distribution
  62. BOOST_MATH_GPU_ENABLED RealType sigma()const
  63. { // Accessor.
  64. return m_sigma;
  65. }
  66. private:
  67. RealType m_sigma;
  68. }; // class rayleigh_distribution
  69. using rayleigh = rayleigh_distribution<double>;
  70. #ifdef __cpp_deduction_guides
  71. template <class RealType>
  72. rayleigh_distribution(RealType)->rayleigh_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  73. #endif
  74. template <class RealType, class Policy>
  75. BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const rayleigh_distribution<RealType, Policy>& /*dist*/)
  76. { // Range of permissible values for random variable x.
  77. using boost::math::tools::max_value;
  78. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), boost::math::numeric_limits<RealType>::has_infinity ? boost::math::numeric_limits<RealType>::infinity() : max_value<RealType>());
  79. }
  80. template <class RealType, class Policy>
  81. BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const rayleigh_distribution<RealType, Policy>& /*dist*/)
  82. { // Range of supported values for random variable x.
  83. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  84. using boost::math::tools::max_value;
  85. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  86. }
  87. template <class RealType, class Policy>
  88. BOOST_MATH_GPU_ENABLED inline RealType pdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
  89. {
  90. BOOST_MATH_STD_USING // for ADL of std function exp.
  91. RealType sigma = dist.sigma();
  92. RealType result = 0;
  93. constexpr auto function = "boost::math::pdf(const rayleigh_distribution<%1%>&, %1%)";
  94. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  95. {
  96. return result;
  97. }
  98. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  99. {
  100. return result;
  101. }
  102. if((boost::math::isinf)(x))
  103. {
  104. return 0;
  105. }
  106. RealType sigmasqr = sigma * sigma;
  107. result = x * (exp(-(x * x) / ( 2 * sigmasqr))) / sigmasqr;
  108. return result;
  109. } // pdf
  110. template <class RealType, class Policy>
  111. BOOST_MATH_GPU_ENABLED inline RealType logpdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
  112. {
  113. BOOST_MATH_STD_USING // for ADL of std function exp.
  114. const RealType sigma = dist.sigma();
  115. RealType result = -boost::math::numeric_limits<RealType>::infinity();
  116. constexpr auto function = "boost::math::logpdf(const rayleigh_distribution<%1%>&, %1%)";
  117. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  118. {
  119. return result;
  120. }
  121. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  122. {
  123. return result;
  124. }
  125. if((boost::math::isinf)(x))
  126. {
  127. return result;
  128. }
  129. result = -(x*x)/(2*sigma*sigma) - 2*log(sigma) + log(x);
  130. return result;
  131. } // logpdf
  132. template <class RealType, class Policy>
  133. BOOST_MATH_GPU_ENABLED inline RealType cdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
  134. {
  135. BOOST_MATH_STD_USING // for ADL of std functions
  136. RealType result = 0;
  137. RealType sigma = dist.sigma();
  138. constexpr auto function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
  139. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  140. {
  141. return result;
  142. }
  143. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  144. {
  145. return result;
  146. }
  147. result = -boost::math::expm1(-x * x / ( 2 * sigma * sigma), Policy());
  148. return result;
  149. } // cdf
  150. template <class RealType, class Policy>
  151. BOOST_MATH_GPU_ENABLED inline RealType logcdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
  152. {
  153. BOOST_MATH_STD_USING // for ADL of std functions
  154. RealType result = 0;
  155. RealType sigma = dist.sigma();
  156. constexpr auto function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
  157. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  158. {
  159. return -boost::math::numeric_limits<RealType>::infinity();
  160. }
  161. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  162. {
  163. return -boost::math::numeric_limits<RealType>::infinity();
  164. }
  165. result = log1p(-exp(-x * x / ( 2 * sigma * sigma)), Policy());
  166. return result;
  167. } // logcdf
  168. template <class RealType, class Policy>
  169. BOOST_MATH_GPU_ENABLED inline RealType quantile(const rayleigh_distribution<RealType, Policy>& dist, const RealType& p)
  170. {
  171. BOOST_MATH_STD_USING // for ADL of std functions
  172. RealType result = 0;
  173. RealType sigma = dist.sigma();
  174. constexpr auto function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
  175. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  176. return result;
  177. if(false == detail::check_probability(function, p, &result, Policy()))
  178. return result;
  179. if(p == 0)
  180. {
  181. return 0;
  182. }
  183. if(p == 1)
  184. {
  185. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  186. }
  187. result = sqrt(-2 * sigma * sigma * boost::math::log1p(-p, Policy()));
  188. return result;
  189. } // quantile
  190. template <class RealType, class Policy>
  191. BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
  192. {
  193. BOOST_MATH_STD_USING // for ADL of std functions
  194. RealType result = 0;
  195. RealType sigma = c.dist.sigma();
  196. constexpr auto function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
  197. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  198. {
  199. return result;
  200. }
  201. RealType x = c.param;
  202. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  203. {
  204. return result;
  205. }
  206. RealType ea = x * x / (2 * sigma * sigma);
  207. // Fix for VC11/12 x64 bug in exp(float):
  208. if (ea >= tools::max_value<RealType>())
  209. return 0;
  210. result = exp(-ea);
  211. return result;
  212. } // cdf complement
  213. template <class RealType, class Policy>
  214. BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
  215. {
  216. BOOST_MATH_STD_USING // for ADL of std functions
  217. RealType result = 0;
  218. RealType sigma = c.dist.sigma();
  219. constexpr auto function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
  220. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  221. {
  222. return -boost::math::numeric_limits<RealType>::infinity();
  223. }
  224. RealType x = c.param;
  225. if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
  226. {
  227. return -boost::math::numeric_limits<RealType>::infinity();
  228. }
  229. RealType ea = x * x / (2 * sigma * sigma);
  230. // Fix for VC11/12 x64 bug in exp(float):
  231. if (ea >= tools::max_value<RealType>())
  232. return 0;
  233. result = -ea;
  234. return result;
  235. } // logcdf complement
  236. template <class RealType, class Policy>
  237. BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
  238. {
  239. BOOST_MATH_STD_USING // for ADL of std functions, log & sqrt.
  240. RealType result = 0;
  241. RealType sigma = c.dist.sigma();
  242. constexpr auto function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
  243. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  244. {
  245. return result;
  246. }
  247. RealType q = c.param;
  248. if(false == detail::check_probability(function, q, &result, Policy()))
  249. {
  250. return result;
  251. }
  252. if(q == 1)
  253. {
  254. return 0;
  255. }
  256. if(q == 0)
  257. {
  258. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  259. }
  260. result = sqrt(-2 * sigma * sigma * log(q));
  261. return result;
  262. } // quantile complement
  263. template <class RealType, class Policy>
  264. BOOST_MATH_GPU_ENABLED inline RealType mean(const rayleigh_distribution<RealType, Policy>& dist)
  265. {
  266. RealType result = 0;
  267. RealType sigma = dist.sigma();
  268. constexpr auto function = "boost::math::mean(const rayleigh_distribution<%1%>&, %1%)";
  269. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  270. {
  271. return result;
  272. }
  273. using boost::math::constants::root_half_pi;
  274. return sigma * root_half_pi<RealType>();
  275. } // mean
  276. template <class RealType, class Policy>
  277. BOOST_MATH_GPU_ENABLED inline RealType variance(const rayleigh_distribution<RealType, Policy>& dist)
  278. {
  279. RealType result = 0;
  280. RealType sigma = dist.sigma();
  281. constexpr auto function = "boost::math::variance(const rayleigh_distribution<%1%>&, %1%)";
  282. if(false == detail::verify_sigma(function, sigma, &result, Policy()))
  283. {
  284. return result;
  285. }
  286. using boost::math::constants::four_minus_pi;
  287. return four_minus_pi<RealType>() * sigma * sigma / 2;
  288. } // variance
  289. template <class RealType, class Policy>
  290. BOOST_MATH_GPU_ENABLED inline RealType mode(const rayleigh_distribution<RealType, Policy>& dist)
  291. {
  292. return dist.sigma();
  293. }
  294. template <class RealType, class Policy>
  295. BOOST_MATH_GPU_ENABLED inline RealType median(const rayleigh_distribution<RealType, Policy>& dist)
  296. {
  297. using boost::math::constants::root_ln_four;
  298. return root_ln_four<RealType>() * dist.sigma();
  299. }
  300. template <class RealType, class Policy>
  301. BOOST_MATH_GPU_ENABLED inline RealType skewness(const rayleigh_distribution<RealType, Policy>& /*dist*/)
  302. {
  303. return static_cast<RealType>(0.63111065781893713819189935154422777984404221106391L);
  304. // Computed using NTL at 150 bit, about 50 decimal digits.
  305. // 2 * sqrt(pi) * (pi-3) / pow(4, 2/3) - pi
  306. }
  307. template <class RealType, class Policy>
  308. BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const rayleigh_distribution<RealType, Policy>& /*dist*/)
  309. {
  310. return static_cast<RealType>(3.2450893006876380628486604106197544154170667057995L);
  311. // Computed using NTL at 150 bit, about 50 decimal digits.
  312. // 3 - (6*pi*pi - 24*pi + 16) / pow(4-pi, 2)
  313. }
  314. template <class RealType, class Policy>
  315. BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const rayleigh_distribution<RealType, Policy>& /*dist*/)
  316. {
  317. return static_cast<RealType>(0.2450893006876380628486604106197544154170667057995L);
  318. // Computed using NTL at 150 bit, about 50 decimal digits.
  319. // -(6*pi*pi - 24*pi + 16) / pow(4-pi,2)
  320. } // kurtosis_excess
  321. template <class RealType, class Policy>
  322. BOOST_MATH_GPU_ENABLED inline RealType entropy(const rayleigh_distribution<RealType, Policy>& dist)
  323. {
  324. BOOST_MATH_STD_USING
  325. return 1 + log(dist.sigma()*constants::one_div_root_two<RealType>()) + constants::euler<RealType>()/2;
  326. }
  327. } // namespace math
  328. } // namespace boost
  329. #ifdef _MSC_VER
  330. # pragma warning(pop)
  331. #endif
  332. // This include must be at the end, *after* the accessors
  333. // for this distribution have been defined, in order to
  334. // keep compilers that support two-phase lookup happy.
  335. #include <boost/math/distributions/detail/derived_accessors.hpp>
  336. #endif // BOOST_STATS_rayleigh_HPP