Skip to main content

bayesian_quadrature/
bayesian_quadrature.rs

1// Infer an integral, and its posterior uncertainty, from a few function evaluations.
2//
3// Run with: cargo run --example bayesian_quadrature
4
5use uncertain_numerics::{BayesianQuadrature, GaussianMeasure, RbfKernel};
6
7fn main() -> Result<(), Box<dyn std::error::Error>> {
8    // Prior over the integrand: zero-mean Gaussian process with an RBF kernel.
9    let kernel = RbfKernel::new(1.0, 1.0)?; // signal variance, length scale
10    // Integration measure p(x) = N(0, 1).
11    let measure = GaussianMeasure::new(0.0, 1.0)?;
12    // Jitter is an explicit, fixed diagonal regularizer. It is never escalated silently.
13    let quadrature = BayesianQuadrature::new(kernel, measure, 1.0e-10);
14
15    // Observe f(x) = cos(x) at seven nodes. Exactly, E[cos X] = exp(-1/2) for X ~ N(0, 1).
16    let nodes = [-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0];
17    let values: Vec<f64> = nodes.iter().copied().map(f64::cos).collect();
18    let exact = (-0.5_f64).exp();
19
20    let posterior = quadrature.posterior(&nodes, &values)?;
21    let standardized_error = (posterior.mean() - exact).abs() / posterior.standard_deviation();
22
23    println!("E[I | y]           = {:.6}", posterior.mean());
24    println!(
25        "sd[I | y]          = {:.3e}",
26        posterior.standard_deviation()
27    );
28    println!("exact integral     = {exact:.6}");
29    println!("standardized error = {standardized_error:.3}");
30
31    // The posterior is honest about its own error here: the exact value lies well
32    // inside the reported uncertainty.
33    assert!((posterior.mean() - exact).abs() < 3.0 * posterior.standard_deviation());
34    Ok(())
35}