inverse_gamma.hpp 17 KB

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