main.rs (1468B)
1 use rand_core::SeedableRng; 2 use rand_distr::{Distribution, LogNormal, Uniform}; 3 use rand_pcg::Pcg64Mcg; 4 5 fn sample_to(low: f64, high: f64, mut rng: &mut Pcg64Mcg) -> f64 { 6 let normal90 = 1.6448536269514727; // change to global const later 7 let loglow = low.ln(); 8 let loghigh = high.ln(); 9 let normal_mean = (loghigh + loglow)/2.0; 10 let normal_std = (loghigh - loglow) / (2.0 * normal90); 11 12 let lognormal = LogNormal::new(normal_mean, normal_std).unwrap(); 13 let x = lognormal.sample(&mut rng); 14 // https://docs.rs/rand_distr/latest/src/rand_distr/normal.rs.html#232-236 15 x; 16 } 17 18 fn model(mut rng: &mut Pcg64Mcg) -> f64 { 19 let p_a = 0.8; 20 let p_b = 0.5; 21 let p_c = p_a * p_b; 22 23 let ws = [1.0 - p_c, p_c / 2.0, p_c / 4.0, p_c / 4.0]; 24 25 let uniform = Uniform::new(0.0, 1.0); /* I don't understand why this doesn't need to be unwrapped, unlike normal below */ 26 let p: f64 = uniform.sample(&mut rng); 27 if p < ws[0] { 28 0.0; 29 } else if p < (ws[0] + ws[1]) { 30 1.0; 31 } else if p < (ws[0] + ws[1] + ws[2]) { 32 sample_to(1.0, 3.0, rng); 33 } else { 34 sample_to(2.0, 10.0, rng); 35 } 36 } 37 38 fn main() { 39 let mut rng = Pcg64Mcg::seed_from_u64(1); 40 41 let mut mean = 0.0; 42 let n_samples = 1_000_000; 43 for _ in 0..n_samples { 44 let x = model(&mut rng); 45 // println!("Sample: {}", x); 46 mean += x; 47 } 48 mean = mean/(n_samples as f64); 49 println!{"Mean: {}", mean}; 50 }