Am/test/hl api noise checks

tfhe/src/core_crypto/algorithms/test/mod.rs

@@ -24,7 +24,7 @@ mod lwe_packing_keyswitch_key_generation;
 mod lwe_private_functional_packing_keyswitch;
 pub(crate) mod lwe_programmable_bootstrapping;
 mod modulus_switch_compression;
-mod noise_distribution;
+pub(crate) mod noise_distribution;
 
 pub struct TestResources {
     pub seeder: Box<dyn Seeder>,

tfhe/src/core_crypto/algorithms/test/noise_distribution/lwe_encryption_noise.rs

@@ -88,7 +88,7 @@ create_parameterized_test!(lwe_encrypt_decrypt_noise_distribution_custom_mod {
     TEST_PARAMS_3_BITS_63_U64
 });
 
-fn lwe_compact_public_key_encryption_expected_variance(
+pub(crate) fn lwe_compact_public_key_encryption_expected_variance(
     input_noise: impl DispersionParameter,
     lwe_dimension: LweDimension,
 ) -> Variance {

tfhe/src/core_crypto/algorithms/test/noise_distribution/mod.rs

@@ -1,6 +1,6 @@
 use super::*;
 
-mod lwe_encryption_noise;
+pub(crate) mod lwe_encryption_noise;
 mod lwe_keyswitch_noise;
 // We are having crashes on aarch64 at the moment, problem is the code paths are not the same
 // between archs, so we disable those on the Apple M1

tfhe/src/core_crypto/commons/math/random/t_uniform.rs

@@ -70,6 +70,22 @@ impl<T: UnsignedInteger> TUniform<T> {
     pub fn max_value_inclusive(&self) -> T::Signed {
         T::Signed::ONE << self.bound_log2 as usize
     }
+
+    // TODO: return the ModularVariance? Essentially the first part of the formula corresponds to
+    // what older iterations of the lib called ModularVariance, this struct does not currently have
+    // the DispersionParameter trait implemented, so here for convenience we compute the Variance
+    // directly by dividing the ModularVariance by the modulus squared to be in the Torus domain.
+    #[cfg(test)]
+    #[allow(clippy::trivially_copy_pass_by_ref)]
+    pub(crate) fn variance(
+        &self,
+        modulus: f64,
+    ) -> crate::core_crypto::commons::dispersion::Variance {
+        // (2^{2 * b + 1} + 1) / 6
+        crate::core_crypto::commons::dispersion::Variance(
+            ((2.0f64.powi(2 * (self.bound_log2 as i32) + 1) + 1.0f64) / 6.0f64) * modulus.powi(-2),
+        )
+    }
 }
 
 macro_rules! implement_t_uniform_uint {

tfhe/src/core_crypto/commons/mod.rs

@@ -26,12 +26,13 @@ pub mod traits;
 #[cfg(test)]
 pub mod test_tools {
     use rand::Rng;
+    use statrs::distribution::{ChiSquared, ContinuousCDF};
 
     pub use crate::core_crypto::algorithms::misc::{
         modular_distance, modular_distance_custom_mod, torus_modular_diff,
     };
     use crate::core_crypto::commons::ciphertext_modulus::CiphertextModulus;
-    use crate::core_crypto::commons::dispersion::{DispersionParameter, Variance};
+    use crate::core_crypto::commons::dispersion::{DispersionParameter, StandardDev, Variance};
     use crate::core_crypto::commons::generators::{
         EncryptionRandomGenerator, SecretRandomGenerator,
     };
@@ -45,14 +46,76 @@ pub mod test_tools {
     use crate::core_crypto::commons::traits::*;
     use tfhe_csprng::seeders::Seed;
 
+    pub fn arithmetic_mean(samples: &[f64]) -> f64 {
+        samples.iter().copied().sum::<f64>() / samples.len() as f64
+    }
+
     pub fn variance(samples: &[f64]) -> Variance {
         let num_samples = samples.len();
-        let mean = samples.iter().sum::<f64>() / (num_samples as f64);
+        let mean = arithmetic_mean(samples);
         Variance(
             samples.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / ((num_samples - 1) as f64),
         )
     }
 
+    #[derive(Debug, Clone, Copy)]
+    pub struct VarianceConfidenceInterval {
+        lower_bound: Variance,
+        upper_bound: Variance,
+    }
+
+    impl VarianceConfidenceInterval {
+        pub fn variance_is_in_interval(&self, variance_to_check: Variance) -> bool {
+            self.lower_bound <= variance_to_check && self.upper_bound >= variance_to_check
+        }
+
+        pub fn lower_bound(&self) -> Variance {
+            self.lower_bound
+        }
+
+        pub fn upper_bound(&self) -> Variance {
+            self.upper_bound
+        }
+    }
+
+    /// Can apply to any distribution as long as the samples are idependant and identically
+    /// distributed.
+    ///
+    /// As this construction relies on the Central Limit Theorem you must ensure that you have
+    /// enough samples for the Central Limit Theorem to apply. Some accepted values can be as low as
+    /// 30 or 50, though if you can afford it, to have tighter intervals and to have a better
+    /// convergence for the CLT you should try to gather more samples.
+    pub fn any_distribution_variance_confidence_interval() {
+        todo!()
+    }
+
+    /// Note that the confidence_level is the probability to be in the interval to be computed.
+    /// Applies only to samples known to come from a normal distribution
+    pub fn normal_distribution_variance_confidence_interval(
+        sample_count: f64,
+        measured_variance: Variance,
+        confidence_level: f64,
+    ) -> VarianceConfidenceInterval {
+        assert!(confidence_level >= 0.0);
+        assert!(confidence_level <= 1.0);
+        let alpha = 1.0 - confidence_level;
+        let degrees_of_freedom = sample_count - 1.0;
+        let chi2 = ChiSquared::new(degrees_of_freedom).unwrap();
+        let chi2_lower = chi2.inverse_cdf(alpha / 2.0);
+        let chi2_upper = chi2.inverse_cdf(1.0 - alpha / 2.0);
+
+        // Lower bound is divided by chi2_upper, upper bound divided by chi2_lower
+        let lower_bound = Variance(degrees_of_freedom * measured_variance.0 / chi2_upper);
+        let upper_bound = Variance(degrees_of_freedom * measured_variance.0 / chi2_lower);
+
+        assert!(lower_bound <= upper_bound);
+
+        VarianceConfidenceInterval {
+            lower_bound,
+            upper_bound,
+        }
+    }
+
     pub fn new_random_generator() -> RandomGenerator<DefaultRandomGenerator> {
         RandomGenerator::new(random_seed())
     }
@@ -346,4 +409,48 @@ pub mod test_tools {
             assert_eq!(distance, -0.25);
         }
     }
+
+    #[test]
+    fn test_confidence_interval() {
+        // https://stats.libretexts.org/Bookshelves/Introductory_Statistics/
+        // Inferential_Statistics_and_Probability_-_A_Holistic_Approach_(Geraghty)/
+        // 09%3A_Point_Estimation_and_Confidence_Intervals/9.03%3A_Confidence_Intervals
+
+        // In performance measurement of investments, standard deviation is a measure of volatility
+        // or risk. Twenty monthly returns from a mutual fund show an average monthly return of
+        // 1 percent and a sample standard deviation of 5 percent.
+        // Find a 95% confidence interval for the monthly standard deviation of the mutual fund.
+
+        // The Chi‐square distribution will have 20‐1 =19 degrees of freedom. Using technology,
+        // we find that the two critical values are  chi2_left=8.90655
+        // and   chi2_right=32.8523
+        // Formula for confidence interval for sigma
+        // is:  sqrt(19 * 5^2 / 32.8523) sqrt(19 * 5^2 / 8.90655) = (3.8,7.3)
+
+        // One can say with 95% confidence that the standard deviation for this mutual fund is
+        // between 3.8 and 7.3 percent per month.
+
+        let measured_std_dev = StandardDev(0.05);
+        let measured_variance = measured_std_dev.get_variance();
+
+        let confidence_level = 0.95;
+
+        let confidence_interval = normal_distribution_variance_confidence_interval(
+            20.,
+            measured_variance,
+            confidence_level,
+        );
+
+        let lower_bound = confidence_interval.lower_bound();
+        let upper_bound = confidence_interval.upper_bound();
+
+        let approx_expected_lower_bound = StandardDev(0.038).get_variance();
+        let approx_expected_upper_bound = StandardDev(0.073).get_variance();
+
+        let lower_bound_abs_diff = (lower_bound.0 - approx_expected_lower_bound.0).abs();
+        let upper_bound_abs_diff = (upper_bound.0 - approx_expected_upper_bound.0).abs();
+
+        assert!(lower_bound_abs_diff / approx_expected_lower_bound.0 < 0.01);
+        assert!(upper_bound_abs_diff / approx_expected_upper_bound.0 < 0.01);
+    }
 }

tfhe/src/high_level_api/tests/mod.rs

@@ -1,3 +1,4 @@
+mod noise_distribution;
 mod tags_on_entities;
 
 use crate::high_level_api::prelude::*;