statrs/distribution/
inverse_gamma.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::gamma;
3use crate::statistics::*;
4use crate::{Result, StatsError};
5use rand::Rng;
6use std::f64;
7
8/// Implements the [Inverse
9/// Gamma](https://en.wikipedia.org/wiki/Inverse-gamma_distribution)
10/// distribution
11///
12/// # Examples
13///
14/// ```
15/// use statrs::distribution::{InverseGamma, Continuous};
16/// use statrs::statistics::Distribution;
17/// use statrs::prec;
18///
19/// let n = InverseGamma::new(1.1, 0.1).unwrap();
20/// assert!(prec::almost_eq(n.mean().unwrap(), 1.0, 1e-14));
21/// assert_eq!(n.pdf(1.0), 0.07554920138253064);
22/// ```
23#[derive(Debug, Copy, Clone, PartialEq)]
24pub struct InverseGamma {
25    shape: f64,
26    rate: f64,
27}
28
29impl InverseGamma {
30    /// Constructs a new inverse gamma distribution with a shape (α)
31    /// of `shape` and a rate (β) of `rate`
32    ///
33    /// # Errors
34    ///
35    /// Returns an error if `shape` or `rate` are `NaN`.
36    /// Also returns an error if `shape` or `rate` are not in `(0, +inf)`
37    ///
38    /// # Examples
39    ///
40    /// ```
41    /// use statrs::distribution::InverseGamma;
42    ///
43    /// let mut result = InverseGamma::new(3.0, 1.0);
44    /// assert!(result.is_ok());
45    ///
46    /// result = InverseGamma::new(0.0, 0.0);
47    /// assert!(result.is_err());
48    /// ```
49    pub fn new(shape: f64, rate: f64) -> Result<InverseGamma> {
50        let is_nan = shape.is_nan() || rate.is_nan();
51        match (shape, rate, is_nan) {
52            (_, _, true) => Err(StatsError::BadParams),
53            (_, _, false) if shape <= 0.0 || rate <= 0.0 => Err(StatsError::BadParams),
54            (_, _, false) if shape.is_infinite() || rate.is_infinite() => {
55                Err(StatsError::BadParams)
56            }
57            (_, _, false) => Ok(InverseGamma { shape, rate }),
58        }
59    }
60
61    /// Returns the shape (α) of the inverse gamma distribution
62    ///
63    /// # Examples
64    ///
65    /// ```
66    /// use statrs::distribution::InverseGamma;
67    ///
68    /// let n = InverseGamma::new(3.0, 1.0).unwrap();
69    /// assert_eq!(n.shape(), 3.0);
70    /// ```
71    pub fn shape(&self) -> f64 {
72        self.shape
73    }
74
75    /// Returns the rate (β) of the inverse gamma distribution
76    ///
77    /// # Examples
78    ///
79    /// ```
80    /// use statrs::distribution::InverseGamma;
81    ///
82    /// let n = InverseGamma::new(3.0, 1.0).unwrap();
83    /// assert_eq!(n.rate(), 1.0);
84    /// ```
85    pub fn rate(&self) -> f64 {
86        self.rate
87    }
88}
89
90impl ::rand::distributions::Distribution<f64> for InverseGamma {
91    fn sample<R: Rng + ?Sized>(&self, r: &mut R) -> f64 {
92        1.0 / super::gamma::sample_unchecked(r, self.shape, self.rate)
93    }
94}
95
96impl ContinuousCDF<f64, f64> for InverseGamma {
97    /// Calculates the cumulative distribution function for the inverse gamma
98    /// distribution at `x`
99    ///
100    /// # Formula
101    ///
102    /// ```ignore
103    /// Γ(α, β / x) / Γ(α)
104    /// ```
105    ///
106    /// where the numerator is the upper incomplete gamma function,
107    /// the denominator is the gamma function, `α` is the shape,
108    /// and `β` is the rate
109    fn cdf(&self, x: f64) -> f64 {
110        if x <= 0.0 {
111            0.0
112        } else if x.is_infinite() {
113            1.0
114        } else {
115            gamma::gamma_ur(self.shape, self.rate / x)
116        }
117    }
118
119    /// Calculates the survival function for the inverse gamma
120    /// distribution at `x`
121    ///
122    /// # Formula
123    ///
124    /// ```ignore
125    /// Γ(α, β / x) / Γ(α)
126    /// ```
127    ///
128    /// where the numerator is the lower incomplete gamma function,
129    /// the denominator is the gamma function, `α` is the shape,
130    /// and `β` is the rate
131    fn sf(&self, x: f64) -> f64 {
132        if x <= 0.0 {
133            1.0
134        } else if x.is_infinite() {
135            0.0
136        } else {
137            gamma::gamma_lr(self.shape, self.rate / x)
138        }
139    }
140}
141
142impl Min<f64> for InverseGamma {
143    /// Returns the minimum value in the domain of the
144    /// inverse gamma distribution representable by a double precision
145    /// float
146    ///
147    /// # Formula
148    ///
149    /// ```ignore
150    /// 0
151    /// ```
152    fn min(&self) -> f64 {
153        0.0
154    }
155}
156
157impl Max<f64> for InverseGamma {
158    /// Returns the maximum value in the domain of the
159    /// inverse gamma distribution representable by a double precision
160    /// float
161    ///
162    /// # Formula
163    ///
164    /// ```ignore
165    /// INF
166    /// ```
167    fn max(&self) -> f64 {
168        f64::INFINITY
169    }
170}
171
172impl Distribution<f64> for InverseGamma {
173    /// Returns the mean of the inverse distribution
174    ///
175    /// # None
176    ///
177    /// If `shape <= 1.0`
178    ///
179    /// # Formula
180    ///
181    /// ```ignore
182    /// β / (α - 1)
183    /// ```
184    ///
185    /// where `α` is the shape and `β` is the rate
186    fn mean(&self) -> Option<f64> {
187        if self.shape <= 1.0 {
188            None
189        } else {
190            Some(self.rate / (self.shape - 1.0))
191        }
192    }
193    /// Returns the variance of the inverse gamma distribution
194    ///
195    /// # None
196    ///
197    /// If `shape <= 2.0`
198    ///
199    /// # Formula
200    ///
201    /// ```ignore
202    /// β^2 / ((α - 1)^2 * (α - 2))
203    /// ```
204    ///
205    /// where `α` is the shape and `β` is the rate
206    fn variance(&self) -> Option<f64> {
207        if self.shape <= 2.0 {
208            None
209        } else {
210            let val = self.rate * self.rate
211                / ((self.shape - 1.0) * (self.shape - 1.0) * (self.shape - 2.0));
212            Some(val)
213        }
214    }
215    /// Returns the entropy of the inverse gamma distribution
216    ///
217    /// # Formula
218    ///
219    /// ```ignore
220    /// α + ln(β * Γ(α)) - (1 + α) * ψ(α)
221    /// ```
222    ///
223    /// where `α` is the shape, `β` is the rate, `Γ` is the gamma function,
224    /// and `ψ` is the digamma function
225    fn entropy(&self) -> Option<f64> {
226        let entr = self.shape + self.rate.ln() + gamma::ln_gamma(self.shape)
227            - (1.0 + self.shape) * gamma::digamma(self.shape);
228        Some(entr)
229    }
230    /// Returns the skewness of the inverse gamma distribution
231    ///
232    /// # None
233    ///
234    /// If `shape <= 3`
235    ///
236    /// # Formula
237    ///
238    /// ```ignore
239    /// 4 * sqrt(α - 2) / (α - 3)
240    /// ```
241    ///
242    /// where `α` is the shape
243    fn skewness(&self) -> Option<f64> {
244        if self.shape <= 3.0 {
245            None
246        } else {
247            Some(4.0 * (self.shape - 2.0).sqrt() / (self.shape - 3.0))
248        }
249    }
250}
251
252impl Mode<Option<f64>> for InverseGamma {
253    /// Returns the mode of the inverse gamma distribution
254    ///
255    /// # Formula
256    ///
257    /// ```ignore
258    /// β / (α + 1)
259    /// ```
260    ///
261    /// /// where `α` is the shape and `β` is the rate
262    fn mode(&self) -> Option<f64> {
263        Some(self.rate / (self.shape + 1.0))
264    }
265}
266
267impl Continuous<f64, f64> for InverseGamma {
268    /// Calculates the probability density function for the
269    /// inverse gamma distribution at `x`
270    ///
271    /// # Formula
272    ///
273    /// ```ignore
274    /// (β^α / Γ(α)) * x^(-α - 1) * e^(-β / x)
275    /// ```
276    ///
277    /// where `α` is the shape, `β` is the rate, and `Γ` is the gamma function
278    fn pdf(&self, x: f64) -> f64 {
279        if x <= 0.0 || x.is_infinite() {
280            0.0
281        } else if ulps_eq!(self.shape, 1.0) {
282            self.rate / (x * x) * (-self.rate / x).exp()
283        } else {
284            self.rate.powf(self.shape) * x.powf(-self.shape - 1.0) * (-self.rate / x).exp()
285                / gamma::gamma(self.shape)
286        }
287    }
288
289    /// Calculates the probability density function for the
290    /// inverse gamma distribution at `x`
291    ///
292    /// # Formula
293    ///
294    /// ```ignore
295    /// ln((β^α / Γ(α)) * x^(-α - 1) * e^(-β / x))
296    /// ```
297    ///
298    /// where `α` is the shape, `β` is the rate, and `Γ` is the gamma function
299    fn ln_pdf(&self, x: f64) -> f64 {
300        self.pdf(x).ln()
301    }
302}
303
304#[rustfmt::skip]
305#[cfg(all(test, feature = "nightly"))]
306mod tests {
307    use crate::statistics::*;
308    use crate::distribution::{ContinuousCDF, Continuous, InverseGamma};
309    use crate::distribution::internal::*;
310    use crate::consts::ACC;
311
312    fn try_create(shape: f64, rate: f64) -> InverseGamma {
313        let n = InverseGamma::new(shape, rate);
314        assert!(n.is_ok());
315        n.unwrap()
316    }
317
318    fn create_case(shape: f64, rate: f64) {
319        let n = try_create(shape, rate);
320        assert_eq!(shape, n.shape());
321        assert_eq!(rate, n.rate());
322    }
323
324    fn bad_create_case(shape: f64, rate: f64) {
325        let n = InverseGamma::new(shape, rate);
326        assert!(n.is_err());
327    }
328
329    fn get_value<F>(shape: f64, rate: f64, eval: F) -> f64
330        where F: Fn(InverseGamma) -> f64
331    {
332        let n = try_create(shape, rate);
333        eval(n)
334    }
335
336    fn test_case<F>(shape: f64, rate: f64, expected: f64, eval: F)
337        where F: Fn(InverseGamma) -> f64
338    {
339        let x = get_value(shape, rate, eval);
340        assert_eq!(expected, x);
341    }
342
343    fn test_almost<F>(shape: f64, rate: f64, expected: f64, acc: f64, eval: F)
344        where F: Fn(InverseGamma) -> f64
345    {
346        let x = get_value(shape, rate, eval);
347        assert_almost_eq!(expected, x, acc);
348    }
349
350    #[test]
351    fn test_create() {
352        create_case(0.1, 0.1);
353        create_case(1.0, 1.0);
354    }
355
356    #[test]
357    fn test_bad_create() {
358        bad_create_case(0.0, 1.0);
359        bad_create_case(-1.0, 1.0);
360        bad_create_case(-100.0, 1.0);
361        bad_create_case(f64::NEG_INFINITY, 1.0);
362        bad_create_case(f64::NAN, 1.0);
363        bad_create_case(1.0, 0.0);
364        bad_create_case(1.0, -1.0);
365        bad_create_case(1.0, -100.0);
366        bad_create_case(1.0, f64::NEG_INFINITY);
367        bad_create_case(1.0, f64::NAN);
368        bad_create_case(f64::INFINITY, 1.0);
369        bad_create_case(1.0, f64::INFINITY);
370        bad_create_case(f64::INFINITY, f64::INFINITY);
371    }
372
373    #[test]
374    fn test_mean() {
375        let mean = |x: InverseGamma| x.mean().unwrap();
376        test_almost(1.1, 0.1, 1.0, 1e-14, mean);
377        test_almost(1.1, 1.0, 10.0, 1e-14, mean);
378    }
379
380    #[test]
381    #[should_panic]
382    fn test_mean_with_shape_lte_1() {
383        let mean = |x: InverseGamma| x.mean().unwrap();
384        get_value(0.1, 0.1, mean);
385    }
386
387    #[test]
388    fn test_variance() {
389        let variance = |x: InverseGamma| x.variance().unwrap();
390        test_almost(2.1, 0.1, 0.08264462809917355371901, 1e-15, variance);
391        test_almost(2.1, 1.0, 8.264462809917355371901, 1e-13, variance);
392    }
393
394    #[test]
395    #[should_panic]
396    fn test_variance_with_shape_lte_2() {
397        let variance = |x: InverseGamma| x.variance().unwrap();
398        get_value(0.1, 0.1, variance);
399    }
400
401    #[test]
402    fn test_entropy() {
403        let entropy = |x: InverseGamma| x.entropy().unwrap();
404        test_almost(0.1, 0.1, 11.51625799319234475054, 1e-14, entropy);
405        test_almost(1.0, 1.0, 2.154431329803065721213, 1e-14, entropy);
406    }
407
408    #[test]
409    fn test_skewness() {
410        let skewness = |x: InverseGamma| x.skewness().unwrap();
411        test_almost(3.1, 0.1, 41.95235392680606187966, 1e-13, skewness);
412        test_almost(3.1, 1.0, 41.95235392680606187966, 1e-13, skewness);
413        test_case(5.0, 0.1, 3.464101615137754587055, skewness);
414    }
415
416    #[test]
417    #[should_panic]
418    fn test_skewness_with_shape_lte_3() {
419        let skewness = |x: InverseGamma| x.skewness().unwrap();
420        get_value(0.1, 0.1, skewness);
421    }
422
423    #[test]
424    fn test_mode() {
425        let mode = |x: InverseGamma| x.mode().unwrap();
426        test_case(0.1, 0.1, 0.09090909090909090909091, mode);
427        test_case(1.0, 1.0, 0.5, mode);
428    }
429
430    #[test]
431    fn test_min_max() {
432        let min = |x: InverseGamma| x.min();
433        let max = |x: InverseGamma| x.max();
434        test_case(1.0, 1.0, 0.0, min);
435        test_case(1.0, 1.0, f64::INFINITY, max);
436    }
437
438    #[test]
439    fn test_pdf() {
440        let pdf = |arg: f64| move |x: InverseGamma| x.pdf(arg);
441        test_almost(0.1, 0.1, 0.0628591853882328004197, 1e-15, pdf(1.2));
442        test_almost(0.1, 1.0, 0.0297426109178248997426, 1e-15, pdf(2.0));
443        test_case(1.0, 0.1, 0.04157808822362745501024, pdf(1.5));
444        test_case(1.0, 1.0, 0.3018043114632487660842, pdf(1.2));
445    }
446
447    #[test]
448    fn test_ln_pdf() {
449        let ln_pdf = |arg: f64| move |x: InverseGamma| x.ln_pdf(arg);
450        test_almost(0.1, 0.1, 0.0628591853882328004197f64.ln(), 1e-15, ln_pdf(1.2));
451        test_almost(0.1, 1.0, 0.0297426109178248997426f64.ln(), 1e-15, ln_pdf(2.0));
452        test_case(1.0, 0.1, 0.04157808822362745501024f64.ln(), ln_pdf(1.5));
453        test_case(1.0, 1.0, 0.3018043114632487660842f64.ln(), ln_pdf(1.2));
454    }
455
456    #[test]
457    fn test_cdf() {
458        let cdf = |arg: f64| move |x: InverseGamma| x.cdf(arg);
459        test_almost(0.1, 0.1, 0.1862151961946054271994, 1e-14, cdf(1.2));
460        test_almost(0.1, 1.0, 0.05859755410986647796141, 1e-14, cdf(2.0));
461        test_case(1.0, 0.1, 0.9355069850316177377304, cdf(1.5));
462        test_almost(1.0, 1.0, 0.4345982085070782231613, 1e-14, cdf(1.2));
463    }
464
465
466    #[test]
467    fn test_sf() {
468        let sf = |arg: f64| move |x: InverseGamma| x.sf(arg);
469        test_almost(0.1, 0.1, 0.8137848038053936, 1e-14, sf(1.2));
470        test_almost(0.1, 1.0, 0.9414024458901327, 1e-14, sf(2.0));
471        test_almost(1.0, 0.1, 0.0644930149683822, 1e-14, sf(1.5));
472        test_almost(1.0, 1.0, 0.565401791492922, 1e-14, sf(1.2));
473    }
474
475    #[test]
476    fn test_continuous() {
477        test::check_continuous_distribution(&try_create(1.0, 0.5), 0.0, 100.0);
478        test::check_continuous_distribution(&try_create(9.0, 2.0), 0.0, 100.0);
479    }
480}
statrs/distribution/inverse_gamma.rs

statrs/distribution/
inverse_gamma.rs
