382 lines · plain
1// Copyright Paul A. Bristow 2007.2// Copyright Matt Borland 2023.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_STATS_rayleigh_HPP8#define BOOST_STATS_rayleigh_HPP9 10#include <boost/math/tools/config.hpp>11#include <boost/math/tools/cstdint.hpp>12#include <boost/math/tools/tuple.hpp>13#include <boost/math/tools/numeric_limits.hpp>14#include <boost/math/distributions/fwd.hpp>15#include <boost/math/constants/constants.hpp>16#include <boost/math/special_functions/log1p.hpp>17#include <boost/math/special_functions/expm1.hpp>18#include <boost/math/distributions/complement.hpp>19#include <boost/math/distributions/detail/common_error_handling.hpp>20 21#ifdef _MSC_VER22# pragma warning(push)23# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).24#endif25 26namespace boost{ namespace math{27 28namespace detail29{ // Error checks:30 template <class RealType, class Policy>31 BOOST_MATH_GPU_ENABLED inline bool verify_sigma(const char* function, RealType sigma, RealType* presult, const Policy& pol)32 {33 if((sigma <= 0) || (!(boost::math::isfinite)(sigma)))34 {35 *presult = policies::raise_domain_error<RealType>(36 function,37 "The scale parameter \"sigma\" must be > 0 and finite, but was: %1%.", sigma, pol);38 return false;39 }40 return true;41 } // bool verify_sigma42 43 template <class RealType, class Policy>44 BOOST_MATH_GPU_ENABLED inline bool verify_rayleigh_x(const char* function, RealType x, RealType* presult, const Policy& pol)45 {46 if((x < 0) || (boost::math::isnan)(x))47 {48 *presult = policies::raise_domain_error<RealType>(49 function,50 "The random variable must be >= 0, but was: %1%.", x, pol);51 return false;52 }53 return true;54 } // bool verify_rayleigh_x55} // namespace detail56 57template <class RealType = double, class Policy = policies::policy<> >58class rayleigh_distribution59{60public:61 using value_type = RealType;62 using policy_type = Policy;63 64 BOOST_MATH_GPU_ENABLED explicit rayleigh_distribution(RealType l_sigma = 1)65 : m_sigma(l_sigma)66 {67 RealType err;68 detail::verify_sigma("boost::math::rayleigh_distribution<%1%>::rayleigh_distribution", l_sigma, &err, Policy());69 } // rayleigh_distribution70 71 BOOST_MATH_GPU_ENABLED RealType sigma()const72 { // Accessor.73 return m_sigma;74 }75 76private:77 RealType m_sigma;78}; // class rayleigh_distribution79 80using rayleigh = rayleigh_distribution<double>;81 82#ifdef __cpp_deduction_guides83template <class RealType>84rayleigh_distribution(RealType)->rayleigh_distribution<typename boost::math::tools::promote_args<RealType>::type>;85#endif86 87template <class RealType, class Policy>88BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const rayleigh_distribution<RealType, Policy>& /*dist*/)89{ // Range of permissible values for random variable x.90 using boost::math::tools::max_value;91 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>());92}93 94template <class RealType, class Policy>95BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const rayleigh_distribution<RealType, Policy>& /*dist*/)96{ // Range of supported values for random variable x.97 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.98 using boost::math::tools::max_value;99 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());100}101 102template <class RealType, class Policy>103BOOST_MATH_GPU_ENABLED inline RealType pdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)104{105 BOOST_MATH_STD_USING // for ADL of std function exp.106 107 RealType sigma = dist.sigma();108 RealType result = 0;109 constexpr auto function = "boost::math::pdf(const rayleigh_distribution<%1%>&, %1%)";110 if(false == detail::verify_sigma(function, sigma, &result, Policy()))111 {112 return result;113 }114 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))115 {116 return result;117 }118 if((boost::math::isinf)(x))119 {120 return 0;121 }122 RealType sigmasqr = sigma * sigma;123 result = x * (exp(-(x * x) / ( 2 * sigmasqr))) / sigmasqr;124 return result;125} // pdf126 127template <class RealType, class Policy>128BOOST_MATH_GPU_ENABLED inline RealType logpdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)129{130 BOOST_MATH_STD_USING // for ADL of std function exp.131 132 const RealType sigma = dist.sigma();133 RealType result = -boost::math::numeric_limits<RealType>::infinity();134 constexpr auto function = "boost::math::logpdf(const rayleigh_distribution<%1%>&, %1%)";135 136 if(false == detail::verify_sigma(function, sigma, &result, Policy()))137 {138 return result;139 }140 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))141 {142 return result;143 }144 if((boost::math::isinf)(x))145 {146 return result;147 }148 149 result = -(x*x)/(2*sigma*sigma) - 2*log(sigma) + log(x);150 return result;151} // logpdf152 153template <class RealType, class Policy>154BOOST_MATH_GPU_ENABLED inline RealType cdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)155{156 BOOST_MATH_STD_USING // for ADL of std functions157 158 RealType result = 0;159 RealType sigma = dist.sigma();160 constexpr auto function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";161 if(false == detail::verify_sigma(function, sigma, &result, Policy()))162 {163 return result;164 }165 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))166 {167 return result;168 }169 result = -boost::math::expm1(-x * x / ( 2 * sigma * sigma), Policy());170 return result;171} // cdf172 173template <class RealType, class Policy>174BOOST_MATH_GPU_ENABLED inline RealType logcdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)175{176 BOOST_MATH_STD_USING // for ADL of std functions177 178 RealType result = 0;179 RealType sigma = dist.sigma();180 constexpr auto function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";181 if(false == detail::verify_sigma(function, sigma, &result, Policy()))182 {183 return -boost::math::numeric_limits<RealType>::infinity();184 }185 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))186 {187 return -boost::math::numeric_limits<RealType>::infinity();188 }189 result = log1p(-exp(-x * x / ( 2 * sigma * sigma)), Policy()); 190 return result;191} // logcdf192 193template <class RealType, class Policy>194BOOST_MATH_GPU_ENABLED inline RealType quantile(const rayleigh_distribution<RealType, Policy>& dist, const RealType& p)195{196 BOOST_MATH_STD_USING // for ADL of std functions197 198 RealType result = 0;199 RealType sigma = dist.sigma();200 constexpr auto function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";201 if(false == detail::verify_sigma(function, sigma, &result, Policy()))202 return result;203 if(false == detail::check_probability(function, p, &result, Policy()))204 return result;205 206 if(p == 0)207 {208 return 0;209 }210 if(p == 1)211 {212 return policies::raise_overflow_error<RealType>(function, 0, Policy());213 }214 result = sqrt(-2 * sigma * sigma * boost::math::log1p(-p, Policy()));215 return result;216} // quantile217 218template <class RealType, class Policy>219BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)220{221 BOOST_MATH_STD_USING // for ADL of std functions222 223 RealType result = 0;224 RealType sigma = c.dist.sigma();225 constexpr auto function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";226 if(false == detail::verify_sigma(function, sigma, &result, Policy()))227 {228 return result;229 }230 RealType x = c.param;231 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))232 {233 return result;234 }235 RealType ea = x * x / (2 * sigma * sigma);236 // Fix for VC11/12 x64 bug in exp(float):237 if (ea >= tools::max_value<RealType>())238 return 0;239 result = exp(-ea);240 return result;241} // cdf complement242 243template <class RealType, class Policy>244BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)245{246 BOOST_MATH_STD_USING // for ADL of std functions247 248 RealType result = 0;249 RealType sigma = c.dist.sigma();250 constexpr auto function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";251 if(false == detail::verify_sigma(function, sigma, &result, Policy()))252 {253 return -boost::math::numeric_limits<RealType>::infinity();254 }255 RealType x = c.param;256 if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))257 {258 return -boost::math::numeric_limits<RealType>::infinity();259 }260 RealType ea = x * x / (2 * sigma * sigma);261 // Fix for VC11/12 x64 bug in exp(float):262 if (ea >= tools::max_value<RealType>())263 return 0;264 result = -ea;265 return result;266} // logcdf complement267 268template <class RealType, class Policy>269BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)270{271 BOOST_MATH_STD_USING // for ADL of std functions, log & sqrt.272 273 RealType result = 0;274 RealType sigma = c.dist.sigma();275 constexpr auto function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";276 if(false == detail::verify_sigma(function, sigma, &result, Policy()))277 {278 return result;279 }280 RealType q = c.param;281 if(false == detail::check_probability(function, q, &result, Policy()))282 {283 return result;284 }285 if(q == 1)286 {287 return 0;288 }289 if(q == 0)290 {291 return policies::raise_overflow_error<RealType>(function, 0, Policy());292 }293 result = sqrt(-2 * sigma * sigma * log(q));294 return result;295} // quantile complement296 297template <class RealType, class Policy>298BOOST_MATH_GPU_ENABLED inline RealType mean(const rayleigh_distribution<RealType, Policy>& dist)299{300 RealType result = 0;301 RealType sigma = dist.sigma();302 constexpr auto function = "boost::math::mean(const rayleigh_distribution<%1%>&, %1%)";303 if(false == detail::verify_sigma(function, sigma, &result, Policy()))304 {305 return result;306 }307 using boost::math::constants::root_half_pi;308 return sigma * root_half_pi<RealType>();309} // mean310 311template <class RealType, class Policy>312BOOST_MATH_GPU_ENABLED inline RealType variance(const rayleigh_distribution<RealType, Policy>& dist)313{314 RealType result = 0;315 RealType sigma = dist.sigma();316 constexpr auto function = "boost::math::variance(const rayleigh_distribution<%1%>&, %1%)";317 if(false == detail::verify_sigma(function, sigma, &result, Policy()))318 {319 return result;320 }321 using boost::math::constants::four_minus_pi;322 return four_minus_pi<RealType>() * sigma * sigma / 2;323} // variance324 325template <class RealType, class Policy>326BOOST_MATH_GPU_ENABLED inline RealType mode(const rayleigh_distribution<RealType, Policy>& dist)327{328 return dist.sigma();329}330 331template <class RealType, class Policy>332BOOST_MATH_GPU_ENABLED inline RealType median(const rayleigh_distribution<RealType, Policy>& dist)333{334 using boost::math::constants::root_ln_four;335 return root_ln_four<RealType>() * dist.sigma();336}337 338template <class RealType, class Policy>339BOOST_MATH_GPU_ENABLED inline RealType skewness(const rayleigh_distribution<RealType, Policy>& /*dist*/)340{341 return static_cast<RealType>(0.63111065781893713819189935154422777984404221106391L);342 // Computed using NTL at 150 bit, about 50 decimal digits.343 // 2 * sqrt(pi) * (pi-3) / pow(4, 2/3) - pi344}345 346template <class RealType, class Policy>347BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const rayleigh_distribution<RealType, Policy>& /*dist*/)348{349 return static_cast<RealType>(3.2450893006876380628486604106197544154170667057995L);350 // Computed using NTL at 150 bit, about 50 decimal digits.351 // 3 - (6*pi*pi - 24*pi + 16) / pow(4-pi, 2)352}353 354template <class RealType, class Policy>355BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const rayleigh_distribution<RealType, Policy>& /*dist*/)356{357 return static_cast<RealType>(0.2450893006876380628486604106197544154170667057995L);358 // Computed using NTL at 150 bit, about 50 decimal digits.359 // -(6*pi*pi - 24*pi + 16) / pow(4-pi,2)360} // kurtosis_excess361 362template <class RealType, class Policy>363BOOST_MATH_GPU_ENABLED inline RealType entropy(const rayleigh_distribution<RealType, Policy>& dist)364{365 BOOST_MATH_STD_USING366 return 1 + log(dist.sigma()*constants::one_div_root_two<RealType>()) + constants::euler<RealType>()/2;367}368 369} // namespace math370} // namespace boost371 372#ifdef _MSC_VER373# pragma warning(pop)374#endif375 376// This include must be at the end, *after* the accessors377// for this distribution have been defined, in order to378// keep compilers that support two-phase lookup happy.379#include <boost/math/distributions/detail/derived_accessors.hpp>380 381#endif // BOOST_STATS_rayleigh_HPP382