probability_rs/dist/
normal.rs

1use crate::dist::{Continuous, DistError, Distribution, Moments};
2use crate::{num, rng::RngCore};
3
4#[derive(Debug, Clone, Copy)]
5pub struct Normal {
6    mu: f64,
7    sigma: f64,
8    inv_sigma: f64,
9    norm: f64, // 1/(sigma*sqrt(2*pi))
10}
11
12impl Normal {
13    pub fn new(mu: f64, sigma: f64) -> Result<Self, DistError> {
14        if !(sigma > 0.0 && sigma.is_finite() && mu.is_finite()) {
15            return Err(DistError::InvalidParameter);
16        }
17        let inv_sigma = 1.0 / sigma;
18        let norm = num::INV_SQRT_2PI * inv_sigma;
19        Ok(Self {
20            mu,
21            sigma,
22            inv_sigma,
23            norm,
24        })
25    }
26    #[inline]
27    pub fn mean_param(&self) -> f64 {
28        self.mu
29    }
30    #[inline]
31    pub fn sigma(&self) -> f64 {
32        self.sigma
33    }
34}
35
36impl Distribution for Normal {
37    type Value = f64;
38    fn cdf(&self, x: f64) -> f64 {
39        let z = (x - self.mu) * self.inv_sigma;
40        num::standard_normal_cdf(z)
41    }
42    fn in_support(&self, x: f64) -> bool {
43        x.is_finite()
44    }
45    fn sample<R: RngCore>(&self, rng: &mut R) -> f64 {
46        // Box-Muller polar (Marsaglia) without external dependencies.
47        loop {
48            let u1 = 2.0 * rng.next_f64() - 1.0; // (-1,1)
49            let u2 = 2.0 * rng.next_f64() - 1.0;
50            let s = u1 * u1 + u2 * u2;
51            if s >= 1.0 || s == 0.0 {
52                continue;
53            }
54            let factor = (-2.0 * s.ln() / s).sqrt();
55            return self.mu + self.sigma * u1 * factor;
56        }
57    }
58}
59
60impl Continuous for Normal {
61    fn pdf(&self, x: f64) -> f64 {
62        if !self.in_support(x) {
63            return 0.0;
64        }
65        let z = (x - self.mu) * self.inv_sigma;
66        self.norm * (-0.5 * z * z).exp()
67    }
68    fn inv_cdf(&self, p: f64) -> f64 {
69        self.mu + self.sigma * num::standard_normal_inv_cdf(p)
70    }
71}
72
73impl Moments for Normal {
74    fn mean(&self) -> f64 {
75        self.mu
76    }
77    fn variance(&self) -> f64 {
78        self.sigma * self.sigma
79    }
80    fn skewness(&self) -> f64 { 0.0 }
81    fn kurtosis(&self) -> f64 { 0.0 }
82    fn entropy(&self) -> f64 {
83        0.5 * (2.0 * std::f64::consts::PI * std::f64::consts::E * self.sigma * self.sigma).ln()
84    }
85}
86
87#[cfg(test)]
88mod tests {
89    use super::*;
90    use crate::rng::SplitMix64;
91
92    #[test]
93    fn normal_basic() {
94        let n = Normal::new(0.0, 1.0).unwrap();
95        assert!((n.pdf(0.0) - 0.3989422804014327).abs() < 1e-12);
96        // CDF approximation via erf has typical error ~1e-7; use generous tolerance.
97        assert!((n.cdf(0.0) - 0.5).abs() < 2e-6);
98        let q = n.inv_cdf(0.975);
99        assert!((q - 1.959963).abs() < 5e-4);
100    }
101
102    #[test]
103    fn cdf_symmetry() {
104        let n = Normal::new(0.0, 1.0).unwrap();
105        let z = 0.7;
106        let f = n.cdf(z);
107        let f_neg = n.cdf(-z);
108        assert!((f + f_neg - 1.0).abs() < 3e-15);
109    }
110
111    #[test]
112    fn moments_higher() {
113        let n = Normal::new(0.0, 2.0).unwrap();
114        assert_eq!(n.skewness(), 0.0);
115        assert_eq!(n.kurtosis(), 0.0);
116        assert_eq!(n.kurtosis_full(), 3.0);
117    }
118
119    #[test]
120    fn sampling_determinism() {
121        let n = Normal::new(0.0, 1.0).unwrap();
122        let mut rng1 = SplitMix64::seed_from_u64(42);
123        let mut rng2 = SplitMix64::seed_from_u64(42);
124        let x1 = n.sample(&mut rng1);
125        let x2 = n.sample(&mut rng2);
126        assert_eq!(x1.to_bits(), x2.to_bits());
127    }
128
129    #[test]
130    fn entropy_normal() {
131        let n = Normal::new(0.0, 2.0).unwrap();
132        let expected = 0.5 * (2.0 * std::f64::consts::PI * std::f64::consts::E * 4.0).ln();
133        assert!((n.entropy() - expected).abs() < 1e-12);
134    }
135}
probability_rs/dist/normal.rs

probability_rs/dist/
normal.rs
