math/distributions/pareto.hpp

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