inverse_gamma.hpp 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461
  1. // inverse_gamma.hpp
  2. // Copyright Paul A. Bristow 2010.
  3. // Copyright John Maddock 2010.
  4. // Use, modification and distribution are subject to the
  5. // Boost Software License, Version 1.0. (See accompanying file
  6. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  7. #ifndef BOOST_STATS_INVERSE_GAMMA_HPP
  8. #define BOOST_STATS_INVERSE_GAMMA_HPP
  9. // Inverse Gamma Distribution is a two-parameter family
  10. // of continuous probability distributions
  11. // on the positive real line, which is the distribution of
  12. // the reciprocal of a variable distributed according to the gamma distribution.
  13. // http://en.wikipedia.org/wiki/Inverse-gamma_distribution
  14. // http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html
  15. // See also gamma distribution at gamma.hpp:
  16. // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
  17. // http://mathworld.wolfram.com/GammaDistribution.html
  18. // http://en.wikipedia.org/wiki/Gamma_distribution
  19. #include <boost/math/distributions/fwd.hpp>
  20. #include <boost/math/special_functions/gamma.hpp>
  21. #include <boost/math/distributions/detail/common_error_handling.hpp>
  22. #include <boost/math/distributions/complement.hpp>
  23. #include <utility>
  24. namespace boost{ namespace math
  25. {
  26. namespace detail
  27. {
  28. template <class RealType, class Policy>
  29. inline bool check_inverse_gamma_shape(
  30. const char* function, // inverse_gamma
  31. RealType shape, // shape aka alpha
  32. RealType* result, // to update, perhaps with NaN
  33. const Policy& pol)
  34. { // Sources say shape argument must be > 0
  35. // but seems logical to allow shape zero as special case,
  36. // returning pdf and cdf zero (but not < 0).
  37. // (Functions like mean, variance with other limits on shape are checked
  38. // in version including an operator & limit below).
  39. if((shape < 0) || !(boost::math::isfinite)(shape))
  40. {
  41. *result = policies::raise_domain_error<RealType>(
  42. function,
  43. "Shape parameter is %1%, but must be >= 0 !", shape, pol);
  44. return false;
  45. }
  46. return true;
  47. } //bool check_inverse_gamma_shape
  48. template <class RealType, class Policy>
  49. inline bool check_inverse_gamma_x(
  50. const char* function,
  51. RealType const& x,
  52. RealType* result, const Policy& pol)
  53. {
  54. if((x < 0) || !(boost::math::isfinite)(x))
  55. {
  56. *result = policies::raise_domain_error<RealType>(
  57. function,
  58. "Random variate is %1% but must be >= 0 !", x, pol);
  59. return false;
  60. }
  61. return true;
  62. }
  63. template <class RealType, class Policy>
  64. inline bool check_inverse_gamma(
  65. const char* function, // TODO swap these over, so shape is first.
  66. RealType scale, // scale aka beta
  67. RealType shape, // shape aka alpha
  68. RealType* result, const Policy& pol)
  69. {
  70. return check_scale(function, scale, result, pol)
  71. && check_inverse_gamma_shape(function, shape, result, pol);
  72. } // bool check_inverse_gamma
  73. } // namespace detail
  74. template <class RealType = double, class Policy = policies::policy<> >
  75. class inverse_gamma_distribution
  76. {
  77. public:
  78. typedef RealType value_type;
  79. typedef Policy policy_type;
  80. inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1)
  81. : m_shape(l_shape), m_scale(l_scale)
  82. {
  83. RealType result;
  84. detail::check_inverse_gamma(
  85. "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution",
  86. l_scale, l_shape, &result, Policy());
  87. }
  88. RealType shape()const
  89. {
  90. return m_shape;
  91. }
  92. RealType scale()const
  93. {
  94. return m_scale;
  95. }
  96. private:
  97. //
  98. // Data members:
  99. //
  100. RealType m_shape; // distribution shape
  101. RealType m_scale; // distribution scale
  102. };
  103. typedef inverse_gamma_distribution<double> inverse_gamma;
  104. // typedef - but potential clash with name of inverse gamma *function*.
  105. // but there is a typedef for gamma
  106. // typedef boost::math::gamma_distribution<Type, Policy> gamma;
  107. // Allow random variable x to be zero, treated as a special case (unlike some definitions).
  108. template <class RealType, class Policy>
  109. inline const std::pair<RealType, RealType> range(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
  110. { // Range of permissible values for random variable x.
  111. using boost::math::tools::max_value;
  112. return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  113. }
  114. template <class RealType, class Policy>
  115. inline const std::pair<RealType, RealType> support(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
  116. { // Range of supported values for random variable x.
  117. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  118. using boost::math::tools::max_value;
  119. using boost::math::tools::min_value;
  120. return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  121. }
  122. template <class RealType, class Policy>
  123. inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
  124. {
  125. BOOST_MATH_STD_USING // for ADL of std functions
  126. static const char* function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)";
  127. RealType shape = dist.shape();
  128. RealType scale = dist.scale();
  129. RealType result = 0;
  130. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  131. { // distribution parameters bad.
  132. return result;
  133. }
  134. if(x == 0)
  135. { // Treat random variate zero as a special case.
  136. return 0;
  137. }
  138. else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
  139. { // x bad.
  140. return result;
  141. }
  142. result = scale / x;
  143. if(result < tools::min_value<RealType>())
  144. return 0; // random variable is infinite or so close as to make no difference.
  145. result = gamma_p_derivative(shape, result, Policy()) * scale;
  146. if(0 != result)
  147. {
  148. if(x < 0)
  149. {
  150. // x * x may under or overflow, likewise our result,
  151. // so be extra careful about the arithmetic:
  152. RealType lim = tools::max_value<RealType>() * x;
  153. if(lim < result)
  154. return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
  155. result /= x;
  156. if(lim < result)
  157. return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
  158. result /= x;
  159. }
  160. result /= (x * x);
  161. }
  162. // better than naive
  163. // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape);
  164. return result;
  165. } // pdf
  166. template <class RealType, class Policy>
  167. inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
  168. {
  169. BOOST_MATH_STD_USING // for ADL of std functions
  170. static const char* function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)";
  171. RealType shape = dist.shape();
  172. RealType scale = dist.scale();
  173. RealType result = 0;
  174. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  175. { // distribution parameters bad.
  176. return result;
  177. }
  178. if (x == 0)
  179. { // Treat zero as a special case.
  180. return 0;
  181. }
  182. else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
  183. { // x bad
  184. return result;
  185. }
  186. result = boost::math::gamma_q(shape, scale / x, Policy());
  187. // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma
  188. return result;
  189. } // cdf
  190. template <class RealType, class Policy>
  191. inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)
  192. {
  193. BOOST_MATH_STD_USING // for ADL of std functions
  194. using boost::math::gamma_q_inv;
  195. static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
  196. RealType shape = dist.shape();
  197. RealType scale = dist.scale();
  198. RealType result = 0;
  199. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  200. return result;
  201. if(false == detail::check_probability(function, p, &result, Policy()))
  202. return result;
  203. if(p == 1)
  204. {
  205. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  206. }
  207. result = gamma_q_inv(shape, p, Policy());
  208. if((result < 1) && (result * tools::max_value<RealType>() < scale))
  209. return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
  210. result = scale / result;
  211. return result;
  212. }
  213. template <class RealType, class Policy>
  214. inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
  215. {
  216. BOOST_MATH_STD_USING // for ADL of std functions
  217. static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
  218. RealType shape = c.dist.shape();
  219. RealType scale = c.dist.scale();
  220. RealType result = 0;
  221. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  222. return result;
  223. if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy()))
  224. return result;
  225. if(c.param == 0)
  226. return 1; // Avoid division by zero
  227. //result = 1. - gamma_q(shape, c.param / scale, Policy());
  228. result = gamma_p(shape, scale/c.param, Policy());
  229. return result;
  230. }
  231. template <class RealType, class Policy>
  232. inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
  233. {
  234. BOOST_MATH_STD_USING // for ADL of std functions
  235. static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
  236. RealType shape = c.dist.shape();
  237. RealType scale = c.dist.scale();
  238. RealType q = c.param;
  239. RealType result = 0;
  240. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  241. return result;
  242. if(false == detail::check_probability(function, q, &result, Policy()))
  243. return result;
  244. if(q == 0)
  245. {
  246. return policies::raise_overflow_error<RealType>(function, 0, Policy());
  247. }
  248. result = gamma_p_inv(shape, q, Policy());
  249. if((result < 1) && (result * tools::max_value<RealType>() < scale))
  250. return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
  251. result = scale / result;
  252. return result;
  253. }
  254. template <class RealType, class Policy>
  255. inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist)
  256. {
  257. BOOST_MATH_STD_USING // for ADL of std functions
  258. static const char* function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)";
  259. RealType shape = dist.shape();
  260. RealType scale = dist.scale();
  261. RealType result = 0;
  262. if(false == detail::check_scale(function, scale, &result, Policy()))
  263. {
  264. return result;
  265. }
  266. if((shape <= 1) || !(boost::math::isfinite)(shape))
  267. {
  268. result = policies::raise_domain_error<RealType>(
  269. function,
  270. "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy());
  271. return result;
  272. }
  273. result = scale / (shape - 1);
  274. return result;
  275. } // mean
  276. template <class RealType, class Policy>
  277. inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist)
  278. {
  279. BOOST_MATH_STD_USING // for ADL of std functions
  280. static const char* function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)";
  281. RealType shape = dist.shape();
  282. RealType scale = dist.scale();
  283. RealType result = 0;
  284. if(false == detail::check_scale(function, scale, &result, Policy()))
  285. {
  286. return result;
  287. }
  288. if((shape <= 2) || !(boost::math::isfinite)(shape))
  289. {
  290. result = policies::raise_domain_error<RealType>(
  291. function,
  292. "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy());
  293. return result;
  294. }
  295. result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2));
  296. return result;
  297. }
  298. template <class RealType, class Policy>
  299. inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist)
  300. {
  301. BOOST_MATH_STD_USING // for ADL of std functions
  302. static const char* function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)";
  303. RealType shape = dist.shape();
  304. RealType scale = dist.scale();
  305. RealType result = 0;
  306. if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
  307. {
  308. return result;
  309. }
  310. // Only defined for shape >= 0, but is checked by check_inverse_gamma.
  311. result = scale / (shape + 1);
  312. return result;
  313. }
  314. //template <class RealType, class Policy>
  315. //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
  316. //{ // Wikipedia does not define median,
  317. // so rely on default definition quantile(0.5) in derived accessors.
  318. // return result.
  319. //}
  320. template <class RealType, class Policy>
  321. inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist)
  322. {
  323. BOOST_MATH_STD_USING // for ADL of std functions
  324. static const char* function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)";
  325. RealType shape = dist.shape();
  326. RealType scale = dist.scale();
  327. RealType result = 0;
  328. if(false == detail::check_scale(function, scale, &result, Policy()))
  329. {
  330. return result;
  331. }
  332. if((shape <= 3) || !(boost::math::isfinite)(shape))
  333. {
  334. result = policies::raise_domain_error<RealType>(
  335. function,
  336. "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy());
  337. return result;
  338. }
  339. result = (4 * sqrt(shape - 2) ) / (shape - 3);
  340. return result;
  341. }
  342. template <class RealType, class Policy>
  343. inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist)
  344. {
  345. BOOST_MATH_STD_USING // for ADL of std functions
  346. static const char* function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)";
  347. RealType shape = dist.shape();
  348. RealType scale = dist.scale();
  349. RealType result = 0;
  350. if(false == detail::check_scale(function, scale, &result, Policy()))
  351. {
  352. return result;
  353. }
  354. if((shape <= 4) || !(boost::math::isfinite)(shape))
  355. {
  356. result = policies::raise_domain_error<RealType>(
  357. function,
  358. "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
  359. return result;
  360. }
  361. result = (30 * shape - 66) / ((shape - 3) * (shape - 4));
  362. return result;
  363. }
  364. template <class RealType, class Policy>
  365. inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist)
  366. {
  367. static const char* function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)";
  368. RealType shape = dist.shape();
  369. RealType scale = dist.scale();
  370. RealType result = 0;
  371. if(false == detail::check_scale(function, scale, &result, Policy()))
  372. {
  373. return result;
  374. }
  375. if((shape <= 4) || !(boost::math::isfinite)(shape))
  376. {
  377. result = policies::raise_domain_error<RealType>(
  378. function,
  379. "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
  380. return result;
  381. }
  382. return kurtosis_excess(dist) + 3;
  383. }
  384. } // namespace math
  385. } // namespace boost
  386. // This include must be at the end, *after* the accessors
  387. // for this distribution have been defined, in order to
  388. // keep compilers that support two-phase lookup happy.
  389. #include <boost/math/distributions/detail/derived_accessors.hpp>
  390. #endif // BOOST_STATS_INVERSE_GAMMA_HPP