gamma.hpp 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359
  1. // Copyright John Maddock 2006.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_STATS_GAMMA_HPP
  6. #define BOOST_STATS_GAMMA_HPP
  7. // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
  8. // http://mathworld.wolfram.com/GammaDistribution.html
  9. // http://en.wikipedia.org/wiki/Gamma_distribution
  10. #include <boost/math/distributions/fwd.hpp>
  11. #include <boost/math/special_functions/gamma.hpp>
  12. #include <boost/math/special_functions/digamma.hpp>
  13. #include <boost/math/distributions/detail/common_error_handling.hpp>
  14. #include <boost/math/distributions/complement.hpp>
  15. #include <utility>
  16. namespace boost{ namespace math
  17. {
  18. namespace detail
  19. {
  20. template <class RealType, class Policy>
  21. inline bool check_gamma_shape(
  22. const char* function,
  23. RealType shape,
  24. RealType* result, const Policy& pol)
  25. {
  26. if((shape <= 0) || !(boost::math::isfinite)(shape))
  27. {
  28. *result = policies::raise_domain_error<RealType>(
  29. function,
  30. "Shape parameter is %1%, but must be > 0 !", shape, pol);
  31. return false;
  32. }
  33. return true;
  34. }
  35. template <class RealType, class Policy>
  36. inline bool check_gamma_x(
  37. const char* function,
  38. RealType const& x,
  39. RealType* result, const Policy& pol)
  40. {
  41. if((x < 0) || !(boost::math::isfinite)(x))
  42. {
  43. *result = policies::raise_domain_error<RealType>(
  44. function,
  45. "Random variate is %1% but must be >= 0 !", x, pol);
  46. return false;
  47. }
  48. return true;
  49. }
  50. template <class RealType, class Policy>
  51. inline bool check_gamma(
  52. const char* function,
  53. RealType scale,
  54. RealType shape,
  55. RealType* result, const Policy& pol)
  56. {
  57. return check_scale(function, scale, result, pol) && check_gamma_shape(function, shape, result, pol);
  58. }
  59. } // namespace detail
  60. template <class RealType = double, class Policy = policies::policy<> >
  61. class gamma_distribution
  62. {
  63. public:
  64. typedef RealType value_type;
  65. typedef Policy policy_type;
  66. gamma_distribution(RealType l_shape, RealType l_scale = 1)
  67. : m_shape(l_shape), m_scale(l_scale)
  68. {
  69. RealType result;
  70. detail::check_gamma("boost::math::gamma_distribution<%1%>::gamma_distribution", l_scale, l_shape, &result, Policy());
  71. }
  72. RealType shape()const
  73. {
  74. return m_shape;
  75. }
  76. RealType scale()const
  77. {
  78. return m_scale;
  79. }
  80. private:
  81. //
  82. // Data members:
  83. //
  84. RealType m_shape; // distribution shape
  85. RealType m_scale; // distribution scale
  86. };
  87. // NO typedef because of clash with name of gamma function.
  88. template <class RealType, class Policy>
  89. inline const std::pair<RealType, RealType> range(const gamma_distribution<RealType, Policy>& /* dist */)
  90. { // Range of permissible values for random variable x.
  91. using boost::math::tools::max_value;
  92. return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  93. }
  94. template <class RealType, class Policy>
  95. inline const std::pair<RealType, RealType> support(const gamma_distribution<RealType, Policy>& /* dist */)
  96. { // Range of supported values for random variable x.
  97. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  98. using boost::math::tools::max_value;
  99. using boost::math::tools::min_value;
  100. return std::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>());
  101. }
  102. template <class RealType, class Policy>
  103. inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
  104. {
  105. BOOST_MATH_STD_USING // for ADL of std functions
  106. static const char* function = "boost::math::pdf(const gamma_distribution<%1%>&, %1%)";
  107. RealType shape = dist.shape();
  108. RealType scale = dist.scale();
  109. RealType result = 0;
  110. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  111. return result;
  112. if(false == detail::check_gamma_x(function, x, &result, Policy()))
  113. return result;
  114. if(x == 0)
  115. {
  116. return 0;
  117. }
  118. result = gamma_p_derivative(shape, x / scale, Policy()) / scale;
  119. return result;
  120. } // pdf
  121. template <class RealType, class Policy>
  122. inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
  123. {
  124. BOOST_MATH_STD_USING // for ADL of std functions
  125. static const char* function = "boost::math::cdf(const gamma_distribution<%1%>&, %1%)";
  126. RealType shape = dist.shape();
  127. RealType scale = dist.scale();
  128. RealType result = 0;
  129. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  130. return result;
  131. if(false == detail::check_gamma_x(function, x, &result, Policy()))
  132. return result;
  133. result = boost::math::gamma_p(shape, x / scale, Policy());
  134. return result;
  135. } // cdf
  136. template <class RealType, class Policy>
  137. inline RealType quantile(const gamma_distribution<RealType, Policy>& dist, const RealType& p)
  138. {
  139. BOOST_MATH_STD_USING // for ADL of std functions
  140. static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
  141. RealType shape = dist.shape();
  142. RealType scale = dist.scale();
  143. RealType result = 0;
  144. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  145. return result;
  146. if(false == detail::check_probability(function, p, &result, Policy()))
  147. return result;
  148. if(p == 1)
  149. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  150. result = gamma_p_inv(shape, p, Policy()) * scale;
  151. return result;
  152. }
  153. template <class RealType, class Policy>
  154. inline RealType cdf(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
  155. {
  156. BOOST_MATH_STD_USING // for ADL of std functions
  157. static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
  158. RealType shape = c.dist.shape();
  159. RealType scale = c.dist.scale();
  160. RealType result = 0;
  161. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  162. return result;
  163. if(false == detail::check_gamma_x(function, c.param, &result, Policy()))
  164. return result;
  165. result = gamma_q(shape, c.param / scale, Policy());
  166. return result;
  167. }
  168. template <class RealType, class Policy>
  169. inline RealType quantile(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
  170. {
  171. BOOST_MATH_STD_USING // for ADL of std functions
  172. static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
  173. RealType shape = c.dist.shape();
  174. RealType scale = c.dist.scale();
  175. RealType q = c.param;
  176. RealType result = 0;
  177. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  178. return result;
  179. if(false == detail::check_probability(function, q, &result, Policy()))
  180. return result;
  181. if(q == 0)
  182. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  183. result = gamma_q_inv(shape, q, Policy()) * scale;
  184. return result;
  185. }
  186. template <class RealType, class Policy>
  187. inline RealType mean(const gamma_distribution<RealType, Policy>& dist)
  188. {
  189. BOOST_MATH_STD_USING // for ADL of std functions
  190. static const char* function = "boost::math::mean(const gamma_distribution<%1%>&)";
  191. RealType shape = dist.shape();
  192. RealType scale = dist.scale();
  193. RealType result = 0;
  194. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  195. return result;
  196. result = shape * scale;
  197. return result;
  198. }
  199. template <class RealType, class Policy>
  200. inline RealType variance(const gamma_distribution<RealType, Policy>& dist)
  201. {
  202. BOOST_MATH_STD_USING // for ADL of std functions
  203. static const char* function = "boost::math::variance(const gamma_distribution<%1%>&)";
  204. RealType shape = dist.shape();
  205. RealType scale = dist.scale();
  206. RealType result = 0;
  207. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  208. return result;
  209. result = shape * scale * scale;
  210. return result;
  211. }
  212. template <class RealType, class Policy>
  213. inline RealType mode(const gamma_distribution<RealType, Policy>& dist)
  214. {
  215. BOOST_MATH_STD_USING // for ADL of std functions
  216. static const char* function = "boost::math::mode(const gamma_distribution<%1%>&)";
  217. RealType shape = dist.shape();
  218. RealType scale = dist.scale();
  219. RealType result = 0;
  220. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  221. return result;
  222. if(shape < 1)
  223. return policies::raise_domain_error<RealType>(
  224. function,
  225. "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.",
  226. shape, Policy());
  227. result = (shape - 1) * scale;
  228. return result;
  229. }
  230. //template <class RealType, class Policy>
  231. //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
  232. //{ // Rely on default definition in derived accessors.
  233. //}
  234. template <class RealType, class Policy>
  235. inline RealType skewness(const gamma_distribution<RealType, Policy>& dist)
  236. {
  237. BOOST_MATH_STD_USING // for ADL of std functions
  238. static const char* function = "boost::math::skewness(const gamma_distribution<%1%>&)";
  239. RealType shape = dist.shape();
  240. RealType scale = dist.scale();
  241. RealType result = 0;
  242. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  243. return result;
  244. result = 2 / sqrt(shape);
  245. return result;
  246. }
  247. template <class RealType, class Policy>
  248. inline RealType kurtosis_excess(const gamma_distribution<RealType, Policy>& dist)
  249. {
  250. BOOST_MATH_STD_USING // for ADL of std functions
  251. static const char* function = "boost::math::kurtosis_excess(const gamma_distribution<%1%>&)";
  252. RealType shape = dist.shape();
  253. RealType scale = dist.scale();
  254. RealType result = 0;
  255. if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
  256. return result;
  257. result = 6 / shape;
  258. return result;
  259. }
  260. template <class RealType, class Policy>
  261. inline RealType kurtosis(const gamma_distribution<RealType, Policy>& dist)
  262. {
  263. return kurtosis_excess(dist) + 3;
  264. }
  265. template <class RealType, class Policy>
  266. inline RealType entropy(const gamma_distribution<RealType, Policy>& dist)
  267. {
  268. RealType k = dist.shape();
  269. RealType theta = dist.scale();
  270. using std::log;
  271. return k + log(theta) + lgamma(k) + (1-k)*digamma(k);
  272. }
  273. } // namespace math
  274. } // namespace boost
  275. // This include must be at the end, *after* the accessors
  276. // for this distribution have been defined, in order to
  277. // keep compilers that support two-phase lookup happy.
  278. #include <boost/math/distributions/detail/derived_accessors.hpp>
  279. #endif // BOOST_STATS_GAMMA_HPP