Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-17 08:35:33 

0001 //  Copyright John Maddock 2007.
0002 //  Copyright Paul A. Bristow 2007, 2009
0003 //  Copyright Matt Borland 2023.
0004 //  Use, modification and distribution are subject to the
0005 //  Boost Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007 
0008 #ifndef BOOST_STATS_PARETO_HPP
0009 #define BOOST_STATS_PARETO_HPP
0010 
0011 // http://en.wikipedia.org/wiki/Pareto_distribution
0012 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
0013 // Also:
0014 // Weisstein, Eric W. "Pareto Distribution."
0015 // From MathWorld--A Wolfram Web Resource.
0016 // http://mathworld.wolfram.com/ParetoDistribution.html
0017 // Handbook of Statistical Distributions with Applications, K Krishnamoorthy, ISBN 1-58488-635-8, Chapter 23, pp 257 - 267.
0018 // Caution KK's a and b are the reverse of Mathworld!
0019 
0020 #include <boost/math/tools/config.hpp>
0021 #include <boost/math/tools/tuple.hpp>
0022 #include <boost/math/tools/numeric_limits.hpp>
0023 #include <boost/math/distributions/fwd.hpp>
0024 #include <boost/math/distributions/complement.hpp>
0025 #include <boost/math/distributions/detail/common_error_handling.hpp>
0026 #include <boost/math/special_functions/powm1.hpp>
0027 #include <boost/math/special_functions/log1p.hpp>
0028 
0029 namespace boost
0030 {
0031   namespace math
0032   {
0033     namespace detail
0034     { // Parameter checking.
0035       template <class RealType, class Policy>
0036       BOOST_MATH_GPU_ENABLED inline bool check_pareto_scale(
0037         const char* function,
0038         RealType scale,
0039         RealType* result, const Policy& pol)
0040       {
0041         if((boost::math::isfinite)(scale))
0042         { // any > 0 finite value is OK.
0043           if (scale > 0)
0044           {
0045             return true;
0046           }
0047           else
0048           {
0049             *result = policies::raise_domain_error<RealType>(
0050               function,
0051               "Scale parameter is %1%, but must be > 0!", scale, pol);
0052             return false;
0053           }
0054         }
0055         else
0056         { // Not finite.
0057           *result = policies::raise_domain_error<RealType>(
0058             function,
0059             "Scale parameter is %1%, but must be finite!", scale, pol);
0060           return false;
0061         }
0062       } // bool check_pareto_scale
0063 
0064       template <class RealType, class Policy>
0065       BOOST_MATH_GPU_ENABLED inline bool check_pareto_shape(
0066         const char* function,
0067         RealType shape,
0068         RealType* result, const Policy& pol)
0069       {
0070         if((boost::math::isfinite)(shape))
0071         { // Any finite value > 0 is OK.
0072           if (shape > 0)
0073           {
0074             return true;
0075           }
0076           else
0077           {
0078             *result = policies::raise_domain_error<RealType>(
0079               function,
0080               "Shape parameter is %1%, but must be > 0!", shape, pol);
0081             return false;
0082           }
0083         }
0084         else
0085         { // Not finite.
0086           *result = policies::raise_domain_error<RealType>(
0087             function,
0088             "Shape parameter is %1%, but must be finite!", shape, pol);
0089           return false;
0090         }
0091       } // bool check_pareto_shape(
0092 
0093       template <class RealType, class Policy>
0094       BOOST_MATH_GPU_ENABLED inline bool check_pareto_x(
0095         const char* function,
0096         RealType const& x,
0097         RealType* result, const Policy& pol)
0098       {
0099         if((boost::math::isfinite)(x))
0100         { //
0101           if (x > 0)
0102           {
0103             return true;
0104           }
0105           else
0106           {
0107             *result = policies::raise_domain_error<RealType>(
0108               function,
0109               "x parameter is %1%, but must be > 0 !", x, pol);
0110             return false;
0111           }
0112         }
0113         else
0114         { // Not finite..
0115           *result = policies::raise_domain_error<RealType>(
0116             function,
0117             "x parameter is %1%, but must be finite!", x, pol);
0118           return false;
0119         }
0120       } // bool check_pareto_x
0121 
0122       template <class RealType, class Policy>
0123       BOOST_MATH_GPU_ENABLED inline bool check_pareto( // distribution parameters.
0124         const char* function,
0125         RealType scale,
0126         RealType shape,
0127         RealType* result, const Policy& pol)
0128       {
0129         return check_pareto_scale(function, scale, result, pol)
0130            && check_pareto_shape(function, shape, result, pol);
0131       } // bool check_pareto(
0132 
0133     } // namespace detail
0134 
0135     template <class RealType = double, class Policy = policies::policy<> >
0136     class pareto_distribution
0137     {
0138     public:
0139       typedef RealType value_type;
0140       typedef Policy policy_type;
0141 
0142       BOOST_MATH_GPU_ENABLED pareto_distribution(RealType l_scale = 1, RealType l_shape = 1)
0143         : m_scale(l_scale), m_shape(l_shape)
0144       { // Constructor.
0145         RealType result = 0;
0146         detail::check_pareto("boost::math::pareto_distribution<%1%>::pareto_distribution", l_scale, l_shape, &result, Policy());
0147       }
0148 
0149       BOOST_MATH_GPU_ENABLED RealType scale()const
0150       { // AKA Xm and Wolfram b and beta
0151         return m_scale;
0152       }
0153 
0154       BOOST_MATH_GPU_ENABLED RealType shape()const
0155       { // AKA k and Wolfram a and alpha
0156         return m_shape;
0157       }
0158     private:
0159       // Data members:
0160       RealType m_scale;  // distribution scale (xm) or beta
0161       RealType m_shape;  // distribution shape (k) or alpha
0162     };
0163 
0164     typedef pareto_distribution<double> pareto; // Convenience to allow pareto(2., 3.);
0165 
0166     #ifdef __cpp_deduction_guides
0167     template <class RealType>
0168     pareto_distribution(RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0169     template <class RealType>
0170     pareto_distribution(RealType,RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0171     #endif
0172 
0173 
0174     template <class RealType, class Policy>
0175     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const pareto_distribution<RealType, Policy>& /*dist*/)
0176     { // Range of permissible values for random variable x.
0177       using boost::math::tools::max_value;
0178       return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // scale zero to + infinity.
0179     } // range
0180 
0181     template <class RealType, class Policy>
0182     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const pareto_distribution<RealType, Policy>& dist)
0183     { // Range of supported values for random variable x.
0184       // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0185       using boost::math::tools::max_value;
0186       return boost::math::pair<RealType, RealType>(dist.scale(), max_value<RealType>() ); // scale to + infinity.
0187     } // support
0188 
0189     template <class RealType, class Policy>
0190     BOOST_MATH_GPU_ENABLED inline RealType pdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
0191     {
0192       BOOST_MATH_STD_USING  // for ADL of std function pow.
0193       constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
0194       RealType scale = dist.scale();
0195       RealType shape = dist.shape();
0196       RealType result = 0;
0197       if(false == (detail::check_pareto_x(function, x, &result, Policy())
0198          && detail::check_pareto(function, scale, shape, &result, Policy())))
0199          return result;
0200       if (x < scale)
0201       { // regardless of shape, pdf is zero (or should be disallow x < scale and throw an exception?).
0202         return 0;
0203       }
0204       result = shape * pow(scale, shape) / pow(x, shape+1);
0205       return result;
0206     } // pdf
0207 
0208     template <class RealType, class Policy>
0209     BOOST_MATH_GPU_ENABLED inline RealType cdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
0210     {
0211       BOOST_MATH_STD_USING  // for ADL of std function pow.
0212       constexpr auto function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";
0213       RealType scale = dist.scale();
0214       RealType shape = dist.shape();
0215       RealType result = 0;
0216 
0217       if(false == (detail::check_pareto_x(function, x, &result, Policy())
0218          && detail::check_pareto(function, scale, shape, &result, Policy())))
0219          return result;
0220 
0221       if (x <= scale)
0222       { // regardless of shape, cdf is zero.
0223         return 0;
0224       }
0225 
0226       // result = RealType(1) - pow((scale / x), shape);
0227       result = -boost::math::powm1(scale/x, shape, Policy()); // should be more accurate.
0228       return result;
0229     } // cdf
0230 
0231     template <class RealType, class Policy>
0232     BOOST_MATH_GPU_ENABLED inline RealType logcdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
0233     {
0234       BOOST_MATH_STD_USING  // for ADL of std function pow.
0235       constexpr auto function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";
0236       RealType scale = dist.scale();
0237       RealType shape = dist.shape();
0238       RealType result = 0;
0239 
0240       if(false == (detail::check_pareto_x(function, x, &result, Policy())
0241          && detail::check_pareto(function, scale, shape, &result, Policy())))
0242          return result;
0243 
0244       if (x <= scale)
0245       { // regardless of shape, cdf is zero.
0246         return -boost::math::numeric_limits<RealType>::infinity();
0247       }
0248 
0249       result = log1p(-pow(scale/x, shape), Policy());
0250       return result;
0251     } // logcdf
0252 
0253     template <class RealType, class Policy>
0254     BOOST_MATH_GPU_ENABLED inline RealType quantile(const pareto_distribution<RealType, Policy>& dist, const RealType& p)
0255     {
0256       BOOST_MATH_STD_USING  // for ADL of std function pow.
0257       constexpr auto function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";
0258       RealType result = 0;
0259       RealType scale = dist.scale();
0260       RealType shape = dist.shape();
0261       if(false == (detail::check_probability(function, p, &result, Policy())
0262            && detail::check_pareto(function, scale, shape, &result, Policy())))
0263       {
0264         return result;
0265       }
0266       if (p == 0)
0267       {
0268         return scale; // x must be scale (or less).
0269       }
0270       if (p == 1)
0271       {
0272         return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.
0273       }
0274       result = scale /
0275         (pow((1 - p), 1 / shape));
0276       // K. Krishnamoorthy,  ISBN 1-58488-635-8 eq 23.1.3
0277       return result;
0278     } // quantile
0279 
0280     template <class RealType, class Policy>
0281     BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
0282     {
0283        BOOST_MATH_STD_USING  // for ADL of std function pow.
0284        constexpr auto function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";
0285        RealType result = 0;
0286        RealType x = c.param;
0287        RealType scale = c.dist.scale();
0288        RealType shape = c.dist.shape();
0289        if(false == (detail::check_pareto_x(function, x, &result, Policy())
0290            && detail::check_pareto(function, scale, shape, &result, Policy())))
0291          return result;
0292 
0293        if (x <= scale)
0294        { // regardless of shape, cdf is zero, and complement is unity.
0295          return 1;
0296        }
0297        result = pow((scale/x), shape);
0298 
0299        return result;
0300     } // cdf complement
0301 
0302     template <class RealType, class Policy>
0303     BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
0304     {
0305        BOOST_MATH_STD_USING  // for ADL of std function pow.
0306        constexpr auto function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";
0307        RealType result = 0;
0308        RealType x = c.param;
0309        RealType scale = c.dist.scale();
0310        RealType shape = c.dist.shape();
0311        if(false == (detail::check_pareto_x(function, x, &result, Policy())
0312            && detail::check_pareto(function, scale, shape, &result, Policy())))
0313          return result;
0314 
0315        if (x <= scale)
0316        { // regardless of shape, cdf is zero, and complement is unity.
0317          return 0;
0318        }
0319        result = log(pow((scale/x), shape));
0320 
0321        return result;
0322     } // logcdf complement
0323 
0324     template <class RealType, class Policy>
0325     BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
0326     {
0327       BOOST_MATH_STD_USING  // for ADL of std function pow.
0328       constexpr auto function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";
0329       RealType result = 0;
0330       RealType q = c.param;
0331       RealType scale = c.dist.scale();
0332       RealType shape = c.dist.shape();
0333       if(false == (detail::check_probability(function, q, &result, Policy())
0334            && detail::check_pareto(function, scale, shape, &result, Policy())))
0335       {
0336         return result;
0337       }
0338       if (q == 1)
0339       {
0340         return scale; // x must be scale (or less).
0341       }
0342       if (q == 0)
0343       {
0344          return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.
0345       }
0346       result = scale / (pow(q, 1 / shape));
0347       // K. Krishnamoorthy,  ISBN 1-58488-635-8 eq 23.1.3
0348       return result;
0349     } // quantile complement
0350 
0351     template <class RealType, class Policy>
0352     BOOST_MATH_GPU_ENABLED inline RealType mean(const pareto_distribution<RealType, Policy>& dist)
0353     {
0354       RealType result = 0;
0355       constexpr auto function = "boost::math::mean(const pareto_distribution<%1%>&, %1%)";
0356       if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))
0357       {
0358         return result;
0359       }
0360       if (dist.shape() > RealType(1))
0361       {
0362         return dist.shape() * dist.scale() / (dist.shape() - 1);
0363       }
0364       else
0365       {
0366         using boost::math::tools::max_value;
0367         return max_value<RealType>(); // +infinity.
0368       }
0369     } // mean
0370 
0371     template <class RealType, class Policy>
0372     BOOST_MATH_GPU_ENABLED inline RealType mode(const pareto_distribution<RealType, Policy>& dist)
0373     {
0374       return dist.scale();
0375     } // mode
0376 
0377     template <class RealType, class Policy>
0378     BOOST_MATH_GPU_ENABLED inline RealType median(const pareto_distribution<RealType, Policy>& dist)
0379     {
0380       RealType result = 0;
0381       constexpr auto function = "boost::math::median(const pareto_distribution<%1%>&, %1%)";
0382       if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))
0383       {
0384         return result;
0385       }
0386       BOOST_MATH_STD_USING
0387       return dist.scale() * pow(RealType(2), (1/dist.shape()));
0388     } // median
0389 
0390     template <class RealType, class Policy>
0391     BOOST_MATH_GPU_ENABLED inline RealType variance(const pareto_distribution<RealType, Policy>& dist)
0392     {
0393       RealType result = 0;
0394       RealType scale = dist.scale();
0395       RealType shape = dist.shape();
0396       constexpr auto function = "boost::math::variance(const pareto_distribution<%1%>&, %1%)";
0397       if(false == detail::check_pareto(function, scale, shape, &result, Policy()))
0398       {
0399         return result;
0400       }
0401       if (shape > 2)
0402       {
0403         result = (scale * scale * shape) /
0404          ((shape - 1) *  (shape - 1) * (shape - 2));
0405       }
0406       else
0407       {
0408         result = policies::raise_domain_error<RealType>(
0409           function,
0410           "variance is undefined for shape <= 2, but got %1%.", dist.shape(), Policy());
0411       }
0412       return result;
0413     } // variance
0414 
0415     template <class RealType, class Policy>
0416     BOOST_MATH_GPU_ENABLED inline RealType skewness(const pareto_distribution<RealType, Policy>& dist)
0417     {
0418       BOOST_MATH_STD_USING
0419       RealType result = 0;
0420       RealType shape = dist.shape();
0421       constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
0422       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
0423       {
0424         return result;
0425       }
0426       if (shape > 3)
0427       {
0428         result = sqrt((shape - 2) / shape) *
0429           2 * (shape + 1) /
0430           (shape - 3);
0431       }
0432       else
0433       {
0434         result = policies::raise_domain_error<RealType>(
0435           function,
0436           "skewness is undefined for shape <= 3, but got %1%.", dist.shape(), Policy());
0437       }
0438       return result;
0439     } // skewness
0440 
0441     template <class RealType, class Policy>
0442     BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const pareto_distribution<RealType, Policy>& dist)
0443     {
0444       RealType result = 0;
0445       RealType shape = dist.shape();
0446       constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
0447       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
0448       {
0449         return result;
0450       }
0451       if (shape > 4)
0452       {
0453         result = 3 * ((shape - 2) * (3 * shape * shape + shape + 2)) /
0454           (shape * (shape - 3) * (shape - 4));
0455       }
0456       else
0457       {
0458         result = policies::raise_domain_error<RealType>(
0459           function,
0460           "kurtosis_excess is undefined for shape <= 4, but got %1%.", shape, Policy());
0461       }
0462       return result;
0463     } // kurtosis
0464 
0465     template <class RealType, class Policy>
0466     BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const pareto_distribution<RealType, Policy>& dist)
0467     {
0468       RealType result = 0;
0469       RealType shape = dist.shape();
0470       constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
0471       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
0472       {
0473         return result;
0474       }
0475       if (shape > 4)
0476       {
0477         result = 6 * ((shape * shape * shape) + (shape * shape) - 6 * shape - 2) /
0478           (shape * (shape - 3) * (shape - 4));
0479       }
0480       else
0481       {
0482         result = policies::raise_domain_error<RealType>(
0483           function,
0484           "kurtosis_excess is undefined for shape <= 4, but got %1%.", dist.shape(), Policy());
0485       }
0486       return result;
0487     } // kurtosis_excess
0488 
0489     template <class RealType, class Policy>
0490     BOOST_MATH_GPU_ENABLED inline RealType entropy(const pareto_distribution<RealType, Policy>& dist)
0491     {
0492       BOOST_MATH_STD_USING
0493       RealType xm = dist.scale();
0494       RealType alpha = dist.shape();
0495       return log(xm/alpha) + 1 + 1/alpha;
0496     }
0497 
0498     } // namespace math
0499   } // namespace boost
0500 
0501   // This include must be at the end, *after* the accessors
0502   // for this distribution have been defined, in order to
0503   // keep compilers that support two-phase lookup happy.
0504 #include <boost/math/distributions/detail/derived_accessors.hpp>
0505 
0506 #endif // BOOST_STATS_PARETO_HPP
0507 
0508
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!