qfall_math/rational/q/sample/
gauss.rs

1// Copyright © 2024 Marvin Beckmann
2//
3// This file is part of qFALL-math.
4//
5// qFALL-math is free software: you can redistribute it and/or modify it under
6// the terms of the Mozilla Public License Version 2.0 as published by the
7// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
8
9//! This module contains sampling algorithms for gaussian distributions over [`Q`].
10
11use crate::{error::MathError, rational::Q};
12use probability::{
13    distribution::{Gaussian, Sample},
14    source,
15};
16use rand::RngCore;
17
18impl Q {
19    /// Chooses a [`Q`] instance according to the continuous Gaussian distribution.
20    ///
21    /// Parameters:
22    /// - `center`: specifies the position of the center
23    /// - `sigma`: specifies the standard deviation
24    ///
25    /// Returns new [`Q`] sample chosen according to the specified continuous Gaussian
26    /// distribution or a [`MathError`] if the specified parameters were not chosen
27    /// appropriately (`sigma > 0`).
28    ///
29    /// # Examples
30    /// ```
31    /// use qfall_math::rational::Q;
32    ///
33    /// let sample = Q::sample_gauss(0, 1).unwrap();
34    /// ```
35    ///
36    /// # Errors and Failures
37    /// - Returns a [`MathError`] of type [`NonPositive`](MathError::NonPositive)
38    ///   if `sigma <= 0`.
39    pub fn sample_gauss(center: impl Into<Q>, sigma: impl Into<f64>) -> Result<Q, MathError> {
40        let center = center.into();
41        let sigma = sigma.into();
42        if sigma <= 0.0 {
43            return Err(MathError::NonPositive(format!(
44                "The sigma has to be positive and not zero, but the provided value is {sigma}."
45            )));
46        }
47        let mut rng = rand::rng();
48        let mut source = source::default(rng.next_u64());
49
50        // Instead of sampling with a center of c, we sample with center 0 and add the
51        // center later. These are equivalent and this way we can sample in larger ranges
52        let sampler = Gaussian::new(0.0, sigma);
53        let sample = center + Q::from(sampler.sample(&mut source));
54
55        Ok(sample)
56    }
57}
58
59#[cfg(test)]
60mod test_sample_gauss {
61    use crate::rational::Q;
62
63    /// Test correct distribution with a confidence level of 99.7% -> 3 standard
64    /// deviations.
65    #[test]
66    fn in_concentration_bound() {
67        let range = 3;
68        for (mu, sigma) in [(i64::MAX, 1), (0, 20), (i64::MIN, 100)] {
69            assert!(range * sigma >= (Q::from(mu) - Q::sample_gauss(mu, sigma).unwrap()).abs())
70        }
71    }
72
73    /// Ensure that an error is returned if `sigma` is not positive
74    #[test]
75    fn non_positive_sigma() {
76        for (mu, sigma) in [(0, 0), (0, -1)] {
77            assert!(Q::sample_gauss(mu, sigma).is_err())
78        }
79    }
80}
qfall_math/rational/q/sample/gauss.rs

qfall_math/rational/q/sample/
gauss.rs
