It seems like it is lacking the functionality R's integrate has for handling infinite boundaries, but I suppose you could implement that yourself on the outside.

For what it's worth,

    use integrate::adaptive_quadrature::simpson::adaptive_simpson_method;
    use statrs::distribution::{Continuous, Normal};

    fn dnorm(x: f64) -> f64 {
        Normal::new(0.0, 1.0).unwrap().pdf(x)
    }
    
    fn main() {
        let result = adaptive_simpson_method(dnorm, -100.0, 100.0, 1e-2, 1e-8);
        println!("Result: {:?}", result);
    }

prints Result: Ok(1.000000000053865)

It does seem to be a usability hazard that the function being integrated is defined as a fn, rather than a Fn, as you can't pass closures that capture variables, requiring the weird dnorm definition

You will be completely blown away, then, from what Wolfram Language (aka Mathematica) can do. (When it comes to numerical integration.)

https://reference.wolfram.com/language/tutorial/NIntegrateOv...

for ]-inf, inf[ integrals, you can use Gauss Hermite method, just keep in mind to multiply your function with exp(x^2).

    use integrate::{
        gauss_quadrature::hermite::gauss_hermite_rule,
    };
    use statrs::distribution::{Continuous, Normal};

    fn dnorm(x: f64) -> f64 {
        Normal::new(0.0, 1.0).unwrap().pdf(x)* x.powi(2).exp()
    }

    fn main() {
        let n: usize = 170;
        let result = gauss_hermite_rule(dnorm, n);
        println!("Result: {:?}", result);
    }

I got Result: 1.0000000183827922.

How many evaluations of the underlying function does it make? (Hoping someone will fire up their R interpreter and find out.)

Or, probably, dnorm is a probability distribution which includes a likeliness function, and a cumulative likeliness function, etc. I bet it doesn't work on arbitrary functions.