]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // inverse_gamma.hpp |
2 | ||
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) | |
8 | ||
9 | #ifndef BOOST_STATS_INVERSE_GAMMA_HPP | |
10 | #define BOOST_STATS_INVERSE_GAMMA_HPP | |
11 | ||
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. | |
16 | ||
17 | // http://en.wikipedia.org/wiki/Inverse-gamma_distribution | |
18 | // http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html | |
19 | ||
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 | |
24 | ||
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> | |
29 | ||
30 | #include <utility> | |
31 | ||
32 | namespace boost{ namespace math | |
33 | { | |
34 | namespace detail | |
35 | { | |
36 | ||
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 | |
42 | const Policy& pol) | |
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)) | |
49 | { | |
50 | *result = policies::raise_domain_error<RealType>( | |
51 | function, | |
52 | "Shape parameter is %1%, but must be >= 0 !", shape, pol); | |
53 | return false; | |
54 | } | |
55 | return true; | |
56 | } //bool check_inverse_gamma_shape | |
57 | ||
58 | template <class RealType, class Policy> | |
59 | inline bool check_inverse_gamma_x( | |
60 | const char* function, | |
61 | RealType const& x, | |
62 | RealType* result, const Policy& pol) | |
63 | { | |
64 | if((x < 0) || !(boost::math::isfinite)(x)) | |
65 | { | |
66 | *result = policies::raise_domain_error<RealType>( | |
67 | function, | |
68 | "Random variate is %1% but must be >= 0 !", x, pol); | |
69 | return false; | |
70 | } | |
71 | return true; | |
72 | } | |
73 | ||
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) | |
80 | { | |
81 | return check_scale(function, scale, result, pol) | |
82 | && check_inverse_gamma_shape(function, shape, result, pol); | |
83 | } // bool check_inverse_gamma | |
84 | ||
85 | } // namespace detail | |
86 | ||
87 | template <class RealType = double, class Policy = policies::policy<> > | |
88 | class inverse_gamma_distribution | |
89 | { | |
90 | public: | |
91 | typedef RealType value_type; | |
92 | typedef Policy policy_type; | |
93 | ||
94 | inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1) | |
95 | : m_shape(l_shape), m_scale(l_scale) | |
96 | { | |
97 | RealType result; | |
98 | detail::check_inverse_gamma( | |
99 | "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution", | |
100 | l_scale, l_shape, &result, Policy()); | |
101 | } | |
102 | ||
103 | RealType shape()const | |
104 | { | |
105 | return m_shape; | |
106 | } | |
107 | ||
108 | RealType scale()const | |
109 | { | |
110 | return m_scale; | |
111 | } | |
112 | private: | |
113 | // | |
114 | // Data members: | |
115 | // | |
116 | RealType m_shape; // distribution shape | |
117 | RealType m_scale; // distribution scale | |
118 | }; | |
119 | ||
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; | |
124 | ||
125 | // Allow random variable x to be zero, treated as a special case (unlike some definitions). | |
126 | ||
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>()); | |
132 | } | |
133 | ||
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>()); | |
141 | } | |
142 | ||
143 | template <class RealType, class Policy> | |
144 | inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x) | |
145 | { | |
146 | BOOST_MATH_STD_USING // for ADL of std functions | |
147 | ||
148 | static const char* function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)"; | |
149 | ||
150 | RealType shape = dist.shape(); | |
151 | RealType scale = dist.scale(); | |
152 | ||
153 | RealType result = 0; | |
154 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
155 | { // distribution parameters bad. | |
156 | return result; | |
157 | } | |
158 | if(x == 0) | |
159 | { // Treat random variate zero as a special case. | |
160 | return 0; | |
161 | } | |
162 | else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy())) | |
163 | { // x bad. | |
164 | return result; | |
165 | } | |
166 | result = scale / x; | |
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; | |
170 | if(0 != result) | |
171 | { | |
172 | if(x < 0) | |
173 | { | |
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; | |
177 | if(lim < result) | |
178 | return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy()); | |
179 | result /= x; | |
180 | if(lim < result) | |
181 | return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy()); | |
182 | result /= x; | |
183 | } | |
184 | result /= (x * x); | |
185 | } | |
186 | // better than naive | |
187 | // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape); | |
188 | return result; | |
189 | ||
190 | ||
191 | template <class RealType, class Policy> | |
192 | inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x) | |
193 | { | |
194 | BOOST_MATH_STD_USING // for ADL of std functions | |
195 | ||
196 | static const char* function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)"; | |
197 | ||
198 | RealType shape = dist.shape(); | |
199 | RealType scale = dist.scale(); | |
200 | ||
201 | RealType result = 0; | |
202 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
203 | { // distribution parameters bad. | |
204 | return result; | |
205 | } | |
206 | if (x == 0) | |
207 | { // Treat zero as a special case. | |
208 | return 0; | |
209 | } | |
210 | else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy())) | |
211 | { // x bad | |
212 | return result; | |
213 | } | |
214 | result = boost::math::gamma_q(shape, scale / x, Policy()); | |
215 | // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma | |
216 | return result; | |
217 | } // cdf | |
218 | ||
219 | template <class RealType, class Policy> | |
220 | inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p) | |
221 | { | |
222 | BOOST_MATH_STD_USING // for ADL of std functions | |
223 | using boost::math::gamma_q_inv; | |
224 | ||
225 | static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)"; | |
226 | ||
227 | RealType shape = dist.shape(); | |
228 | RealType scale = dist.scale(); | |
229 | ||
230 | RealType result = 0; | |
231 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
232 | return result; | |
233 | if(false == detail::check_probability(function, p, &result, Policy())) | |
234 | return result; | |
235 | if(p == 1) | |
236 | { | |
237 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
238 | } | |
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; | |
243 | return result; | |
244 | } | |
245 | ||
246 | template <class RealType, class Policy> | |
247 | inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c) | |
248 | { | |
249 | BOOST_MATH_STD_USING // for ADL of std functions | |
250 | ||
251 | static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)"; | |
252 | ||
253 | RealType shape = c.dist.shape(); | |
254 | RealType scale = c.dist.scale(); | |
255 | ||
256 | RealType result = 0; | |
257 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
258 | return result; | |
259 | if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy())) | |
260 | return result; | |
261 | ||
262 | if(c.param == 0) | |
263 | return 1; // Avoid division by zero | |
264 | ||
265 | //result = 1. - gamma_q(shape, c.param / scale, Policy()); | |
266 | result = gamma_p(shape, scale/c.param, Policy()); | |
267 | return result; | |
268 | } | |
269 | ||
270 | template <class RealType, class Policy> | |
271 | inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c) | |
272 | { | |
273 | BOOST_MATH_STD_USING // for ADL of std functions | |
274 | ||
275 | static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)"; | |
276 | ||
277 | RealType shape = c.dist.shape(); | |
278 | RealType scale = c.dist.scale(); | |
279 | RealType q = c.param; | |
280 | ||
281 | RealType result = 0; | |
282 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
283 | return result; | |
284 | if(false == detail::check_probability(function, q, &result, Policy())) | |
285 | return result; | |
286 | ||
287 | if(q == 0) | |
288 | { | |
289 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
290 | } | |
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; | |
295 | return result; | |
296 | } | |
297 | ||
298 | template <class RealType, class Policy> | |
299 | inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist) | |
300 | { | |
301 | BOOST_MATH_STD_USING // for ADL of std functions | |
302 | ||
303 | static const char* function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)"; | |
304 | ||
305 | RealType shape = dist.shape(); | |
306 | RealType scale = dist.scale(); | |
307 | ||
308 | RealType result = 0; | |
309 | ||
310 | if(false == detail::check_scale(function, scale, &result, Policy())) | |
311 | { | |
312 | return result; | |
313 | } | |
314 | if((shape <= 1) || !(boost::math::isfinite)(shape)) | |
315 | { | |
316 | result = policies::raise_domain_error<RealType>( | |
317 | function, | |
318 | "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy()); | |
319 | return result; | |
320 | } | |
321 | result = scale / (shape - 1); | |
322 | return result; | |
323 | } // mean | |
324 | ||
325 | template <class RealType, class Policy> | |
326 | inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist) | |
327 | { | |
328 | BOOST_MATH_STD_USING // for ADL of std functions | |
329 | ||
330 | static const char* function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)"; | |
331 | ||
332 | RealType shape = dist.shape(); | |
333 | RealType scale = dist.scale(); | |
334 | ||
335 | RealType result = 0; | |
336 | if(false == detail::check_scale(function, scale, &result, Policy())) | |
337 | { | |
338 | return result; | |
339 | } | |
340 | if((shape <= 2) || !(boost::math::isfinite)(shape)) | |
341 | { | |
342 | result = policies::raise_domain_error<RealType>( | |
343 | function, | |
344 | "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy()); | |
345 | return result; | |
346 | } | |
347 | result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2)); | |
348 | return result; | |
349 | } | |
350 | ||
351 | template <class RealType, class Policy> | |
352 | inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist) | |
353 | { | |
354 | BOOST_MATH_STD_USING // for ADL of std functions | |
355 | ||
356 | static const char* function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)"; | |
357 | ||
358 | RealType shape = dist.shape(); | |
359 | RealType scale = dist.scale(); | |
360 | ||
361 | RealType result = 0; | |
362 | if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy())) | |
363 | { | |
364 | return result; | |
365 | } | |
366 | // Only defined for shape >= 0, but is checked by check_inverse_gamma. | |
367 | result = scale / (shape + 1); | |
368 | return result; | |
369 | } | |
370 | ||
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. | |
375 | // return result. | |
376 | //} | |
377 | ||
378 | template <class RealType, class Policy> | |
379 | inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist) | |
380 | { | |
381 | BOOST_MATH_STD_USING // for ADL of std functions | |
382 | ||
383 | static const char* function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)"; | |
384 | ||
385 | RealType shape = dist.shape(); | |
386 | RealType scale = dist.scale(); | |
387 | RealType result = 0; | |
388 | ||
389 | if(false == detail::check_scale(function, scale, &result, Policy())) | |
390 | { | |
391 | return result; | |
392 | } | |
393 | if((shape <= 3) || !(boost::math::isfinite)(shape)) | |
394 | { | |
395 | result = policies::raise_domain_error<RealType>( | |
396 | function, | |
397 | "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy()); | |
398 | return result; | |
399 | } | |
400 | result = (4 * sqrt(shape - 2) ) / (shape - 3); | |
401 | return result; | |
402 | } | |
403 | ||
404 | template <class RealType, class Policy> | |
405 | inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist) | |
406 | { | |
407 | BOOST_MATH_STD_USING // for ADL of std functions | |
408 | ||
409 | static const char* function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)"; | |
410 | ||
411 | RealType shape = dist.shape(); | |
412 | RealType scale = dist.scale(); | |
413 | ||
414 | RealType result = 0; | |
415 | if(false == detail::check_scale(function, scale, &result, Policy())) | |
416 | { | |
417 | return result; | |
418 | } | |
419 | if((shape <= 4) || !(boost::math::isfinite)(shape)) | |
420 | { | |
421 | result = policies::raise_domain_error<RealType>( | |
422 | function, | |
423 | "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy()); | |
424 | return result; | |
425 | } | |
426 | result = (30 * shape - 66) / ((shape - 3) * (shape - 4)); | |
427 | return result; | |
428 | } | |
429 | ||
430 | template <class RealType, class Policy> | |
431 | inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist) | |
432 | { | |
433 | static const char* function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)"; | |
434 | RealType shape = dist.shape(); | |
435 | RealType scale = dist.scale(); | |
436 | ||
437 | RealType result = 0; | |
438 | ||
439 | if(false == detail::check_scale(function, scale, &result, Policy())) | |
440 | { | |
441 | return result; | |
442 | } | |
443 | if((shape <= 4) || !(boost::math::isfinite)(shape)) | |
444 | { | |
445 | result = policies::raise_domain_error<RealType>( | |
446 | function, | |
447 | "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy()); | |
448 | return result; | |
449 | } | |
450 | return kurtosis_excess(dist) + 3; | |
451 | } | |
452 | ||
453 | } // namespace math | |
454 | } // namespace boost | |
455 | ||
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> | |
460 | ||
461 | #endif // BOOST_STATS_INVERSE_GAMMA_HPP |