use-monte-carlo 0.0.1

Primitive Monte Carlo helpers for repeatable simulation runs
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

#![forbid(unsafe_code)]
//! Primitive Monte Carlo helpers for repeatable simulation runs.
//!
//! The crate evaluates a caller-provided observation function over deterministic
//! unit samples derived from a `SimulationSeed`.
//!
//! # Examples
//!
//! ```rust
//! use use_monte_carlo::monte_carlo_mean;
//! use use_seed::SimulationSeed;
//!
//! let estimate = monte_carlo_mean(SimulationSeed::new(4), 8, |sample| sample * 2.0).unwrap();
//!
//! assert_eq!(estimate.samples, 8);
//! assert!(estimate.mean.is_finite());
//! assert!(estimate.variance.is_finite());
//! ```

use use_seed::SimulationSeed;

#[derive(Debug, Clone, Copy, PartialEq)]
pub struct MonteCarloEstimate {
    pub mean: f64,
    pub variance: f64,
    pub samples: usize,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum MonteCarloError {
    InvalidSampleCount,
    NonFiniteObservation,
}

pub fn monte_carlo_mean<F>(
    seed: SimulationSeed,
    samples: usize,
    observation: F,
) -> Result<MonteCarloEstimate, MonteCarloError>
where
    F: Fn(f64) -> f64,
{
    if samples == 0 {
        return Err(MonteCarloError::InvalidSampleCount);
    }

    let mut values = Vec::with_capacity(samples);
    for index in 0..samples {
        let sample = seed.mix(index + 1).to_unit_f64();
        let value = observation(sample);
        if !value.is_finite() {
            return Err(MonteCarloError::NonFiniteObservation);
        }

        values.push(value);
    }

    let mean = values.iter().sum::<f64>() / samples as f64;
    let variance = values
        .iter()
        .map(|value| {
            let diff = *value - mean;
            diff * diff
        })
        .sum::<f64>()
        / samples as f64;

    Ok(MonteCarloEstimate {
        mean,
        variance,
        samples,
    })
}

#[cfg(test)]
mod tests {
    use super::{MonteCarloError, monte_carlo_mean};
    use use_seed::SimulationSeed;

    #[test]
    fn estimates_constant_observations() {
        let estimate = monte_carlo_mean(SimulationSeed::new(3), 5, |_sample| 2.0).unwrap();

        assert_eq!(estimate.mean, 2.0);
        assert_eq!(estimate.variance, 0.0);
        assert_eq!(estimate.samples, 5);
    }

    #[test]
    fn stays_repeatable_for_the_same_seed() {
        let first = monte_carlo_mean(SimulationSeed::new(11), 8, |sample| sample * 2.0).unwrap();
        let second = monte_carlo_mean(SimulationSeed::new(11), 8, |sample| sample * 2.0).unwrap();

        assert_eq!(first, second);
    }

    #[test]
    fn rejects_invalid_inputs() {
        assert_eq!(
            monte_carlo_mean(SimulationSeed::new(1), 0, |sample| sample),
            Err(MonteCarloError::InvalidSampleCount)
        );
        assert_eq!(
            monte_carlo_mean(SimulationSeed::new(1), 2, |_sample| f64::NAN),
            Err(MonteCarloError::NonFiniteObservation)
        );
    }
}