ceph/src/boost/boost/math/distributions/pareto.hpp

   1 //  Copyright John Maddock 2007.
   2 //  Copyright Paul A. Bristow 2007, 2009
   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
   7 #ifndef BOOST_STATS_PARETO_HPP
   8 #define BOOST_STATS_PARETO_HPP
   9
  10 // http://en.wikipedia.org/wiki/Pareto_distribution
  11 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
  12 // Also:
  13 // Weisstein, Eric W. "Pareto Distribution."
  14 // From MathWorld--A Wolfram Web Resource.
  15 // http://mathworld.wolfram.com/ParetoDistribution.html
  16 // Handbook of Statistical Distributions with Applications, K Krishnamoorthy, ISBN 1-58488-635-8, Chapter 23, pp 257 - 267.
  17 // Caution KK's a and b are the reverse of Mathworld!
  18
  19 #include <boost/math/distributions/fwd.hpp>
  20 #include <boost/math/distributions/complement.hpp>
  21 #include <boost/math/distributions/detail/common_error_handling.hpp>
  22 #include <boost/math/special_functions/powm1.hpp>
  23
  24 #include <utility> // for BOOST_CURRENT_VALUE?
  25
  26 namespace boost
  27 {
  28   namespace math
  29   {
  30     namespace detail
  31     { // Parameter checking.
  32       template <class RealType, class Policy>
  33       inline bool check_pareto_scale(
  34         const char* function,
  35         RealType scale,
  36         RealType* result, const Policy& pol)
  37       {
  38         if((boost::math::isfinite)(scale))
  39         { // any > 0 finite value is OK.
  40           if (scale > 0)
  41           {
  42             return true;
  43           }
  44           else
  45           {
  46             *result = policies::raise_domain_error<RealType>(
  47               function,
  48               "Scale parameter is %1%, but must be > 0!", scale, pol);
  49             return false;
  50           }
  51         }
  52         else
  53         { // Not finite.
  54           *result = policies::raise_domain_error<RealType>(
  55             function,
  56             "Scale parameter is %1%, but must be finite!", scale, pol);
  57           return false;
  58         }
  59       } // bool check_pareto_scale
  60
  61       template <class RealType, class Policy>
  62       inline bool check_pareto_shape(
  63         const char* function,
  64         RealType shape,
  65         RealType* result, const Policy& pol)
  66       {
  67         if((boost::math::isfinite)(shape))
  68         { // Any finite value > 0 is OK.
  69           if (shape > 0)
  70           {
  71             return true;
  72           }
  73           else
  74           {
  75             *result = policies::raise_domain_error<RealType>(
  76               function,
  77               "Shape parameter is %1%, but must be > 0!", shape, pol);
  78             return false;
  79           }
  80         }
  81         else
  82         { // Not finite.
  83           *result = policies::raise_domain_error<RealType>(
  84             function,
  85             "Shape parameter is %1%, but must be finite!", shape, pol);
  86           return false;
  87         }
  88       } // bool check_pareto_shape(
  89
  90       template <class RealType, class Policy>
  91       inline bool check_pareto_x(
  92         const char* function,
  93         RealType const& x,
  94         RealType* result, const Policy& pol)
  95       {
  96         if((boost::math::isfinite)(x))
  97         { //
  98           if (x > 0)
  99           {
 100             return true;
 101           }
 102           else
 103           {
 104             *result = policies::raise_domain_error<RealType>(
 105               function,
 106               "x parameter is %1%, but must be > 0 !", x, pol);
 107             return false;
 108           }
 109         }
 110         else
 111         { // Not finite..
 112           *result = policies::raise_domain_error<RealType>(
 113             function,
 114             "x parameter is %1%, but must be finite!", x, pol);
 115           return false;
 116         }
 117       } // bool check_pareto_x
 118
 119       template <class RealType, class Policy>
 120       inline bool check_pareto( // distribution parameters.
 121         const char* function,
 122         RealType scale,
 123         RealType shape,
 124         RealType* result, const Policy& pol)
 125       {
 126         return check_pareto_scale(function, scale, result, pol)
 127            && check_pareto_shape(function, shape, result, pol);
 128       } // bool check_pareto(
 129
 130     } // namespace detail
 131
 132     template <class RealType = double, class Policy = policies::policy<> >
 133     class pareto_distribution
 134     {
 135     public:
 136       typedef RealType value_type;
 137       typedef Policy policy_type;
 138
 139       pareto_distribution(RealType l_scale = 1, RealType l_shape = 1)
 140         : m_scale(l_scale), m_shape(l_shape)
 141       { // Constructor.
 142         RealType result = 0;
 143         detail::check_pareto("boost::math::pareto_distribution<%1%>::pareto_distribution", l_scale, l_shape, &result, Policy());
 144       }
 145
 146       RealType scale()const
 147       { // AKA Xm and Wolfram b and beta
 148         return m_scale;
 149       }
 150
 151       RealType shape()const
 152       { // AKA k and Wolfram a and alpha
 153         return m_shape;
 154       }
 155     private:
 156       // Data members:
 157       RealType m_scale;  // distribution scale (xm) or beta
 158       RealType m_shape;  // distribution shape (k) or alpha
 159     };
 160
 161     typedef pareto_distribution<double> pareto; // Convenience to allow pareto(2., 3.);
 162
 163     #ifdef __cpp_deduction_guides
 164     template <class RealType>
 165     pareto_distribution(RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;
 166     template <class RealType>
 167     pareto_distribution(RealType,RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;
 168     #endif
 169
 170
 171     template <class RealType, class Policy>
 172     inline const std::pair<RealType, RealType> range(const pareto_distribution<RealType, Policy>& /*dist*/)
 173     { // Range of permissible values for random variable x.
 174       using boost::math::tools::max_value;
 175       return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // scale zero to + infinity.
 176     } // range
 177
 178     template <class RealType, class Policy>
 179     inline const std::pair<RealType, RealType> support(const pareto_distribution<RealType, Policy>& dist)
 180     { // Range of supported values for random variable x.
 181       // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 182       using boost::math::tools::max_value;
 183       return std::pair<RealType, RealType>(dist.scale(), max_value<RealType>() ); // scale to + infinity.
 184     } // support
 185
 186     template <class RealType, class Policy>
 187     inline RealType pdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
 188     {
 189       BOOST_MATH_STD_USING  // for ADL of std function pow.
 190       static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
 191       RealType scale = dist.scale();
 192       RealType shape = dist.shape();
 193       RealType result = 0;
 194       if(false == (detail::check_pareto_x(function, x, &result, Policy())
 195          && detail::check_pareto(function, scale, shape, &result, Policy())))
 196          return result;
 197       if (x < scale)
 198       { // regardless of shape, pdf is zero (or should be disallow x < scale and throw an exception?).
 199         return 0;
 200       }
 201       result = shape * pow(scale, shape) / pow(x, shape+1);
 202       return result;
 203     } // pdf
 204
 205     template <class RealType, class Policy>
 206     inline RealType cdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
 207     {
 208       BOOST_MATH_STD_USING  // for ADL of std function pow.
 209       static const char* function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";
 210       RealType scale = dist.scale();
 211       RealType shape = dist.shape();
 212       RealType result = 0;
 213
 214       if(false == (detail::check_pareto_x(function, x, &result, Policy())
 215          && detail::check_pareto(function, scale, shape, &result, Policy())))
 216          return result;
 217
 218       if (x <= scale)
 219       { // regardless of shape, cdf is zero.
 220         return 0;
 221       }
 222
 223       // result = RealType(1) - pow((scale / x), shape);
 224       result = -boost::math::powm1(scale/x, shape, Policy()); // should be more accurate.
 225       return result;
 226     } // cdf
 227
 228     template <class RealType, class Policy>
 229     inline RealType quantile(const pareto_distribution<RealType, Policy>& dist, const RealType& p)
 230     {
 231       BOOST_MATH_STD_USING  // for ADL of std function pow.
 232       static const char* function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";
 233       RealType result = 0;
 234       RealType scale = dist.scale();
 235       RealType shape = dist.shape();
 236       if(false == (detail::check_probability(function, p, &result, Policy())
 237            && detail::check_pareto(function, scale, shape, &result, Policy())))
 238       {
 239         return result;
 240       }
 241       if (p == 0)
 242       {
 243         return scale; // x must be scale (or less).
 244       }
 245       if (p == 1)
 246       {
 247         return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.
 248       }
 249       result = scale /
 250         (pow((1 - p), 1 / shape));
 251       // K. Krishnamoorthy,  ISBN 1-58488-635-8 eq 23.1.3
 252       return result;
 253     } // quantile
 254
 255     template <class RealType, class Policy>
 256     inline RealType cdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
 257     {
 258        BOOST_MATH_STD_USING  // for ADL of std function pow.
 259        static const char* function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";
 260        RealType result = 0;
 261        RealType x = c.param;
 262        RealType scale = c.dist.scale();
 263        RealType shape = c.dist.shape();
 264        if(false == (detail::check_pareto_x(function, x, &result, Policy())
 265            && detail::check_pareto(function, scale, shape, &result, Policy())))
 266          return result;
 267
 268        if (x <= scale)
 269        { // regardless of shape, cdf is zero, and complement is unity.
 270          return 1;
 271        }
 272        result = pow((scale/x), shape);
 273
 274        return result;
 275     } // cdf complement
 276
 277     template <class RealType, class Policy>
 278     inline RealType quantile(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
 279     {
 280       BOOST_MATH_STD_USING  // for ADL of std function pow.
 281       static const char* function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";
 282       RealType result = 0;
 283       RealType q = c.param;
 284       RealType scale = c.dist.scale();
 285       RealType shape = c.dist.shape();
 286       if(false == (detail::check_probability(function, q, &result, Policy())
 287            && detail::check_pareto(function, scale, shape, &result, Policy())))
 288       {
 289         return result;
 290       }
 291       if (q == 1)
 292       {
 293         return scale; // x must be scale (or less).
 294       }
 295       if (q == 0)
 296       {
 297          return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.
 298       }
 299       result = scale / (pow(q, 1 / shape));
 300       // K. Krishnamoorthy,  ISBN 1-58488-635-8 eq 23.1.3
 301       return result;
 302     } // quantile complement
 303
 304     template <class RealType, class Policy>
 305     inline RealType mean(const pareto_distribution<RealType, Policy>& dist)
 306     {
 307       RealType result = 0;
 308       static const char* function = "boost::math::mean(const pareto_distribution<%1%>&, %1%)";
 309       if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))
 310       {
 311         return result;
 312       }
 313       if (dist.shape() > RealType(1))
 314       {
 315         return dist.shape() * dist.scale() / (dist.shape() - 1);
 316       }
 317       else
 318       {
 319         using boost::math::tools::max_value;
 320         return max_value<RealType>(); // +infinity.
 321       }
 322     } // mean
 323
 324     template <class RealType, class Policy>
 325     inline RealType mode(const pareto_distribution<RealType, Policy>& dist)
 326     {
 327       return dist.scale();
 328     } // mode
 329
 330     template <class RealType, class Policy>
 331     inline RealType median(const pareto_distribution<RealType, Policy>& dist)
 332     {
 333       RealType result = 0;
 334       static const char* function = "boost::math::median(const pareto_distribution<%1%>&, %1%)";
 335       if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))
 336       {
 337         return result;
 338       }
 339       BOOST_MATH_STD_USING
 340       return dist.scale() * pow(RealType(2), (1/dist.shape()));
 341     } // median
 342
 343     template <class RealType, class Policy>
 344     inline RealType variance(const pareto_distribution<RealType, Policy>& dist)
 345     {
 346       RealType result = 0;
 347       RealType scale = dist.scale();
 348       RealType shape = dist.shape();
 349       static const char* function = "boost::math::variance(const pareto_distribution<%1%>&, %1%)";
 350       if(false == detail::check_pareto(function, scale, shape, &result, Policy()))
 351       {
 352         return result;
 353       }
 354       if (shape > 2)
 355       {
 356         result = (scale * scale * shape) /
 357          ((shape - 1) *  (shape - 1) * (shape - 2));
 358       }
 359       else
 360       {
 361         result = policies::raise_domain_error<RealType>(
 362           function,
 363           "variance is undefined for shape <= 2, but got %1%.", dist.shape(), Policy());
 364       }
 365       return result;
 366     } // variance
 367
 368     template <class RealType, class Policy>
 369     inline RealType skewness(const pareto_distribution<RealType, Policy>& dist)
 370     {
 371       BOOST_MATH_STD_USING
 372       RealType result = 0;
 373       RealType shape = dist.shape();
 374       static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
 375       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
 376       {
 377         return result;
 378       }
 379       if (shape > 3)
 380       {
 381         result = sqrt((shape - 2) / shape) *
 382           2 * (shape + 1) /
 383           (shape - 3);
 384       }
 385       else
 386       {
 387         result = policies::raise_domain_error<RealType>(
 388           function,
 389           "skewness is undefined for shape <= 3, but got %1%.", dist.shape(), Policy());
 390       }
 391       return result;
 392     } // skewness
 393
 394     template <class RealType, class Policy>
 395     inline RealType kurtosis(const pareto_distribution<RealType, Policy>& dist)
 396     {
 397       RealType result = 0;
 398       RealType shape = dist.shape();
 399       static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
 400       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
 401       {
 402         return result;
 403       }
 404       if (shape > 4)
 405       {
 406         result = 3 * ((shape - 2) * (3 * shape * shape + shape + 2)) /
 407           (shape * (shape - 3) * (shape - 4));
 408       }
 409       else
 410       {
 411         result = policies::raise_domain_error<RealType>(
 412           function,
 413           "kurtosis_excess is undefined for shape <= 4, but got %1%.", shape, Policy());
 414       }
 415       return result;
 416     } // kurtosis
 417
 418     template <class RealType, class Policy>
 419     inline RealType kurtosis_excess(const pareto_distribution<RealType, Policy>& dist)
 420     {
 421       RealType result = 0;
 422       RealType shape = dist.shape();
 423       static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";
 424       if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))
 425       {
 426         return result;
 427       }
 428       if (shape > 4)
 429       {
 430         result = 6 * ((shape * shape * shape) + (shape * shape) - 6 * shape - 2) /
 431           (shape * (shape - 3) * (shape - 4));
 432       }
 433       else
 434       {
 435         result = policies::raise_domain_error<RealType>(
 436           function,
 437           "kurtosis_excess is undefined for shape <= 4, but got %1%.", dist.shape(), Policy());
 438       }
 439       return result;
 440     } // kurtosis_excess
 441
 442     template <class RealType, class Policy>
 443     inline RealType entropy(const pareto_distribution<RealType, Policy>& dist)
 444     {
 445       using std::log;
 446       RealType xm = dist.scale();
 447       RealType alpha = dist.shape();
 448       return log(xm/alpha) + 1 + 1/alpha;
 449     }
 450
 451     } // namespace math
 452   } // namespace boost
 453
 454   // This include must be at the end, *after* the accessors
 455   // for this distribution have been defined, in order to
 456   // keep compilers that support two-phase lookup happy.
 457 #include <boost/math/distributions/detail/derived_accessors.hpp>
 458
 459 #endif // BOOST_STATS_PARETO_HPP
 460
 461