rand_distr/
normal_inverse_gaussian.rs

1use crate::{Distribution, InverseGaussian, InverseGaussianError, StandardNormal, StandardUniform};
2use core::fmt;
3use num_traits::Float;
4use rand::{Rng, RngExt};
5
6/// Error type returned from [`NormalInverseGaussian::new`]
7#[derive(Debug, Clone, Copy, PartialEq, Eq)]
8pub enum Error {
9    /// `alpha <= 0` or `nan`.
10    AlphaNegativeOrNull,
11    /// `alpha` is `inf` (or, if subnormals are disabled, too close to the maximum finite float value)
12    AlphaInfinite,
13    /// `|beta| >= alpha` or `nan`.
14    AbsoluteBetaNotLessThanAlpha,
15}
16
17impl fmt::Display for Error {
18    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
19        f.write_str(match self {
20            Error::AlphaNegativeOrNull => {
21                "alpha <= 0 or is NaN in normal inverse Gaussian distribution"
22            }
23            Error::AlphaInfinite => {
24                "alpha is +infinity (or too close to the maximum finite value, if subnormal numbers are not supported) in normal inverse Gaussian distribution"
25            }
26            Error::AbsoluteBetaNotLessThanAlpha => {
27                "|beta| >= alpha or is NaN in normal inverse Gaussian distribution"
28            }
29        })
30    }
31}
32
33#[cfg(feature = "std")]
34impl std::error::Error for Error {}
35
36/// The [normal-inverse Gaussian distribution](https://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution) `NIG(α, β)`.
37///
38/// This is a continuous probability distribution with two parameters,
39/// `α` (`alpha`) and `β` (`beta`), defined in `(-∞, ∞)`.
40/// It is also known as the normal-Wald distribution.
41///
42/// # Plot
43///
44/// The following plot shows the normal-inverse Gaussian distribution with various values of `α` and `β`.
45///
46/// ![Normal-inverse Gaussian distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/normal_inverse_gaussian.svg)
47///
48/// # Example
49/// ```
50/// use rand_distr::{NormalInverseGaussian, Distribution};
51///
52/// let norm_inv_gauss = NormalInverseGaussian::new(2.0, 1.0).unwrap();
53/// let v = norm_inv_gauss.sample(&mut rand::rng());
54/// println!("{} is from a normal-inverse Gaussian(2, 1) distribution", v);
55/// ```
56#[derive(Debug, Clone, Copy, PartialEq)]
57#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
58pub struct NormalInverseGaussian<F>
59where
60    F: Float,
61    StandardNormal: Distribution<F>,
62    StandardUniform: Distribution<F>,
63{
64    beta: F,
65    inverse_gaussian: InverseGaussian<F>,
66}
67
68impl<F> NormalInverseGaussian<F>
69where
70    F: Float,
71    StandardNormal: Distribution<F>,
72    StandardUniform: Distribution<F>,
73{
74    /// Construct a new `NormalInverseGaussian` distribution with the given alpha (tail heaviness) and
75    /// beta (asymmetry) parameters.
76    ///
77    /// Note: If subnormal numbers are not supported or are disabled, this sampler may panic or produce
78    /// incorrect output if alpha is close to the maximum finite float value.
79    pub fn new(alpha: F, beta: F) -> Result<NormalInverseGaussian<F>, Error> {
80        if !(alpha > F::zero()) {
81            return Err(Error::AlphaNegativeOrNull);
82        }
83
84        if !(beta.abs() < alpha) {
85            return Err(Error::AbsoluteBetaNotLessThanAlpha);
86        }
87        // Note: this calculation method for gamma = sqrt(alpha * alpha - beta * beta)
88        // avoids overflow if alpha is large, ensuring gamma <= alpha, which implies
89        // (assuming IEEE754 with subnormals) mu = 1.0 / gamma >= 1 / F::max_value() > 0.
90        let r = beta / alpha;
91        let gamma = alpha * (F::one() - r * r).sqrt();
92        let mu = F::one() / gamma;
93        let inverse_gaussian = InverseGaussian::new(mu, F::one()).map_err(|x| match x {
94            InverseGaussianError::MeanNegativeOrNull => Error::AlphaInfinite,
95            InverseGaussianError::ShapeNegativeOrNull => unreachable!(),
96        })?;
97
98        Ok(Self {
99            beta,
100            inverse_gaussian,
101        })
102    }
103}
104
105impl<F> Distribution<F> for NormalInverseGaussian<F>
106where
107    F: Float,
108    StandardNormal: Distribution<F>,
109    StandardUniform: Distribution<F>,
110{
111    fn sample<R>(&self, rng: &mut R) -> F
112    where
113        R: Rng + ?Sized,
114    {
115        let inv_gauss = rng.sample(self.inverse_gaussian);
116
117        self.beta * inv_gauss + inv_gauss.sqrt() * rng.sample(StandardNormal)
118    }
119}
120
121#[cfg(test)]
122mod tests {
123    use super::*;
124
125    #[test]
126    fn test_normal_inverse_gaussian() {
127        let norm_inv_gauss = NormalInverseGaussian::new(2.0, 1.0).unwrap();
128        let mut rng = crate::test::rng(210);
129        for _ in 0..1000 {
130            norm_inv_gauss.sample(&mut rng);
131        }
132    }
133
134    #[test]
135    fn test_normal_inverse_gaussian_invalid_param() {
136        assert!(NormalInverseGaussian::new(-1.0, 1.0).is_err());
137        assert!(NormalInverseGaussian::new(-1.0, -1.0).is_err());
138        assert!(NormalInverseGaussian::new(1.0, 2.0).is_err());
139        assert!(NormalInverseGaussian::new(f64::INFINITY, 2.0).is_err());
140        assert!(NormalInverseGaussian::new(2.0, 1.0).is_ok());
141    }
142
143    #[test]
144    fn normal_inverse_gaussian_distributions_can_be_compared() {
145        assert_eq!(
146            NormalInverseGaussian::new(1.0, 2.0),
147            NormalInverseGaussian::new(1.0, 2.0)
148        );
149    }
150}
rand_distr/normal_inverse_gaussian.rs

rand_distr/
normal_inverse_gaussian.rs
