prodef 0.2.2

A simple Rust crate for handling probability distributions, primarily intended for use with Bayesian inference.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

use crate::{Density, SamplingMode, domain::Domain, macros::tval};
use nalgebra::{Dim, OVector, RealField, SVector, Scalar, U1, VectorView};
use rand::RngExt;
use rand_distr::{Uniform, uniform::SampleUniform};
use serde::{Deserialize, Serialize};

/// A loguniform PDF.
#[derive(Clone, Debug, Deserialize, PartialEq, Serialize)]
pub struct LogUniformDensity<T>(Domain<T, U1>)
where
    T: Scalar;

impl<T> LogUniformDensity<T>
where
    T: RealField,
{
    /// Create a new [`LogUniformDensity`].
    ///
    /// Returns `None` if:
    /// - `a >= b` (inverted or degenerate bounds)
    /// - `a <= 0` or `b <= 0` (log domain requires positive values)
    pub fn new(a: T, b: T) -> Option<Self> {
        if a >= b || a <= T::zero() || b <= T::zero() {
            None
        } else {
            Some(Self(Domain::new_mdomain(OVector::from_element_generic(
                U1,
                U1,
                (Some(a), Some(b)),
            ))))
        }
    }

    /// Returns the maximum value of the domain.
    pub fn maximum(&self) -> T {
        match &self.0.inner().unwrap() {
            (_, Some(max)) => max.clone(),
            // Safe: LogUniformDensity constructor validates a < b condition
            // and creates MDomain with explicit bounds, guaranteeing (min, max) both Some.
            _ => unreachable!("MDomain always has explicit bounds in LogUniformDensity"),
        }
    }

    /// Returns the minimum value of the domain.
    pub fn minimum(&self) -> T {
        match &self.0.inner().unwrap() {
            (Some(min), _) => min.clone(),
            // Safe: LogUniformDensity constructor validates a < b condition
            // and creates MDomain with explicit bounds, guaranteeing (min, max) both Some.
            _ => unreachable!("MDomain always has explicit bounds in LogUniformDensity"),
        }
    }
}

impl<T> Density<T, U1> for LogUniformDensity<T>
where
    T: RealField + SampleUniform,
{
    fn density<RStride: Dim, CStride: Dim>(
        &self,
        sample: &VectorView<T, U1, RStride, CStride>,
    ) -> Option<T> {
        if !self.0.contains(sample) {
            return None;
        }

        Some(
            T::one()
                / (sample[0].clone()
                    * (self.0.maximum_values()[0].clone().unwrap().ln()
                        - self.0.minimum_values()[0].clone().unwrap().ln())),
        )
    }

    fn domain(&self) -> Domain<T, U1> {
        self.0.clone()
    }

    fn mean(&self) -> SVector<T, 1> {
        let a = self.minimum();
        let b = self.maximum();
        SVector::from([(b.clone() - a.clone()) / (b.ln() - a.ln())])
    }

    fn sample(&self, rng: &mut impl RngExt, _mode: &SamplingMode) -> Option<SVector<T, 1>> {
        let uniform = Uniform::new_inclusive(
            self.0.minimum_values()[0].clone().unwrap().ln(),
            self.0.maximum_values()[0].clone().unwrap().ln(),
        )
        .unwrap();

        Some(SVector::from([rng.sample(uniform).exp()]))
    }

    fn sample_iter(&self, rng: &mut impl RngExt) -> impl Iterator<Item = Option<SVector<T, 1>>> {
        let uniform = Uniform::new_inclusive(
            self.0.minimum_values()[0].clone().unwrap().ln(),
            self.0.maximum_values()[0].clone().unwrap().ln(),
        )
        .unwrap();

        rng.sample_iter(uniform)
            .map(|value| Some(OVector::from_element_generic(U1, U1, value.exp())))
    }

    fn variance(&self) -> SVector<T, 1> {
        // For a log-uniform distribution on [a, b], the variance is:
        // variance = E[X²] - (E[X])²
        // where E[X²] = (b² - a²) / (2 * (ln(b) - ln(a)))
        let a = self.minimum();
        let b = self.maximum();
        let ln_ratio = b.clone().ln() - a.clone().ln();

        let a_sq = a.clone() * a;
        let b_sq = b.clone() * b;
        let e_x_sq = (b_sq - a_sq) / (tval!(2.0, f64) * ln_ratio);

        let mean = self.mean()[0].clone();
        let mean_sq = mean.clone() * mean;

        SVector::from([e_x_sq - mean_sq])
    }
}

impl<T> TryFrom<crate::univariate::UnivariateDensity<T>> for LogUniformDensity<T>
where
    T: RealField,
{
    type Error = ();

    fn try_from(value: crate::univariate::UnivariateDensity<T>) -> Result<Self, Self::Error> {
        match value {
            crate::univariate::UnivariateDensity::Loguniform(pdf) => Ok(pdf),
            _ => Err(()),
        }
    }
}