509 lines · plain
1// Copyright John Maddock 2007.2// Copyright Paul A. Bristow 2007, 20093// Copyright Matt Borland 2023.4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_STATS_PARETO_HPP9#define BOOST_STATS_PARETO_HPP10 11// http://en.wikipedia.org/wiki/Pareto_distribution12// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm13// Also:14// Weisstein, Eric W. "Pareto Distribution."15// From MathWorld--A Wolfram Web Resource.16// http://mathworld.wolfram.com/ParetoDistribution.html17// Handbook of Statistical Distributions with Applications, K Krishnamoorthy, ISBN 1-58488-635-8, Chapter 23, pp 257 - 267.18// Caution KK's a and b are the reverse of Mathworld!19 20#include <boost/math/tools/config.hpp>21#include <boost/math/tools/tuple.hpp>22#include <boost/math/tools/numeric_limits.hpp>23#include <boost/math/distributions/fwd.hpp>24#include <boost/math/distributions/complement.hpp>25#include <boost/math/distributions/detail/common_error_handling.hpp>26#include <boost/math/special_functions/powm1.hpp>27#include <boost/math/special_functions/log1p.hpp>28 29namespace boost30{31 namespace math32 {33 namespace detail34 { // Parameter checking.35 template <class RealType, class Policy>36 BOOST_MATH_GPU_ENABLED inline bool check_pareto_scale(37 const char* function,38 RealType scale,39 RealType* result, const Policy& pol)40 {41 if((boost::math::isfinite)(scale))42 { // any > 0 finite value is OK.43 if (scale > 0)44 {45 return true;46 }47 else48 {49 *result = policies::raise_domain_error<RealType>(50 function,51 "Scale parameter is %1%, but must be > 0!", scale, pol);52 return false;53 }54 }55 else56 { // Not finite.57 *result = policies::raise_domain_error<RealType>(58 function,59 "Scale parameter is %1%, but must be finite!", scale, pol);60 return false;61 }62 } // bool check_pareto_scale63 64 template <class RealType, class Policy>65 BOOST_MATH_GPU_ENABLED inline bool check_pareto_shape(66 const char* function,67 RealType shape,68 RealType* result, const Policy& pol)69 {70 if((boost::math::isfinite)(shape))71 { // Any finite value > 0 is OK.72 if (shape > 0)73 {74 return true;75 }76 else77 {78 *result = policies::raise_domain_error<RealType>(79 function,80 "Shape parameter is %1%, but must be > 0!", shape, pol);81 return false;82 }83 }84 else85 { // Not finite.86 *result = policies::raise_domain_error<RealType>(87 function,88 "Shape parameter is %1%, but must be finite!", shape, pol);89 return false;90 }91 } // bool check_pareto_shape(92 93 template <class RealType, class Policy>94 BOOST_MATH_GPU_ENABLED inline bool check_pareto_x(95 const char* function,96 RealType const& x,97 RealType* result, const Policy& pol)98 {99 if((boost::math::isfinite)(x))100 { //101 if (x > 0)102 {103 return true;104 }105 else106 {107 *result = policies::raise_domain_error<RealType>(108 function,109 "x parameter is %1%, but must be > 0 !", x, pol);110 return false;111 }112 }113 else114 { // Not finite..115 *result = policies::raise_domain_error<RealType>(116 function,117 "x parameter is %1%, but must be finite!", x, pol);118 return false;119 }120 } // bool check_pareto_x121 122 template <class RealType, class Policy>123 BOOST_MATH_GPU_ENABLED inline bool check_pareto( // distribution parameters.124 const char* function,125 RealType scale,126 RealType shape,127 RealType* result, const Policy& pol)128 {129 return check_pareto_scale(function, scale, result, pol)130 && check_pareto_shape(function, shape, result, pol);131 } // bool check_pareto(132 133 } // namespace detail134 135 template <class RealType = double, class Policy = policies::policy<> >136 class pareto_distribution137 {138 public:139 typedef RealType value_type;140 typedef Policy policy_type;141 142 BOOST_MATH_GPU_ENABLED pareto_distribution(RealType l_scale = 1, RealType l_shape = 1)143 : m_scale(l_scale), m_shape(l_shape)144 { // Constructor.145 RealType result = 0;146 detail::check_pareto("boost::math::pareto_distribution<%1%>::pareto_distribution", l_scale, l_shape, &result, Policy());147 }148 149 BOOST_MATH_GPU_ENABLED RealType scale()const150 { // AKA Xm and Wolfram b and beta151 return m_scale;152 }153 154 BOOST_MATH_GPU_ENABLED RealType shape()const155 { // AKA k and Wolfram a and alpha156 return m_shape;157 }158 private:159 // Data members:160 RealType m_scale; // distribution scale (xm) or beta161 RealType m_shape; // distribution shape (k) or alpha162 };163 164 typedef pareto_distribution<double> pareto; // Convenience to allow pareto(2., 3.);165 166 #ifdef __cpp_deduction_guides167 template <class RealType>168 pareto_distribution(RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;169 template <class RealType>170 pareto_distribution(RealType,RealType)->pareto_distribution<typename boost::math::tools::promote_args<RealType>::type>;171 #endif172 173 174 template <class RealType, class Policy>175 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const pareto_distribution<RealType, Policy>& /*dist*/)176 { // Range of permissible values for random variable x.177 using boost::math::tools::max_value;178 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // scale zero to + infinity.179 } // range180 181 template <class RealType, class Policy>182 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const pareto_distribution<RealType, Policy>& dist)183 { // Range of supported values for random variable x.184 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.185 using boost::math::tools::max_value;186 return boost::math::pair<RealType, RealType>(dist.scale(), max_value<RealType>() ); // scale to + infinity.187 } // support188 189 template <class RealType, class Policy>190 BOOST_MATH_GPU_ENABLED inline RealType pdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)191 {192 BOOST_MATH_STD_USING // for ADL of std function pow.193 constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";194 RealType scale = dist.scale();195 RealType shape = dist.shape();196 RealType result = 0;197 if(false == (detail::check_pareto_x(function, x, &result, Policy())198 && detail::check_pareto(function, scale, shape, &result, Policy())))199 return result;200 if (x < scale)201 { // regardless of shape, pdf is zero (or should be disallow x < scale and throw an exception?).202 return 0;203 }204 result = shape * pow(scale, shape) / pow(x, shape+1);205 return result;206 } // pdf207 208 template <class RealType, class Policy>209 BOOST_MATH_GPU_ENABLED inline RealType cdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)210 {211 BOOST_MATH_STD_USING // for ADL of std function pow.212 constexpr auto function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";213 RealType scale = dist.scale();214 RealType shape = dist.shape();215 RealType result = 0;216 217 if(false == (detail::check_pareto_x(function, x, &result, Policy())218 && detail::check_pareto(function, scale, shape, &result, Policy())))219 return result;220 221 if (x <= scale)222 { // regardless of shape, cdf is zero.223 return 0;224 }225 226 // result = RealType(1) - pow((scale / x), shape);227 result = -boost::math::powm1(scale/x, shape, Policy()); // should be more accurate.228 return result;229 } // cdf230 231 template <class RealType, class Policy>232 BOOST_MATH_GPU_ENABLED inline RealType logcdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)233 {234 BOOST_MATH_STD_USING // for ADL of std function pow.235 constexpr auto function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";236 RealType scale = dist.scale();237 RealType shape = dist.shape();238 RealType result = 0;239 240 if(false == (detail::check_pareto_x(function, x, &result, Policy())241 && detail::check_pareto(function, scale, shape, &result, Policy())))242 return result;243 244 if (x <= scale)245 { // regardless of shape, cdf is zero.246 return -boost::math::numeric_limits<RealType>::infinity();247 }248 249 result = log1p(-pow(scale/x, shape), Policy());250 return result;251 } // logcdf252 253 template <class RealType, class Policy>254 BOOST_MATH_GPU_ENABLED inline RealType quantile(const pareto_distribution<RealType, Policy>& dist, const RealType& p)255 {256 BOOST_MATH_STD_USING // for ADL of std function pow.257 constexpr auto function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";258 RealType result = 0;259 RealType scale = dist.scale();260 RealType shape = dist.shape();261 if(false == (detail::check_probability(function, p, &result, Policy())262 && detail::check_pareto(function, scale, shape, &result, Policy())))263 {264 return result;265 }266 if (p == 0)267 {268 return scale; // x must be scale (or less).269 }270 if (p == 1)271 {272 return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.273 }274 result = scale /275 (pow((1 - p), 1 / shape));276 // K. Krishnamoorthy, ISBN 1-58488-635-8 eq 23.1.3277 return result;278 } // quantile279 280 template <class RealType, class Policy>281 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)282 {283 BOOST_MATH_STD_USING // for ADL of std function pow.284 constexpr auto function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)";285 RealType result = 0;286 RealType x = c.param;287 RealType scale = c.dist.scale();288 RealType shape = c.dist.shape();289 if(false == (detail::check_pareto_x(function, x, &result, Policy())290 && detail::check_pareto(function, scale, shape, &result, Policy())))291 return result;292 293 if (x <= scale)294 { // regardless of shape, cdf is zero, and complement is unity.295 return 1;296 }297 result = pow((scale/x), shape);298 299 return result;300 } // cdf complement301 302 template <class RealType, class Policy>303 BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)304 {305 BOOST_MATH_STD_USING // for ADL of std function pow.306 constexpr auto function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";307 RealType result = 0;308 RealType x = c.param;309 RealType scale = c.dist.scale();310 RealType shape = c.dist.shape();311 if(false == (detail::check_pareto_x(function, x, &result, Policy())312 && detail::check_pareto(function, scale, shape, &result, Policy())))313 return result;314 315 if (x <= scale)316 { // regardless of shape, cdf is zero, and complement is unity.317 return 0;318 }319 result = log(pow((scale/x), shape));320 321 return result;322 } // logcdf complement323 324 template <class RealType, class Policy>325 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)326 {327 BOOST_MATH_STD_USING // for ADL of std function pow.328 constexpr auto function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)";329 RealType result = 0;330 RealType q = c.param;331 RealType scale = c.dist.scale();332 RealType shape = c.dist.shape();333 if(false == (detail::check_probability(function, q, &result, Policy())334 && detail::check_pareto(function, scale, shape, &result, Policy())))335 {336 return result;337 }338 if (q == 1)339 {340 return scale; // x must be scale (or less).341 }342 if (q == 0)343 {344 return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity.345 }346 result = scale / (pow(q, 1 / shape));347 // K. Krishnamoorthy, ISBN 1-58488-635-8 eq 23.1.3348 return result;349 } // quantile complement350 351 template <class RealType, class Policy>352 BOOST_MATH_GPU_ENABLED inline RealType mean(const pareto_distribution<RealType, Policy>& dist)353 {354 RealType result = 0;355 constexpr auto function = "boost::math::mean(const pareto_distribution<%1%>&, %1%)";356 if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))357 {358 return result;359 }360 if (dist.shape() > RealType(1))361 {362 return dist.shape() * dist.scale() / (dist.shape() - 1);363 }364 else365 {366 using boost::math::tools::max_value;367 return max_value<RealType>(); // +infinity.368 }369 } // mean370 371 template <class RealType, class Policy>372 BOOST_MATH_GPU_ENABLED inline RealType mode(const pareto_distribution<RealType, Policy>& dist)373 {374 return dist.scale();375 } // mode376 377 template <class RealType, class Policy>378 BOOST_MATH_GPU_ENABLED inline RealType median(const pareto_distribution<RealType, Policy>& dist)379 {380 RealType result = 0;381 constexpr auto function = "boost::math::median(const pareto_distribution<%1%>&, %1%)";382 if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy()))383 {384 return result;385 }386 BOOST_MATH_STD_USING387 return dist.scale() * pow(RealType(2), (1/dist.shape()));388 } // median389 390 template <class RealType, class Policy>391 BOOST_MATH_GPU_ENABLED inline RealType variance(const pareto_distribution<RealType, Policy>& dist)392 {393 RealType result = 0;394 RealType scale = dist.scale();395 RealType shape = dist.shape();396 constexpr auto function = "boost::math::variance(const pareto_distribution<%1%>&, %1%)";397 if(false == detail::check_pareto(function, scale, shape, &result, Policy()))398 {399 return result;400 }401 if (shape > 2)402 {403 result = (scale * scale * shape) /404 ((shape - 1) * (shape - 1) * (shape - 2));405 }406 else407 {408 result = policies::raise_domain_error<RealType>(409 function,410 "variance is undefined for shape <= 2, but got %1%.", dist.shape(), Policy());411 }412 return result;413 } // variance414 415 template <class RealType, class Policy>416 BOOST_MATH_GPU_ENABLED inline RealType skewness(const pareto_distribution<RealType, Policy>& dist)417 {418 BOOST_MATH_STD_USING419 RealType result = 0;420 RealType shape = dist.shape();421 constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";422 if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))423 {424 return result;425 }426 if (shape > 3)427 {428 result = sqrt((shape - 2) / shape) *429 2 * (shape + 1) /430 (shape - 3);431 }432 else433 {434 result = policies::raise_domain_error<RealType>(435 function,436 "skewness is undefined for shape <= 3, but got %1%.", dist.shape(), Policy());437 }438 return result;439 } // skewness440 441 template <class RealType, class Policy>442 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const pareto_distribution<RealType, Policy>& dist)443 {444 RealType result = 0;445 RealType shape = dist.shape();446 constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";447 if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))448 {449 return result;450 }451 if (shape > 4)452 {453 result = 3 * ((shape - 2) * (3 * shape * shape + shape + 2)) /454 (shape * (shape - 3) * (shape - 4));455 }456 else457 {458 result = policies::raise_domain_error<RealType>(459 function,460 "kurtosis_excess is undefined for shape <= 4, but got %1%.", shape, Policy());461 }462 return result;463 } // kurtosis464 465 template <class RealType, class Policy>466 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const pareto_distribution<RealType, Policy>& dist)467 {468 RealType result = 0;469 RealType shape = dist.shape();470 constexpr auto function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)";471 if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy()))472 {473 return result;474 }475 if (shape > 4)476 {477 result = 6 * ((shape * shape * shape) + (shape * shape) - 6 * shape - 2) /478 (shape * (shape - 3) * (shape - 4));479 }480 else481 {482 result = policies::raise_domain_error<RealType>(483 function,484 "kurtosis_excess is undefined for shape <= 4, but got %1%.", dist.shape(), Policy());485 }486 return result;487 } // kurtosis_excess488 489 template <class RealType, class Policy>490 BOOST_MATH_GPU_ENABLED inline RealType entropy(const pareto_distribution<RealType, Policy>& dist)491 {492 BOOST_MATH_STD_USING493 RealType xm = dist.scale();494 RealType alpha = dist.shape();495 return log(xm/alpha) + 1 + 1/alpha;496 }497 498 } // namespace math499 } // namespace boost500 501 // This include must be at the end, *after* the accessors502 // for this distribution have been defined, in order to503 // keep compilers that support two-phase lookup happy.504#include <boost/math/distributions/detail/derived_accessors.hpp>505 506#endif // BOOST_STATS_PARETO_HPP507 508 509