support integrand function

[CarloLucibello]

Hi, thanks for this package which I only recently discovered. It is often the case with some types of calculations I do that I have the necessity to integrate some function g(x) with a Gaussian measure:

\int_{a_1}^{b_1} dx_1\dots\int_{a_n}^{b_n} dx_n\  \frac{1}{\sqrt{(2\pi)^n \det\Sigma}} e^{-\frac{x^T \Sigma^1 x}{2}}\ g(x)

Genz's paper discusses a minimal modification of his algorithm to deal with an integrated function at the end of Section 3.

Taking inspiration from this package I implemented the method, trying to keep the code as simple as possible and relying on SpecialFunctions.jl for erf and erfinv.

The code can be found at https://gist.github.com/CarloLucibello/c3f3196f3ed89bbc0f296151f32dba0e

I thought it could be useful to share it here.

It can be used as follows:

julia> Σ = [4 3 2 1 
            3 5 -1 1 
            2 -1 4 2 
            1 1 2 5];

julia> a = [-Inf, -Inf ,-Inf, -Inf];

julia> b = [-1, -2, 4, 1];

julia> g1(x) = 1
f1 (generic function with 1 method)

julia> val, err = ∫D(g1, Σ, a, b; nevals=100_000)
(0.10550910246523065, 1.3200680774522703e-5)

julia> g(x) = log(1 + sum(abs2, x))

julia> val, err = ∫D(g, Σ, a, b; nevals=100_000)
(0.3434052358363517, 4.420996703598769e-5)

Results are consistent with MvNormcalCDF.mvnormcdf when g1(x) = 1 is used, although my function appears to be 3x slower.

I hope this will be helpful to someone.

Best, Carlo

[CarloLucibello]

Maybe since it is a simple extension to the code already here, this could be integrated in the package.

[PharmCat]

Maybe since it is a simple extension to the code already here, this could be integrated in the package.

Hi! Thank you! I think it can bu useful. Pull request with tests welcome)