stats_claw/distributions/positive/
beta.rs

1//! Beta distribution numerics, for the [`BetaDistribution`] parameter struct.
2//!
3//! Equivalent to `scipy.stats.beta(alpha_parameter, beta_parameter)` on the
4//! support `[0, 1]`. The density uses `ln_beta` for stability, the CDF is the
5//! regularized incomplete beta `betai`, the quantile inverts it by bracketed
6//! bisection, and sampling forms `X / (X + Y)` from two unit-scale gamma draws.
7//!
8//! # Examples
9//!
10//! ```
11//! use stats_claw::distributions::{Cdf, Pdf};
12//! use stats_claw::distributions::BetaDistribution;
13//!
14//! // Beta(1, 1) is the standard uniform.
15//! let d = BetaDistribution { alpha_parameter: 1.0, beta_parameter: 1.0, ..Default::default() };
16//! assert!((d.pdf(0.3) - 1.0).abs() < 1e-12, "pdf was {}", d.pdf(0.3));
17//! assert!((d.cdf(0.4) - 0.4).abs() < 1e-12);
18//! ```
19
20use super::super::{Cdf, Moments, Pdf, Quantile, Sample, bisection_quantile};
21use super::gamma::marsaglia_tsang;
22use crate::distributions::BetaDistribution;
23use crate::rng::SplitMix64;
24use crate::special::{betai, ln_beta};
25
26impl Pdf for BetaDistribution {
27    fn pdf(&self, x: f64) -> f64 {
28        if !(0.0..=1.0).contains(&x) {
29            return 0.0;
30        }
31        let (a, b) = (self.alpha_parameter, self.beta_parameter);
32        let ln_density = (a - 1.0).mul_add(x.ln(), (b - 1.0) * (1.0 - x).ln()) - ln_beta(a, b);
33        ln_density.exp()
34    }
35}
36
37impl Cdf for BetaDistribution {
38    fn cdf(&self, x: f64) -> f64 {
39        betai(self.alpha_parameter, self.beta_parameter, x)
40    }
41}
42
43impl Quantile for BetaDistribution {
44    fn quantile(&self, p: f64) -> f64 {
45        bisection_quantile(p, 0.0, 1.0, |x| self.cdf(x))
46    }
47}
48
49impl Moments for BetaDistribution {
50    fn mean(&self) -> Option<f64> {
51        let (a, b) = (self.alpha_parameter, self.beta_parameter);
52        Some(a / (a + b))
53    }
54    fn variance(&self) -> Option<f64> {
55        let (a, b) = (self.alpha_parameter, self.beta_parameter);
56        let s = a + b;
57        let denom = s * s * (s + 1.0);
58        Some((a * b) / denom)
59    }
60}
61
62impl Sample for BetaDistribution {
63    fn sample(&self, rng: &mut SplitMix64) -> f64 {
64        let x = marsaglia_tsang(self.alpha_parameter, rng);
65        let y = marsaglia_tsang(self.beta_parameter, rng);
66        x / (x + y)
67    }
68}
69
70#[cfg(test)]
71mod tests {
72    use super::*;
73
74    /// Beta(1, 1) is the standard uniform: density 1 everywhere on [0, 1].
75    #[test]
76    fn beta_one_one_is_uniform() {
77        let d = BetaDistribution {
78            alpha_parameter: 1.0,
79            beta_parameter: 1.0,
80            ..Default::default()
81        };
82        assert!((d.pdf(0.3) - 1.0).abs() < 1e-12, "was {}", d.pdf(0.3));
83        assert!((d.cdf(0.4) - 0.4).abs() < 1e-12);
84    }
85
86    /// The mean is `a / (a + b)`.
87    #[test]
88    fn mean_is_a_over_a_plus_b() {
89        let d = BetaDistribution {
90            alpha_parameter: 2.0,
91            beta_parameter: 5.0,
92            ..Default::default()
93        };
94        assert!(matches!(d.mean(), Some(m) if (m - 2.0 / 7.0).abs() < 1e-12));
95    }
96}
stats_claw/distributions/positive/beta.rs

stats_claw/distributions/positive/
beta.rs
