math/distributions/rayleigh.hpp

0001 //  Copyright Paul A. Bristow 2007.
0002 //  Copyright Matt Borland 2023.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_STATS_rayleigh_HPP
0008 #define BOOST_STATS_rayleigh_HPP
0009
0010 #include <boost/math/distributions/fwd.hpp>
0011 #include <boost/math/constants/constants.hpp>
0012 #include <boost/math/special_functions/log1p.hpp>
0013 #include <boost/math/special_functions/expm1.hpp>
0014 #include <boost/math/distributions/complement.hpp>
0015 #include <boost/math/distributions/detail/common_error_handling.hpp>
0016
0017 #ifdef _MSC_VER
0018 # pragma warning(push)
0019 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
0020 #endif
0021
0022 #include <utility>
0023 #include <limits>
0024 #include <cmath>
0025
0026 namespace boost{ namespace math{
0027
0028 namespace detail
0029 { // Error checks:
0030   template <class RealType, class Policy>
0031   inline bool verify_sigma(const char* function, RealType sigma, RealType* presult, const Policy& pol)
0032   {
0033      if((sigma <= 0) || (!(boost::math::isfinite)(sigma)))
0034      {
0035         *presult = policies::raise_domain_error<RealType>(
0036            function,
0037            "The scale parameter \"sigma\" must be > 0 and finite, but was: %1%.", sigma, pol);
0038         return false;
0039      }
0040      return true;
0041   } // bool verify_sigma
0042
0043   template <class RealType, class Policy>
0044   inline bool verify_rayleigh_x(const char* function, RealType x, RealType* presult, const Policy& pol)
0045   {
0046      if((x < 0) || (boost::math::isnan)(x))
0047      {
0048         *presult = policies::raise_domain_error<RealType>(
0049            function,
0050            "The random variable must be >= 0, but was: %1%.", x, pol);
0051         return false;
0052      }
0053      return true;
0054   } // bool verify_rayleigh_x
0055 } // namespace detail
0056
0057 template <class RealType = double, class Policy = policies::policy<> >
0058 class rayleigh_distribution
0059 {
0060 public:
0061    using value_type = RealType;
0062    using policy_type = Policy;
0063
0064    explicit rayleigh_distribution(RealType l_sigma = 1)
0065       : m_sigma(l_sigma)
0066    {
0067       RealType err;
0068       detail::verify_sigma("boost::math::rayleigh_distribution<%1%>::rayleigh_distribution", l_sigma, &err, Policy());
0069    } // rayleigh_distribution
0070
0071    RealType sigma()const
0072    { // Accessor.
0073      return m_sigma;
0074    }
0075
0076 private:
0077    RealType m_sigma;
0078 }; // class rayleigh_distribution
0079
0080 using rayleigh = rayleigh_distribution<double>;
0081
0082 #ifdef __cpp_deduction_guides
0083 template <class RealType>
0084 rayleigh_distribution(RealType)->rayleigh_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0085 #endif
0086
0087 template <class RealType, class Policy>
0088 inline std::pair<RealType, RealType> range(const rayleigh_distribution<RealType, Policy>& /*dist*/)
0089 { // Range of permissible values for random variable x.
0090    using boost::math::tools::max_value;
0091    return std::pair<RealType, RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>());
0092 }
0093
0094 template <class RealType, class Policy>
0095 inline std::pair<RealType, RealType> support(const rayleigh_distribution<RealType, Policy>& /*dist*/)
0096 { // Range of supported values for random variable x.
0097    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0098    using boost::math::tools::max_value;
0099    return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
0100 }
0101
0102 template <class RealType, class Policy>
0103 inline RealType pdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
0104 {
0105    BOOST_MATH_STD_USING // for ADL of std function exp.
0106
0107    RealType sigma = dist.sigma();
0108    RealType result = 0;
0109    static const char* function = "boost::math::pdf(const rayleigh_distribution<%1%>&, %1%)";
0110    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0111    {
0112       return result;
0113    }
0114    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0115    {
0116       return result;
0117    }
0118    if((boost::math::isinf)(x))
0119    {
0120       return 0;
0121    }
0122    RealType sigmasqr = sigma * sigma;
0123    result = x * (exp(-(x * x) / ( 2 * sigmasqr))) / sigmasqr;
0124    return result;
0125 } // pdf
0126
0127 template <class RealType, class Policy>
0128 inline RealType logpdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
0129 {
0130    BOOST_MATH_STD_USING // for ADL of std function exp.
0131
0132    const RealType sigma = dist.sigma();
0133    RealType result = -std::numeric_limits<RealType>::infinity();
0134    static const char* function = "boost::math::logpdf(const rayleigh_distribution<%1%>&, %1%)";
0135
0136    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0137    {
0138       return result;
0139    }
0140    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0141    {
0142       return result;
0143    }
0144    if((boost::math::isinf)(x))
0145    {
0146       return result;
0147    }
0148
0149    result = -(x*x)/(2*sigma*sigma) - 2*log(sigma) + log(x);
0150    return result;
0151 } // logpdf
0152
0153 template <class RealType, class Policy>
0154 inline RealType cdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
0155 {
0156    BOOST_MATH_STD_USING // for ADL of std functions
0157
0158    RealType result = 0;
0159    RealType sigma = dist.sigma();
0160    static const char* function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
0161    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0162    {
0163       return result;
0164    }
0165    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0166    {
0167       return result;
0168    }
0169    result = -boost::math::expm1(-x * x / ( 2 * sigma * sigma), Policy());
0170    return result;
0171 } // cdf
0172
0173 template <class RealType, class Policy>
0174 inline RealType logcdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
0175 {
0176    BOOST_MATH_STD_USING // for ADL of std functions
0177
0178    RealType result = 0;
0179    RealType sigma = dist.sigma();
0180    static const char* function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
0181    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0182    {
0183       return -std::numeric_limits<RealType>::infinity();
0184    }
0185    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0186    {
0187       return -std::numeric_limits<RealType>::infinity();
0188    }
0189    result = log1p(-exp(-x * x / ( 2 * sigma * sigma)), Policy());
0190    return result;
0191 } // logcdf
0192
0193 template <class RealType, class Policy>
0194 inline RealType quantile(const rayleigh_distribution<RealType, Policy>& dist, const RealType& p)
0195 {
0196    BOOST_MATH_STD_USING // for ADL of std functions
0197
0198    RealType result = 0;
0199    RealType sigma = dist.sigma();
0200    static const char* function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
0201    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0202       return result;
0203    if(false == detail::check_probability(function, p, &result, Policy()))
0204       return result;
0205
0206    if(p == 0)
0207    {
0208       return 0;
0209    }
0210    if(p == 1)
0211    {
0212      return policies::raise_overflow_error<RealType>(function, 0, Policy());
0213    }
0214    result = sqrt(-2 * sigma * sigma * boost::math::log1p(-p, Policy()));
0215    return result;
0216 } // quantile
0217
0218 template <class RealType, class Policy>
0219 inline RealType cdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
0220 {
0221    BOOST_MATH_STD_USING // for ADL of std functions
0222
0223    RealType result = 0;
0224    RealType sigma = c.dist.sigma();
0225    static const char* function = "boost::math::cdf(const rayleigh_distribution<%1%>&, %1%)";
0226    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0227    {
0228       return result;
0229    }
0230    RealType x = c.param;
0231    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0232    {
0233       return result;
0234    }
0235    RealType ea = x * x / (2 * sigma * sigma);
0236    // Fix for VC11/12 x64 bug in exp(float):
0237    if (ea >= tools::max_value<RealType>())
0238       return 0;
0239    result =  exp(-ea);
0240    return result;
0241 } // cdf complement
0242
0243 template <class RealType, class Policy>
0244 inline RealType logcdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
0245 {
0246    BOOST_MATH_STD_USING // for ADL of std functions
0247
0248    RealType result = 0;
0249    RealType sigma = c.dist.sigma();
0250    static const char* function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
0251    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0252    {
0253       return -std::numeric_limits<RealType>::infinity();
0254    }
0255    RealType x = c.param;
0256    if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
0257    {
0258       return -std::numeric_limits<RealType>::infinity();
0259    }
0260    RealType ea = x * x / (2 * sigma * sigma);
0261    // Fix for VC11/12 x64 bug in exp(float):
0262    if (ea >= tools::max_value<RealType>())
0263       return 0;
0264    result = -ea;
0265    return result;
0266 } // logcdf complement
0267
0268 template <class RealType, class Policy>
0269 inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
0270 {
0271    BOOST_MATH_STD_USING // for ADL of std functions, log & sqrt.
0272
0273    RealType result = 0;
0274    RealType sigma = c.dist.sigma();
0275    static const char* function = "boost::math::quantile(const rayleigh_distribution<%1%>&, %1%)";
0276    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0277    {
0278       return result;
0279    }
0280    RealType q = c.param;
0281    if(false == detail::check_probability(function, q, &result, Policy()))
0282    {
0283       return result;
0284    }
0285    if(q == 1)
0286    {
0287       return 0;
0288    }
0289    if(q == 0)
0290    {
0291      return policies::raise_overflow_error<RealType>(function, 0, Policy());
0292    }
0293    result = sqrt(-2 * sigma * sigma * log(q));
0294    return result;
0295 } // quantile complement
0296
0297 template <class RealType, class Policy>
0298 inline RealType mean(const rayleigh_distribution<RealType, Policy>& dist)
0299 {
0300    RealType result = 0;
0301    RealType sigma = dist.sigma();
0302    static const char* function = "boost::math::mean(const rayleigh_distribution<%1%>&, %1%)";
0303    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0304    {
0305       return result;
0306    }
0307    using boost::math::constants::root_half_pi;
0308    return sigma * root_half_pi<RealType>();
0309 } // mean
0310
0311 template <class RealType, class Policy>
0312 inline RealType variance(const rayleigh_distribution<RealType, Policy>& dist)
0313 {
0314    RealType result = 0;
0315    RealType sigma = dist.sigma();
0316    static const char* function = "boost::math::variance(const rayleigh_distribution<%1%>&, %1%)";
0317    if(false == detail::verify_sigma(function, sigma, &result, Policy()))
0318    {
0319       return result;
0320    }
0321    using boost::math::constants::four_minus_pi;
0322    return four_minus_pi<RealType>() * sigma * sigma / 2;
0323 } // variance
0324
0325 template <class RealType, class Policy>
0326 inline RealType mode(const rayleigh_distribution<RealType, Policy>& dist)
0327 {
0328    return dist.sigma();
0329 }
0330
0331 template <class RealType, class Policy>
0332 inline RealType median(const rayleigh_distribution<RealType, Policy>& dist)
0333 {
0334    using boost::math::constants::root_ln_four;
0335    return root_ln_four<RealType>() * dist.sigma();
0336 }
0337
0338 template <class RealType, class Policy>
0339 inline RealType skewness(const rayleigh_distribution<RealType, Policy>& /*dist*/)
0340 {
0341   return static_cast<RealType>(0.63111065781893713819189935154422777984404221106391L);
0342   // Computed using NTL at 150 bit, about 50 decimal digits.
0343   // 2 * sqrt(pi) * (pi-3) / pow(4, 2/3) - pi
0344 }
0345
0346 template <class RealType, class Policy>
0347 inline RealType kurtosis(const rayleigh_distribution<RealType, Policy>& /*dist*/)
0348 {
0349   return static_cast<RealType>(3.2450893006876380628486604106197544154170667057995L);
0350   // Computed using NTL at 150 bit, about 50 decimal digits.
0351   // 3 - (6*pi*pi - 24*pi + 16) / pow(4-pi, 2)
0352 }
0353
0354 template <class RealType, class Policy>
0355 inline RealType kurtosis_excess(const rayleigh_distribution<RealType, Policy>& /*dist*/)
0356 {
0357   return static_cast<RealType>(0.2450893006876380628486604106197544154170667057995L);
0358   // Computed using NTL at 150 bit, about 50 decimal digits.
0359   // -(6*pi*pi - 24*pi + 16) / pow(4-pi,2)
0360 } // kurtosis_excess
0361
0362 template <class RealType, class Policy>
0363 inline RealType entropy(const rayleigh_distribution<RealType, Policy>& dist)
0364 {
0365    using std::log;
0366    return 1 + log(dist.sigma()*constants::one_div_root_two<RealType>()) + constants::euler<RealType>()/2;
0367 }
0368
0369 } // namespace math
0370 } // namespace boost
0371
0372 #ifdef _MSC_VER
0373 # pragma warning(pop)
0374 #endif
0375
0376 // This include must be at the end, *after* the accessors
0377 // for this distribution have been defined, in order to
0378 // keep compilers that support two-phase lookup happy.
0379 #include <boost/math/distributions/detail/derived_accessors.hpp>
0380
0381 #endif // BOOST_STATS_rayleigh_HPP