pareto.hpp 15 KB

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