statrs/distribution/
inverse_gamma.rs

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

statrs/distribution/
inverse_gamma.rs
