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

General Usage
, , ,
2

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
end

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)
https://jverzani.github.io/CalculusWithJuliaNotes.jl/integral_vector_calculus/double_triple_integrals.html

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)
    end

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

# `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).

4 Likes