3 // Copyright Paul A. Bristow 2010.
4 // Copyright John Maddock 2010.
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)
9 #ifndef BOOST_STATS_INVERSE_GAMMA_HPP
10 #define BOOST_STATS_INVERSE_GAMMA_HPP
12 // Inverse Gamma Distribution is a two-parameter family
13 // of continuous probability distributions
14 // on the positive real line, which is the distribution of
15 // the reciprocal of a variable distributed according to the gamma distribution.
17 // http://en.wikipedia.org/wiki/Inverse-gamma_distribution
18 // http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html
20 // See also gamma distribution at gamma.hpp:
21 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
22 // http://mathworld.wolfram.com/GammaDistribution.html
23 // http://en.wikipedia.org/wiki/Gamma_distribution
25 #include <boost/math/distributions/fwd.hpp>
26 #include <boost/math/special_functions/gamma.hpp>
27 #include <boost/math/distributions/detail/common_error_handling.hpp>
28 #include <boost/math/distributions/complement.hpp>
32 namespace boost{ namespace math
37 template <class RealType, class Policy>
38 inline bool check_inverse_gamma_shape(
39 const char* function, // inverse_gamma
40 RealType shape, // shape aka alpha
41 RealType* result, // to update, perhaps with NaN
43 { // Sources say shape argument must be > 0
44 // but seems logical to allow shape zero as special case,
45 // returning pdf and cdf zero (but not < 0).
46 // (Functions like mean, variance with other limits on shape are checked
47 // in version including an operator & limit below).
48 if((shape < 0) || !(boost::math::isfinite)(shape))
50 *result = policies::raise_domain_error<RealType>(
52 "Shape parameter is %1%, but must be >= 0 !", shape, pol);
56 } //bool check_inverse_gamma_shape
58 template <class RealType, class Policy>
59 inline bool check_inverse_gamma_x(
62 RealType* result, const Policy& pol)
64 if((x < 0) || !(boost::math::isfinite)(x))
66 *result = policies::raise_domain_error<RealType>(
68 "Random variate is %1% but must be >= 0 !", x, pol);
74 template <class RealType, class Policy>
75 inline bool check_inverse_gamma(
76 const char* function, // TODO swap these over, so shape is first.
77 RealType scale, // scale aka beta
78 RealType shape, // shape aka alpha
79 RealType* result, const Policy& pol)
81 return check_scale(function, scale, result, pol)
82 && check_inverse_gamma_shape(function, shape, result, pol);
83 } // bool check_inverse_gamma
87 template <class RealType = double, class Policy = policies::policy<> >
88 class inverse_gamma_distribution
91 typedef RealType value_type;
92 typedef Policy policy_type;
94 inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1)
95 : m_shape(l_shape), m_scale(l_scale)
98 detail::check_inverse_gamma(
99 "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution",
100 l_scale, l_shape, &result, Policy());
103 RealType shape()const
108 RealType scale()const
116 RealType m_shape; // distribution shape
117 RealType m_scale; // distribution scale
120 typedef inverse_gamma_distribution<double> inverse_gamma;
121 // typedef - but potential clash with name of inverse gamma *function*.
122 // but there is a typedef for gamma
123 // typedef boost::math::gamma_distribution<Type, Policy> gamma;
125 // Allow random variable x to be zero, treated as a special case (unlike some definitions).
127 template <class RealType, class Policy>
128 inline const std::pair<RealType, RealType> range(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
129 { // Range of permissible values for random variable x.
130 using boost::math::tools::max_value;
131 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
134 template <class RealType, class Policy>
135 inline const std::pair<RealType, RealType> support(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
136 { // Range of supported values for random variable x.
137 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
138 using boost::math::tools::max_value;
139 using boost::math::tools::min_value;
140 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
143 template <class RealType, class Policy>
144 inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
146 BOOST_MATH_STD_USING // for ADL of std functions
148 static const char* function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)";
150 RealType shape = dist.shape();
151 RealType scale = dist.scale();
154 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
155 { // distribution parameters bad.
159 { // Treat random variate zero as a special case.
162 else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
167 if(result < tools::min_value<RealType>())
168 return 0; // random variable is infinite or so close as to make no difference.
169 result = gamma_p_derivative(shape, result, Policy()) * scale;
174 // x * x may under or overflow, likewise our result,
175 // so be extra careful about the arithmetic:
176 RealType lim = tools::max_value<RealType>() * x;
178 return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
181 return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
187 // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape);
191 template <class RealType, class Policy>
192 inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
194 BOOST_MATH_STD_USING // for ADL of std functions
196 static const char* function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)";
198 RealType shape = dist.shape();
199 RealType scale = dist.scale();
202 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
203 { // distribution parameters bad.
207 { // Treat zero as a special case.
210 else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
214 result = boost::math::gamma_q(shape, scale / x, Policy());
215 // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma
219 template <class RealType, class Policy>
220 inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)
222 BOOST_MATH_STD_USING // for ADL of std functions
223 using boost::math::gamma_q_inv;
225 static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
227 RealType shape = dist.shape();
228 RealType scale = dist.scale();
231 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
233 if(false == detail::check_probability(function, p, &result, Policy()))
237 return policies::raise_overflow_error<RealType>(function, 0, Policy());
239 result = gamma_q_inv(shape, p, Policy());
240 if((result < 1) && (result * tools::max_value<RealType>() < scale))
241 return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
242 result = scale / result;
246 template <class RealType, class Policy>
247 inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
249 BOOST_MATH_STD_USING // for ADL of std functions
251 static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
253 RealType shape = c.dist.shape();
254 RealType scale = c.dist.scale();
257 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
259 if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy()))
263 return 1; // Avoid division by zero
265 //result = 1. - gamma_q(shape, c.param / scale, Policy());
266 result = gamma_p(shape, scale/c.param, Policy());
270 template <class RealType, class Policy>
271 inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
273 BOOST_MATH_STD_USING // for ADL of std functions
275 static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
277 RealType shape = c.dist.shape();
278 RealType scale = c.dist.scale();
279 RealType q = c.param;
282 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
284 if(false == detail::check_probability(function, q, &result, Policy()))
289 return policies::raise_overflow_error<RealType>(function, 0, Policy());
291 result = gamma_p_inv(shape, q, Policy());
292 if((result < 1) && (result * tools::max_value<RealType>() < scale))
293 return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
294 result = scale / result;
298 template <class RealType, class Policy>
299 inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist)
301 BOOST_MATH_STD_USING // for ADL of std functions
303 static const char* function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)";
305 RealType shape = dist.shape();
306 RealType scale = dist.scale();
310 if(false == detail::check_scale(function, scale, &result, Policy()))
314 if((shape <= 1) || !(boost::math::isfinite)(shape))
316 result = policies::raise_domain_error<RealType>(
318 "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy());
321 result = scale / (shape - 1);
325 template <class RealType, class Policy>
326 inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist)
328 BOOST_MATH_STD_USING // for ADL of std functions
330 static const char* function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)";
332 RealType shape = dist.shape();
333 RealType scale = dist.scale();
336 if(false == detail::check_scale(function, scale, &result, Policy()))
340 if((shape <= 2) || !(boost::math::isfinite)(shape))
342 result = policies::raise_domain_error<RealType>(
344 "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy());
347 result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2));
351 template <class RealType, class Policy>
352 inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist)
354 BOOST_MATH_STD_USING // for ADL of std functions
356 static const char* function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)";
358 RealType shape = dist.shape();
359 RealType scale = dist.scale();
362 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
366 // Only defined for shape >= 0, but is checked by check_inverse_gamma.
367 result = scale / (shape + 1);
371 //template <class RealType, class Policy>
372 //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
373 //{ // Wikipedia does not define median,
374 // so rely on default definition quantile(0.5) in derived accessors.
378 template <class RealType, class Policy>
379 inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist)
381 BOOST_MATH_STD_USING // for ADL of std functions
383 static const char* function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)";
385 RealType shape = dist.shape();
386 RealType scale = dist.scale();
389 if(false == detail::check_scale(function, scale, &result, Policy()))
393 if((shape <= 3) || !(boost::math::isfinite)(shape))
395 result = policies::raise_domain_error<RealType>(
397 "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy());
400 result = (4 * sqrt(shape - 2) ) / (shape - 3);
404 template <class RealType, class Policy>
405 inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist)
407 BOOST_MATH_STD_USING // for ADL of std functions
409 static const char* function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)";
411 RealType shape = dist.shape();
412 RealType scale = dist.scale();
415 if(false == detail::check_scale(function, scale, &result, Policy()))
419 if((shape <= 4) || !(boost::math::isfinite)(shape))
421 result = policies::raise_domain_error<RealType>(
423 "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
426 result = (30 * shape - 66) / ((shape - 3) * (shape - 4));
430 template <class RealType, class Policy>
431 inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist)
433 static const char* function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)";
434 RealType shape = dist.shape();
435 RealType scale = dist.scale();
439 if(false == detail::check_scale(function, scale, &result, Policy()))
443 if((shape <= 4) || !(boost::math::isfinite)(shape))
445 result = policies::raise_domain_error<RealType>(
447 "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
450 return kurtosis_excess(dist) + 3;
456 // This include must be at the end, *after* the accessors
457 // for this distribution have been defined, in order to
458 // keep compilers that support two-phase lookup happy.
459 #include <boost/math/distributions/detail/derived_accessors.hpp>
461 #endif // BOOST_STATS_INVERSE_GAMMA_HPP