ferray_random/distributions/
exponential.rs

1// ferray-random: Exponential distribution sampling — standard_exponential, exponential
2
3use ferray_core::{Array, FerrayError, Ix1};
4
5use crate::bitgen::BitGenerator;
6use crate::generator::{Generator, generate_vec, vec_to_array1};
7
8/// Generate a single standard exponential variate (rate=1) via inverse CDF.
9pub(crate) fn standard_exponential_single<B: BitGenerator>(bg: &mut B) -> f64 {
10    loop {
11        let u = bg.next_f64();
12        if u > f64::EPSILON {
13            return -u.ln();
14        }
15    }
16}
17
18impl<B: BitGenerator> Generator<B> {
19    /// Generate an array of standard exponential (rate=1, scale=1) variates.
20    ///
21    /// Uses the inverse CDF method: -ln(U) where U ~ Uniform(0,1).
22    ///
23    /// # Arguments
24    /// * `size` - Number of values to generate.
25    ///
26    /// # Errors
27    /// Returns `FerrayError::InvalidValue` if `size` is zero.
28    pub fn standard_exponential(&mut self, size: usize) -> Result<Array<f64, Ix1>, FerrayError> {
29        if size == 0 {
30            return Err(FerrayError::invalid_value("size must be > 0"));
31        }
32        let data = generate_vec(self, size, standard_exponential_single);
33        vec_to_array1(data)
34    }
35
36    /// Generate an array of exponential variates with the given scale.
37    ///
38    /// The exponential distribution has PDF: f(x) = (1/scale) * exp(-x/scale).
39    ///
40    /// # Arguments
41    /// * `scale` - Scale parameter (1/rate), must be positive.
42    /// * `size` - Number of values to generate.
43    ///
44    /// # Errors
45    /// Returns `FerrayError::InvalidValue` if `scale <= 0` or `size` is zero.
46    pub fn exponential(&mut self, scale: f64, size: usize) -> Result<Array<f64, Ix1>, FerrayError> {
47        if size == 0 {
48            return Err(FerrayError::invalid_value("size must be > 0"));
49        }
50        if scale <= 0.0 {
51            return Err(FerrayError::invalid_value(format!(
52                "scale must be positive, got {scale}"
53            )));
54        }
55        let data = generate_vec(self, size, |bg| scale * standard_exponential_single(bg));
56        vec_to_array1(data)
57    }
58}
59
60#[cfg(test)]
61mod tests {
62    use crate::default_rng_seeded;
63
64    #[test]
65    fn standard_exponential_positive() {
66        let mut rng = default_rng_seeded(42);
67        let arr = rng.standard_exponential(10_000).unwrap();
68        let slice = arr.as_slice().unwrap();
69        for &v in slice {
70            assert!(
71                v > 0.0,
72                "standard_exponential produced non-positive value: {v}"
73            );
74        }
75    }
76
77    #[test]
78    fn standard_exponential_mean_variance() {
79        let mut rng = default_rng_seeded(42);
80        let n = 100_000;
81        let arr = rng.standard_exponential(n).unwrap();
82        let slice = arr.as_slice().unwrap();
83        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
84        let var: f64 = slice.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n as f64;
85        // Exp(1): mean=1, var=1
86        let se = (1.0 / n as f64).sqrt();
87        assert!((mean - 1.0).abs() < 3.0 * se, "mean {mean} too far from 1");
88        assert!((var - 1.0).abs() < 0.05, "variance {var} too far from 1");
89    }
90
91    #[test]
92    fn exponential_mean() {
93        let mut rng = default_rng_seeded(42);
94        let n = 100_000;
95        let scale = 3.0;
96        let arr = rng.exponential(scale, n).unwrap();
97        let slice = arr.as_slice().unwrap();
98        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
99        let se = (scale * scale / n as f64).sqrt();
100        assert!(
101            (mean - scale).abs() < 3.0 * se,
102            "mean {mean} too far from {scale}"
103        );
104    }
105
106    #[test]
107    fn exponential_bad_scale() {
108        let mut rng = default_rng_seeded(42);
109        assert!(rng.exponential(0.0, 100).is_err());
110        assert!(rng.exponential(-1.0, 100).is_err());
111    }
112
113    #[test]
114    fn exponential_deterministic() {
115        let mut rng1 = default_rng_seeded(42);
116        let mut rng2 = default_rng_seeded(42);
117        let a = rng1.exponential(2.0, 100).unwrap();
118        let b = rng2.exponential(2.0, 100).unwrap();
119        assert_eq!(a.as_slice().unwrap(), b.as_slice().unwrap());
120    }
121}
ferray_random/distributions/exponential.rs

ferray_random/distributions/
exponential.rs
