sine/
sine.rs

1//! Explicit univariate biparametric sinusoidal model.
2
3use oefpil::{
4    Algorithm, Covariance, Distribution, Logfile, Model, OefpilError, Parameter, Variable,
5};
6use rand_chacha::{ChaChaRng, rand_core::SeedableRng};
7
8#[allow(clippy::float_cmp)]
9fn main() -> Result<(), OefpilError> {
10    // Define explicit model of one independent variable and two parameters.
11    const MODEL: Model = Model {
12        fx: |x, p| p[1] * (x[0] + p[0]).to_radians().sin(),
13        dfdx: &[|x, p| p[1] * (x[0] + p[0]).to_radians().cos().to_radians()],
14        dfdp: &[
15            |x, p| p[1] * (x[0] + p[0]).to_radians().cos().to_radians(),
16            |x, p| (x[0] + p[0]).to_radians().sin(),
17        ],
18        implicit: false,
19    };
20    // Define true values.
21    let p_mean = [8f64, 10.0];
22    let x_mean = [-60.0, -40.0, -20.0, -10.0, 10.0, 20.0, 40.0, 60.0];
23    let y_mean = x_mean.map(|x_mean| (MODEL.fx)(&[x_mean], &p_mean));
24    // Define covariance matrix with constant Pearson `correlation` `coefficient` (PCC).
25    let x_deviation = x_mean.map(|x| x.abs() * 2e-2);
26    let y_deviation = y_mean.map(|y| y.abs() * 8e-2);
27    let coefficient = 0.5;
28    let correlation = x_deviation
29        .iter()
30        .zip(&y_deviation)
31        .map(|(x, y)| x * y * coefficient)
32        .collect::<Vec<f64>>();
33    let covariance = &Covariance::new_diagonals(8, 2)
34        .with_tile(0, 0, &x_deviation.map(|x| x * x))?
35        .with_tile(1, 1, &y_deviation.map(|y| y * y))?
36        .with_tile(0, 1, &correlation)?
37        .with_decomposition()?;
38    // Seed of 256 bits for deterministic high-quality randomness.
39    let seed = b"85a4d500a53d81b65255dc718d7df63a3cd3d6bb15e6f724af04a23d74cc02d3";
40    // Cryptographically secure pseudorandom number generator (CSPRNG).
41    let mut csprng = ChaChaRng::from_seed(const_hex::decode_to_array(seed).unwrap());
42    // Define correlated normal distributions of variables.
43    let mut distribution = Distribution {
44        sample: &mut vec![0f64; x_mean.len() + y_mean.len()],
45        mean: &mut x_mean.iter().chain(&y_mean).copied().collect::<Vec<f64>>(),
46        covariance,
47    };
48    // Simulate measurement by sampling observations from distributions of variables.
49    distribution.sample(&mut csprng)?;
50    // Pretend mediocre estimates of variable means (empty slices default to samples).
51    let mut variable = Variable {
52        sample: distribution.sample,
53        mean: &mut [],
54        covariance,
55    };
56    // Pretend mediocre estimates of parameter means and assign output slices.
57    let p_mean_initial = [9.5f64, 12.0];
58    let mut parameter = Parameter {
59        mean: &mut p_mean_initial.clone(),
60        deviation: &mut [0f64; 2],
61        covariance: &mut [0f64; 4],
62        correlation: &mut [0f64; 4],
63    };
64    // Perform fitting with default settings.
65    let oefpil = Algorithm::default();
66    let report = oefpil.fit(MODEL, &mut variable, &mut parameter, Logfile::default())?;
67    // Assert expected results.
68    let round = |value: f64| (value * 1e2).round() / 1e2;
69    assert_eq!(7, report.iterations);
70    assert_eq!(0.73, report.chi_squared_p_value().map(round)?);
71    // `parameter.mean` deviates from `p_mean` in agreement with `parameter.deviation`.
72    assert_eq!(8.16, round(parameter.mean[0]));
73    assert_eq!(0.19, round(parameter.deviation[0]));
74    assert_eq!(9.95, round(parameter.mean[1]));
75    assert_eq!(0.28, round(parameter.deviation[1]));
76    assert_eq!(0.30, round(parameter.correlation[1]));
77    // Confirm `parameter.deviation` via Monte Carlo simulation.
78    let samples = 10_000;
79    let mut p_sample = vec![0f64; samples * p_mean.len()];
80    let mut p_correlation = vec![0f64; samples];
81    for sample in 0..samples {
82        distribution.sample(&mut csprng)?;
83        let mut variable = Variable {
84            sample: distribution.sample,
85            mean: &mut [],
86            covariance,
87        };
88        let mut parameter = Parameter {
89            mean: &mut p_mean_initial.clone(),
90            covariance: &mut [0f64; 4],
91            deviation: &mut [0f64; 2],
92            correlation: &mut [0f64; 4],
93        };
94        oefpil.fit(MODEL, &mut variable, &mut parameter, Logfile::default())?;
95        p_sample[sample] = parameter.mean[0];
96        p_sample[sample + samples] = parameter.mean[1];
97        p_correlation[sample] = parameter.correlation[1];
98    }
99    let mean_deviation = |s: &[f64]| {
100        let n = f64::from(u32::try_from(s.len()).unwrap());
101        let m = s.iter().sum::<f64>() / n;
102        let v = s.iter().map(|s| s - m).map(|d| d * d).sum::<f64>() / (n - 1.0);
103        [m, v.sqrt()]
104    };
105    let [p_mean_0, p_deviation_0] = mean_deviation(&p_sample[..samples]);
106    let [p_mean_1, p_deviation_1] = mean_deviation(&p_sample[samples..]);
107    let [p_correlation_mean, p_correlation_deviation] = mean_deviation(&p_correlation);
108    // Simulation agrees well with `p_mean` and `parameter.deviation`.
109    assert_eq!(8.00, round(p_mean_0));
110    assert_eq!(0.19, round(p_deviation_0));
111    assert_eq!(10.0, round(p_mean_1));
112    assert_eq!(0.28, round(p_deviation_1));
113    assert_eq!(0.30, round(p_correlation_mean));
114    assert_eq!(0.02, round(p_correlation_deviation));
115    Ok(())
116}
sine/sine.rs

sine/
sine.rs
