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

**URL:** <https://discourse.julialang.org/t/how-to-compute-integral-of-multivariate-normal-with-bounds-defined-by-lines/104160>\
**Category:** General Usage\
**Tags:** question, calculus, distributions, integral\
**Created:** [September 22, 2023, 5:40pm UTC](https://discourse.julialang.org/t/how-to-compute-integral-of-multivariate-normal-with-bounds-defined-by-lines/104160 "2023-09-22T17:40:00Z")\
**Posts on this page:** 1\
**Showing post:** 2

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [September 23, 2023, 8:14pm UTC](https://discourse.julialang.org/t/how-to-compute-integral-of-multivariate-normal-with-bounds-defined-by-lines/104160/2 "2023-09-23T20:14:28Z")

</div>

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](https://discourse.julialang.org/t/2d-integration-over-non-rectangular-domain-using-cubature/2991) ).

> **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!
> 
> ```julia
> 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)

> **[63  Multi-dimensional integrals – Calculus with Julia](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](https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals):

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**
>
> ```julia
> 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).

---

_[View the full topic](https://discourse.julialang.org/t/how-to-compute-integral-of-multivariate-normal-with-bounds-defined-by-lines/104160)._