numeris/stats/
beta_dist.rs

1use crate::FloatScalar;
2use crate::special::{lbeta, betainc};
3use super::{ContinuousDistribution, StatsError, quantile_newton};
4
5/// Beta distribution with shape parameters α and β on [0, 1].
6///
7/// f(x) = x^{α−1} (1−x)^{β−1} / B(α, β) for 0 ≤ x ≤ 1.
8///
9/// # Example
10///
11/// ```
12/// use numeris::stats::{Beta, ContinuousDistribution};
13///
14/// let b = Beta::new(2.0_f64, 5.0).unwrap();
15/// assert!((b.mean() - 2.0/7.0).abs() < 1e-14);
16/// ```
17#[derive(Debug, Clone, Copy)]
18pub struct Beta<T> {
19    alpha: T,
20    beta: T,
21}
22
23impl<T: FloatScalar> Beta<T> {
24    /// Create a Beta distribution with shape parameters `alpha` and `beta`.
25    /// Requires both > 0.
26    pub fn new(alpha: T, beta: T) -> Result<Self, StatsError> {
27        if alpha <= T::zero() || beta <= T::zero() {
28            return Err(StatsError::InvalidParameter);
29        }
30        Ok(Self { alpha, beta })
31    }
32}
33
34impl<T: FloatScalar> Beta<T> {
35    /// Draw a random sample from this distribution.
36    ///
37    /// Samples X ~ Gamma(alpha, 1), Y ~ Gamma(beta, 1), returns X / (X + Y).
38    pub fn sample(&self, rng: &mut super::Rng) -> T {
39        let x = rng.next_gamma(self.alpha);
40        let y = rng.next_gamma(self.beta);
41        x / (x + y)
42    }
43
44    /// Fill a fixed-size array with independent samples.
45    pub fn sample_array<const K: usize>(&self, rng: &mut super::Rng) -> [T; K] {
46        let mut out = [T::zero(); K];
47        for v in out.iter_mut() {
48            *v = self.sample(rng);
49        }
50        out
51    }
52}
53
54impl<T: FloatScalar> ContinuousDistribution<T> for Beta<T> {
55    fn pdf(&self, x: T) -> T {
56        if x < T::zero() || x > T::one() {
57            return T::zero();
58        }
59        self.ln_pdf(x).exp()
60    }
61
62    fn ln_pdf(&self, x: T) -> T {
63        if x < T::zero() || x > T::one() {
64            return T::neg_infinity();
65        }
66        let one = T::one();
67        (self.alpha - one) * x.ln() + (self.beta - one) * (one - x).ln()
68            - lbeta(self.alpha, self.beta)
69    }
70
71    fn cdf(&self, x: T) -> T {
72        if x <= T::zero() {
73            return T::zero();
74        }
75        if x >= T::one() {
76            return T::one();
77        }
78        betainc(self.alpha, self.beta, x).unwrap_or(T::nan())
79    }
80
81    fn quantile(&self, p: T) -> T {
82        let x0 = self.mean();
83        let eps = T::epsilon();
84        quantile_newton(
85            |x| self.cdf(x),
86            |x| self.pdf(x),
87            p,
88            x0,
89            eps,
90            T::one() - eps,
91        )
92    }
93
94    fn mean(&self) -> T {
95        self.alpha / (self.alpha + self.beta)
96    }
97
98    fn variance(&self) -> T {
99        let ab = self.alpha + self.beta;
100        self.alpha * self.beta / (ab * ab * (ab + T::one()))
101    }
102}
numeris/stats/beta_dist.rs

numeris/stats/
beta_dist.rs
