ringkernel_montecarlo/variance/
antithetic.rs

1//! Antithetic variates for variance reduction.
2//!
3//! Antithetic variates use negatively correlated pairs of samples to reduce
4//! variance. If U is uniform on [0,1], then 1-U is also uniform on [0,1]
5//! but negatively correlated with U.
6
7use crate::rng::GpuRng;
8
9/// Generate an antithetic pair of uniform variates.
10///
11/// Returns (u, 1-u) where u is uniform in [0, 1).
12/// The pair is negatively correlated, which can reduce variance
13/// when the estimator is monotonic.
14///
15/// # Example
16///
17/// ```ignore
18/// let (u1, u2) = antithetic_pair(&mut state);
19/// // Estimate E[f(U)] using both samples
20/// let estimate = 0.5 * (f(u1) + f(u2));
21/// ```
22pub fn antithetic_pair<R: GpuRng>(state: &mut R::State) -> (f32, f32) {
23    let u = R::next_uniform(state);
24    (u, 1.0 - u)
25}
26
27/// Generate an antithetic pair of normal variates.
28///
29/// Returns (z, -z) where z is standard normal.
30pub fn antithetic_normal_pair<R: GpuRng>(state: &mut R::State) -> (f32, f32) {
31    let z = R::next_normal(state);
32    (z, -z)
33}
34
35/// Configuration for antithetic variates estimator.
36#[derive(Debug, Clone)]
37pub struct AntitheticVariates {
38    /// Number of antithetic pairs to generate.
39    pub n_pairs: usize,
40}
41
42impl AntitheticVariates {
43    /// Create new antithetic variates configuration.
44    pub fn new(n_pairs: usize) -> Self {
45        Self { n_pairs }
46    }
47
48    /// Estimate E[f(U)] using antithetic variates.
49    ///
50    /// The estimator is:
51    /// `(1/n) * sum_i (f(U_i) + f(1-U_i)) / 2`
52    ///
53    /// This has lower variance than naive estimation when f is monotonic.
54    pub fn estimate<R: GpuRng, F>(&self, state: &mut R::State, f: F) -> f32
55    where
56        F: Fn(f32) -> f32,
57    {
58        let mut sum = 0.0;
59        for _ in 0..self.n_pairs {
60            let (u1, u2) = antithetic_pair::<R>(state);
61            sum += 0.5 * (f(u1) + f(u2));
62        }
63        sum / self.n_pairs as f32
64    }
65
66    /// Estimate E[f(Z)] where Z ~ N(0,1) using antithetic variates.
67    pub fn estimate_normal<R: GpuRng, F>(&self, state: &mut R::State, f: F) -> f32
68    where
69        F: Fn(f32) -> f32,
70    {
71        let mut sum = 0.0;
72        for _ in 0..self.n_pairs {
73            let (z1, z2) = antithetic_normal_pair::<R>(state);
74            sum += 0.5 * (f(z1) + f(z2));
75        }
76        sum / self.n_pairs as f32
77    }
78}
79
80impl Default for AntitheticVariates {
81    fn default() -> Self {
82        Self::new(1000)
83    }
84}
85
86#[cfg(test)]
87mod tests {
88    use super::*;
89    use crate::rng::PhiloxRng;
90
91    #[test]
92    fn test_antithetic_pair_range() {
93        let mut state = PhiloxRng::seed(42, 0);
94
95        for _ in 0..100 {
96            let (u1, u2) = antithetic_pair::<PhiloxRng>(&mut state);
97            assert!((0.0..1.0).contains(&u1));
98            assert!((0.0..1.0).contains(&u2));
99            assert!((u1 + u2 - 1.0).abs() < 1e-6, "u1 + u2 should equal 1");
100        }
101    }
102
103    #[test]
104    fn test_antithetic_normal_pair() {
105        let mut state = PhiloxRng::seed(42, 0);
106
107        for _ in 0..100 {
108            let (z1, z2) = antithetic_normal_pair::<PhiloxRng>(&mut state);
109            assert!((z1 + z2).abs() < 1e-6, "z1 + z2 should equal 0");
110        }
111    }
112
113    #[test]
114    fn test_antithetic_variance_reduction() {
115        // Test on f(x) = x^2, estimating E[X^2] where X ~ U[0,1]
116        // True value is 1/3
117        let n = 10000;
118
119        // Naive estimation
120        let mut state = PhiloxRng::seed(123, 0);
121        let naive_sum: f32 = (0..n)
122            .map(|_| {
123                let u = PhiloxRng::next_uniform(&mut state);
124                u * u
125            })
126            .sum();
127        let naive_estimate = naive_sum / n as f32;
128
129        // Antithetic estimation
130        let mut state = PhiloxRng::seed(123, 0);
131        let av = AntitheticVariates::new(n / 2);
132        let av_estimate = av.estimate::<PhiloxRng, _>(&mut state, |u| u * u);
133
134        let true_value = 1.0 / 3.0;
135
136        // Both should be close to true value
137        assert!(
138            (naive_estimate - true_value).abs() < 0.05,
139            "Naive estimate {} far from {}",
140            naive_estimate,
141            true_value
142        );
143        assert!(
144            (av_estimate - true_value).abs() < 0.05,
145            "AV estimate {} far from {}",
146            av_estimate,
147            true_value
148        );
149    }
150}
ringkernel_montecarlo/variance/antithetic.rs

ringkernel_montecarlo/variance/
antithetic.rs
