blob: 252a319d877e9c3020c83af795686939211147e9 [file] [log] [blame]
 use crate::{Distribution, InverseGaussian, Standard, StandardNormal}; use num_traits::Float; use rand::Rng; /// Error type returned from `NormalInverseGaussian::new` #[derive(Debug, PartialEq)] pub enum Error { /// `alpha <= 0` or `nan`. AlphaNegativeOrNull, /// `|beta| >= alpha` or `nan`. AbsoluteBetaNotLessThanAlpha, } /// The [normal-inverse Gaussian distribution](https://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution) #[derive(Debug)] pub struct NormalInverseGaussian where F: Float, StandardNormal: Distribution, Standard: Distribution, { alpha: F, beta: F, inverse_gaussian: InverseGaussian, } impl NormalInverseGaussian where F: Float, StandardNormal: Distribution, Standard: Distribution, { /// Construct a new `NormalInverseGaussian` distribution with the given alpha (tail heaviness) and /// beta (asymmetry) parameters. pub fn new(alpha: F, beta: F) -> Result, Error> { if !(alpha > F::zero()) { return Err(Error::AlphaNegativeOrNull); } if !(beta.abs() < alpha) { return Err(Error::AbsoluteBetaNotLessThanAlpha); } let gamma = (alpha * alpha - beta * beta).sqrt(); let mu = F::one() / gamma; let inverse_gaussian = InverseGaussian::new(mu, F::one()).unwrap(); Ok(Self { alpha, beta, inverse_gaussian, }) } } impl Distribution for NormalInverseGaussian where F: Float, StandardNormal: Distribution, Standard: Distribution, { fn sample(&self, rng: &mut R) -> F where R: Rng + ?Sized { let inv_gauss = rng.sample(&self.inverse_gaussian); self.beta * inv_gauss + inv_gauss.sqrt() * rng.sample(StandardNormal) } } #[cfg(test)] mod tests { use super::*; #[test] fn test_normal_inverse_gaussian() { let norm_inv_gauss = NormalInverseGaussian::new(2.0, 1.0).unwrap(); let mut rng = crate::test::rng(210); for _ in 0..1000 { norm_inv_gauss.sample(&mut rng); } } #[test] fn test_normal_inverse_gaussian_invalid_param() { assert!(NormalInverseGaussian::new(-1.0, 1.0).is_err()); assert!(NormalInverseGaussian::new(-1.0, -1.0).is_err()); assert!(NormalInverseGaussian::new(1.0, 2.0).is_err()); assert!(NormalInverseGaussian::new(2.0, 1.0).is_ok()); } }