How to compute integral of multivariate normal with bounds defined by lines

Hi all,

I would like to use Distributions.jl and HCubature.jl (or suitable alternative) to compute an integral over part of the domain of a multivariate normal distribution. The image below provides an example. In the image, there are two lines and a contour of the density of a multivariate normal. How would I compute the integral of, for example, the top right quadrant (above the green line and to the right of the red line)? I don’t know if there is a clever transformation I should do, or if there is a clever technique I should use with HCubature.

Thank you

using Distributions 
using Plots 

line(x, β0, β1) = β0 + β1 * x 

μ = [0.0,0.0]
Σ = [1.0 0.0; 0.0 1.0]
dist = MvNormal(μ, Σ)

y1 = line.(x, .5, 3)
y2 = line.(x, 0, -.5)
v = range(-4, 4, length=200)
densities = [pdf(dist, [x,y]) for x ∈ v, y ∈ v]

contour(v, v, densities, leg=false, grid=false)
plot!(x, y1)
plot!(x, y2)


1 Like

I’ll try to put the description of the problem into formulas, to see if I understand the problem correctly, and provide a possible solution.

The problem

We have a 2-dimensional (probability-density) function p(x,y), defined on the whole plane \mathbb{R}^2 and we have two lines L_a and L_b, which we assume not to be parallel, such that they cross in one point and divide the plane in four quadrants. Let’s call the intersection point P = (x_p, y_p).

Each of the two line is either vertical(*) or we can represent it as a linear function:

  • a(x) = m_a (x - x_p) + y_p
  • b(x) = m_b (x - x_p) + y_p

In the original Cartesian coordinates of the (x,y)-plane, an integral over one quadrant can now be expressed as

I = \int_{x_p}^{\infty} \mathrm{d}x \int_{a(x)}^{b(x)} \mathrm{d}y ~ p(x, y)

where we assume that the quadrant we are loooking at is bounded in y-direction from below by L_a and from above by L_b. If we consider a different quadrant, we might have to switch the roles of x and y (e.g. integrate \mathrm{d}y from y_p to \infty) or switch a(x) and b(x).

(*) If one of the bounding lines is a vertical line, we can just replace the corresponding bound in the integral with \pm \infty and if it is horizontal, we can replace it with y_p (and analogously for the other quadrants).

“Quick and dirty”

In principle, this should define the integral in a way that can be used with numerical integration packages. For Cubature, we could do an integration in two consecutive steps (although this is quite inefficient and not recommended, see e.g. this discussion ).

Inefficient/hacky example

Note that I didn’t transform the x-integral here, which should be done for (semi-)infinite domains, since Cubature can only handle finite domains (see below). Just for demonstration purposes, I just chose a large value for the upper bound. If we choose it too large, however, the result will be wrong!

using Cubature

p(x, y) = exp(- (x^2 + y^2))
line(x, m, x_p, y_p) = m * (x - x_p) + y_p

function integrate_top_right_quadrant(x_p, y_p, m_a, m_b)

    a(x) = line(x, m_a, x_p, y_p)
    b(x) = line(x, m_b, x_p, y_p)

    yIntegral(x) = hquadrature(y -> p(x, y), a(x), b(x)) |> first
    I, err = hquadrature(yIntegral, x_p, 1.0e3)
    return I, err

m_a = -150.0
m_b = 150.0
x_p = 0.0
y_p = 0.0

# Very steep lines which should give a result close to the integral over the right half-plane, i.e. pi/2
integrate_top_right_quadrant(x_p, y_p, m_a, m_b)
# (1.5641297588910283, 1.5432516101686943e-10)

Transforming the integral

To deal with the integral properly (i.p. to use Cubature which requires finite, rectangular domains), we have to do some transformations, like you already mentioned.

There are some excellent illustrations on this page (I haven’t used the package myself, but it might be useful for you in general as well)

First, let’s shift the origin of the coordinate system, such that is is now at the intersection point P = (x_p, y_p). This only changes the argument of the function and the integral bounds (to keep the notation a bit simpler, we don’t bother with new variable names):

I = \int_{0}^{\infty} \mathrm{d}x \int_{a(x) - y_p}^{b(x) - y_p} \mathrm{d}y ~ p(x + x_p, y + y_p)

Next, we could look for a rotation and shear to transform the triangular-looking region of the plane such that the two lines L_a and L_b are horizontal and vertical, respectively. Then we wouldn’t have to worry about the x-dependent integral bounds anymore. However, this seems pretty cumbersome and since we’re going to use numerical integration anyway, we want to rather look for a simple transformation, regardless of whether the integrand stays simple or not.

After shifting the origin, the region to integrate is just a “slice of a pie”, so polar coordinates are very convenient to represent it, i.e. instead of denoting a point’s coordinates as (x-y)-pair, we give its distance from the origin and its angle (measured counter-clockwise) to the x-axis: (r, \varphi)
The relation between the coordinates is x = r \cos \varphi and y = r \sin \varphi and thus the integral then becomes (note the Jacobian determinant we have to include):

I = \int_{0}^{\infty} \mathrm{d}r \int_{\varphi_a}^{\varphi_b} \mathrm{d}\varphi ~ r p(r \cos \varphi + x_p, r \sin \varphi + y_p)

The angles \varphi_a = \arctan m_a and \varphi_b = \arctan m_b are the angles of the lines with respect to the x-axis and define the boundary of the angular integration.
Most importantly, they are now constants and don’t depend on the value of the second variable r anymore, so indeed have a rectangular domain of integration now and can use the Cubature package.

The only thing left for us to do is to get rid of the infinity in the r-domain (otherwise Cubature cannot handle the integral out of the box). There are many transformations that would achieve this, some might be better (performance-wise) than others, but let’s just try the one suggested in the docs of Cubature:

r(t) = \frac{t}{1-t}

The new integral with the adjusted bound and the Jacobian determinant is now

I = \int_{0}^{1} \mathrm{d}t \int_{\varphi_a}^{\varphi_b} \mathrm{d}\varphi ~ \frac{t}{(1-t)^3} p \left( \frac{t}{1-t} \cos \varphi + x_p, \frac{t}{1-t} \sin \varphi + y_p \right) \, .

This looks rather nasty when trying to solve it on paper, but numerically it should be much more convenient (minding the many possible opportunities for typos and the right expressions for the angles).

Code for transformed integral
using Cubature

p(x, y) = exp(- (x^2 + y^2))

function integrate_top_right_quadrant_polar(x_p, y_p, m_a, m_b)

    φ_a = atan(m_a)
    φ_b = atan(m_b)

    function integrand(point)
        t, φ = point
        r = t / (1 - t)
        return r / (1 - t)^2 * p(r * cos(φ) + x_p, r * sin(φ) + y_p)

    I, err = hcubature(integrand, [0.0, φ_a], [1.0, φ_b])
    return I, err

# `atan` can also handle `Inf`, so we can calculate the exact integral of half a Gaussian, for example
m_a = -Inf
m_b = Inf
x_p = 0.0
y_p = 0.0

integrate_top_right_quadrant_polar(x_p, y_p, m_a, m_b)
# (1.5707963268678242, 1.3329275146601752e-8)
#  1.5707963267948966 == π/2

Of course, whether the same solution works for you depends on how much of the actual problem stays the same (is it always just to lines in 2D or do you want to generalize it, etc), but I hope this illustrates some possible solutions.

Also note that I haven’t tested the code thoroughly and there are definitely some parameter regimes where results of the integral are wrong (because the slopes of the lines are too big/small for example or because some assumptions from above are violated).


My take would be :

  1. use a change of variables to reduce your integral to \mathbb P(X>x, Y> y), where (X,Y) is a bivariate gaussian random vector with a given mean and given variance (which both come from the change of variables), and given x,y.

  2. MvNormalCdf.jl has very good routines.


Good point! As long as the integration region is given by two straight lines, the transformed distribution should still be Gaussian.

1 Like

Thank you both for your responses. This is very helpful.

For my use case, the dimensionality will be 2 and the two boundaries will always be linear. The only thing that can change is which region to integrate. It looks like I should be able to figure that out from what you have.

One question: the function integrate_top_right_quadrant_polar returns 1.57, but I expect the solution to be <= 1 because MvNormal is a probability distribution. Is that because I need to apply an inverse transformation or were you demonstrating the approach for a simpler problem?

1 Like

Ah yeah, sorry. I was lazy and just used the exponential without the normalization and with standard deviation equal to 1. The Gaußian integral in 2D is \pi, so if the two lines are chosen such that the region covers half the plane, the result should be roughly \pi / 2.

If you use an actual PDF for the function p(x,y) you should not need to add any additional factors. They already come out of the variable transformation (as the Jacobian determinant).

1 Like

Fleshing out what Oskar suggests, note that the region of integration between the two lines is the image of the upper right quadrant \mathbb{R}^2_{\geq 0} under a liner map, i.e. I \ni x = Su with S the matrix whose columns are the vectors that span the wedge (here you need to be a bit careful with orientation).

Given a multivariate normal \mathcal{N}(\mu, \Sigma) with pdf p, the integral of interest is

\int_I p(x)\, dx = \int_{\mathbb{R}^2_{\geq 0}} (p\circ S)(u) |\det{S}| \, du
  1. The new integrand is the pdf of a normal distribution \mathcal{N}(S^{-1}\mu, S^{-1}\Sigma (S^{-1})^t) which is easily constructed. (The new covariance matrix is the inverse of S^t\Sigma^{-1}S )
  2. The new domain of integration is rectangular, albeit unbounded.

The second point poses a problem in theory since one cannot strictly integrate numerically to infinity. One certainly clean solution is to use a range transform as was implemented in Sevi’s answer. However, it is perfectly sufficient to integrate up to some finite distance from the origin; maybe 6 or more standard variations to be safe, or even a lot more. It won’t incur much of a cost, because an adaptive integration algorithm is going to take huge steps far out, because the normal distribution decays very fast.

julia> S = [ [cos(pi/6), sin(pi/6)];; [cos(pi/2-0.2), sin(pi/2-0.2)] ]
2×2 Matrix{Float64}:
 0.866025  0.198669
 0.5       0.980067

julia> Sinv = inv(S);

julia> D2 = MvNormal(Sinv*Sinv')
dim: 2
μ: Zeros(2)
Σ: [1.7804931170451757 -1.1788395237690619; -1.1788395237690619 1.7804931170451754]

julia> @time int_upper_right(D2)
  0.002839 seconds (27.10 k allocations: 1.291 MiB)
(0.1348356780494405, 1.996183543386641e-9)

Much appreciated!

I don’t have a strong calculus background. So it will take me a bit to fully understand, but I think I understand the gist. The explanations and code you guys have provided help.

1 Like

Yes, these approaches are difficult to grasp fully without a bit of multivariate calculus and matrix algebra.

Note my edit! I’ve originally made a somewhat critical blunder and gave the wrong covariance matrix for the transformed distribution… :grimacing:

Another way is to use an indicator function, i.e. a function that is 1 if its argument is within the domain of interest, and 0 else. HCubature on its own seems to struggle, at least it was taking forever in my test.
But it works decently with DomainIntegrals.jl and DomainSets.jl

# slopes of the bounding lines; same as before
julia> m1 = tan(pi/6); m2 = tan(pi/2-0.2);

julia> domain = Domain(y>=m1*x && y<=m2*x for (x,y) in Rectangle(0..100, 0..100))
indicator function bounded by: (0 .. 100) × (0 .. 100)

julia> D = MvNormal(diagm([1.0,1.0]))
dim: 2
μ: Zeros(2)
Σ: [1.0 0.0; 0.0 1.0]

julia> h(x) = pdf(D, x)
h (generic function with 1 method)

julia> integrate(h, domain)
(0.13483828322767, 0.00112992476806625)

It gives a larger error than the previous approach, but is certainly more tangible.

1 Like

I used this approach naively as my first attempt. I also found that HCubature was struggling to compute the integral. In fact, I let it run for about 20 seconds before manually interrupting. My guess it failed to adapt the correct step size.

For larger numbers of lines and/or higher dimensions, you may find Monte Carlo simulation to be more efficient for this sort of task. The integral over all of \mathbb{R}^n is always 1, so the value of the integral in question is just the probability that a draw from the normal distribution falls in between the lines of interest. Simulating a large number of draws and computing the sample probability will give you a good estimate of the value of the integral.


I also wanted to suggest this initially, but at least you cannot make the upper bound of the region arbitrarily large with the default settings of Cubature, as (I suspect) the algorithm at some point thinks the function is zero everywhere in the region and just silently gives the wrong result:

julia> using Cubature

julia> hquadrature(x -> exp(-x^2), 0.0, 10.0)
(0.8862269254527582, 1.2023665832673908e-9)

julia> hquadrature(x -> exp(-x^2), 0.0, 100.0)
(0.8862269254527579, 2.3374315648411976e-9)

julia> hquadrature(x -> exp(-x^2), 0.0, 1000.0)
(0.886226925452758, 9.287832678601159e-10)

julia> hquadrature(x -> exp(-x^2), 0.0, 10000.0)
(0.0, 0.0)

The region I had to chose here to run into the issue is definitely a lot larger than six standard deviations, but I still feel a bit uncomfortable, knowing that the result might be nonsensical if I chose a bad value for the bound. In any case, O(10) standard deviations should be safe! :sweat_smile:

Not sure if the same thing happens with the other integration packages, but I think anything that uses quadrature/cubature internally will run into this issue?

I have some fuzzy recollection of this behavior occurring with Quadgk.jl. Selecting the right bounds was precarious because I wasn’t sure what would trigger a silent error.

1 Like

Yes, you are right. If one makes the integration rectangle too large, the algorithm will fail. This can be mitigated to some degree by choosing a larger initdiv in hcubature, but ultimately one has to chose upper bounds based on the behavior of the integrand.
A surefire way to tell that the integration result is nonsensical, is when the result and reported error are exactly zero. But I fully agree that a range transformation is the cleaner solution.