# Calculating multivariate integrals between lower and upper bounds

Suppose $$\vec{X}=(x_1,x_2,...,x_n)$$ follows some continuous multivariate distribution, such that $$x_i\in{\rm I\!R}, i=1,...,n$$.

Suppose also that I have access to the following functions:

• $$\phi(\vec{x})$$, which gives me the pdf at point $$\vec{x}=(x_1,x_2,...,x_n)$$
• $$\Phi^{upper}(\vec{x})$$, which calculates the CDF from $$-\infty$$ to $$\vec{x}$$. More specifically: $$\Phi^{upper}(\vec{x})=\int_{\vec{k}=(-\infty,...,-\infty)}^{(x_1,...,x_n)}\phi(\vec{k})d\vec{k}$$

My main question is: how can I use $$\Phi^{upper}(\vec{x})$$ to calculate the CDF between lower and upper bounds?

In other words, how can I calculate integral of $$\phi$$ from $$\vec{x}_{lower}$$ to $$\vec{x}_{upper}$$ using $$\Phi^{upper}$$?

In mathematical terms, how can I use $$\Phi^{upper}$$ to calculate the following integral?

$$\Phi^{upper}_{lower}(\vec{x}_{lower},\vec{x}_{upper})=\int_{\vec{k}=(x_1^{lower},...,x_n^{lower})}^{(x_1^{upper},...,x_n^{upper})}\phi(\vec{k})d\vec{k}$$

I know that if there are only 2 dimensions in $$\vec{X}$$, I can calculate $$\Phi^{upper}_{lower}(\vec{x}_{lower},\vec{x}_{upper})$$ like this:

$$\Phi^{upper}_{lower}(\vec{x}_{lower},\vec{x}_{upper}) = \Phi^{upper}(x_1^{upper},x_2^{upper}) - \Phi^{upper}(x_1^{upper},x_2^{lower}) - \Phi^{upper}(x_1^{lower},x_2^{upper}) + \Phi^{upper}(x_1^{lower},x_2^{lower})$$

However, I don't know how to generalize this result for an N-dimensional case.

Does anyone have any insight on how to do this?

Thank you very much!

• Not exactly, no. Consider the 2 dimensional example I gave. If you just do $\Phi^{upper}(\vec{x}_{upper}) - \Phi^{upper}(\vec{x}_{lower})$, you're leaving in a few areas that shouldn't be considered. Here is an image to illustrate my point. – Felipe D. Jul 22 '18 at 3:35
• Assuming you're talking about a rectangular area, you probably want to apply the inclusion-exclusion rule here. – Ben Jul 22 '18 at 3:38
• I can see how this touches on what I'm getting at, but I'm having a bit of a hard time understanding how I would apply this concept from set theory to this particular problem. Could you maybe elaborate a little bit further on how to do so? Thanks! – Felipe D. Jul 22 '18 at 3:45

This problem involves the application of a probability measure over a union of non-disjoint sets, and so it can be solved by application of the inclusion-exclusion rule.

To facilitate this analysis, we will let $\mathcal{A}(x) \equiv (-\infty, x]$ denote the real numbers up to $x$ so that $\mathcal{A}(x_1) -\mathcal{A}(x_0) = (x_0,x_1]$ is a bounded interval. (If you are dealing with a continuous random variable you do not have to worry about whether the ends of the intervals are open or closed.) We will also simplify the problem by using a slight abuse of notation, treating each $\mathcal{A}(x_i)$ as a subset of $\mathbb{R}^n$ with an upper bound in dimension $i$ and free in the other dimensions. This means that for an input vector $\mathbf{x} = (x_1,...,x_n)$ you have CDF values of the form:

$$\Phi(\mathbf{x}) = \mathbb{P} (\mathbf{X} \leqslant \mathbf{x}) = \mathbb{P} \Bigg( \bigcap_{k=1}^n \mathcal{A}(x_k) \Bigg).$$

You want to calculate the probability of a rectangular area, which can be written as an intersection of bounded intervals as $\mathcal{R}_n = \bigcap_{k=1}^n (\mathcal{A}(\bar{x}_{k}) -\mathcal{A}(\underline{x}_{k}))$, where you have lower and upper bounds $\underline{\mathbf{x}} < \bar{\mathbf{x}}$. You want to be able to write the probability of this rectangular event using your CDF $\Phi$. To apply the inclusion-exclusion rule we will let $\mathfrak{N}_k$ denote the class of all subsets of $\{ 1,...,n \}$ with exactly $k$ elements. Using this rule, and some other set algebra, we have:

\begin{equation} \begin{aligned} \mathbb{P}(\mathcal{R}_n) &= \mathbb{P} \Bigg( \bigcap_{k=1}^n (\mathcal{A}(\bar{x}_{k}) -\mathcal{A}(\underline{x}_{k})) \Bigg) \\[6pt] &= \mathbb{P} \Bigg( \bigcap_{k=1}^n \mathcal{A}(\bar{x}_{k}) \Bigg) - \mathbb{P} \Bigg( \bigcap_{k=1}^n \mathcal{A}(\bar{x}_{k}) \cap \bigcup_{k=1}^n \mathcal{A}(\underline{x}_{k}) \Bigg) \\[6pt] &= \mathbb{P} \Bigg( \bigcap_{k=1}^n \mathcal{A}(\bar{x}_{k}) \Bigg) - \sum_{k=1}^n \Bigg[ (-1)^{k-1} \sum_{\mathcal{D} \in \mathfrak{N}_k} \mathbb{P} \Bigg( \bigcap_{i \notin \mathcal{D}} \mathcal{A}(\bar{x}_{i}) \cap \bigcap_{i \in \mathcal{D}} \mathcal{A}(\underline{x}_i) \Bigg) \Bigg] \\[6pt] &= \Phi (\bar{\mathbf{x}}) - \sum_{k=1}^n (-1)^{k-1} \sum_{\mathcal{D} \in \mathfrak{N}_k} \Phi(\mathbf{x_\mathcal{D}}), \\[6pt] \end{aligned} \end{equation}

where the data vector $\mathbf{x}_\mathcal{D}$ uses the lower bounds $\underline{x}_i$ for all $i \in \mathcal{D}$ and the upper bounds $\bar{x}_i$ for all $i \notin \mathcal{D}$. This gives you a general mathematical form for calculating the probability of a rectangle directly using a multivariate CDF.

Special cases: Application of the general rule yields and expression with $2^n$ terms. For small $n$ we can write this out explicitly without using summation notation and it is not too long. For $n=2$ we get the special case:

$$\mathbb{P}(\mathcal{R_2}) = \Phi(\bar{x}_1, \bar{x}_2) - \Phi(\underline{x}_1, \bar{x}_2) - \Phi(\bar{x}_1, \underline{x}_2) + \Phi(\underline{x}_1, \underline{x}_2).$$

For $n=3$ we get the special case:

\begin{equation} \begin{aligned} \mathbb{P}(\mathcal{R_3}) &= \Phi (\bar{x}_1, \bar{x}_2, \bar{x}_3) - \Phi(\underline{x}_1, \bar{x}_2, \bar{x}_3) - \Phi(\bar{x}_1, \underline{x}_2, \bar{x}_3) - \Phi(\bar{x}_1, \bar{x}_2, \underline{x}_3) \\[4pt] &\quad + \Phi(\underline{x}_1, \underline{x}_2, \bar{x}_3) + \Phi(\underline{x}_1, \bar{x}_2, \underline{x}_3) + \Phi(\bar{x}_1, \underline{x}_2, \underline{x}_3) - \Phi(\underline{x}_1, \underline{x}_2, \underline{x}_3). \end{aligned} \end{equation}

For larger $n$ the number of terms increases exponentially and so it becomes cumbersome to write the expression out without using the summation notation in the general form.

• Thanks for the clear example. The last special case really helped me see what was going on in those nested summations. Just got a python code up and running and spitting out the correct results. Thanks again! – Felipe D. Jul 22 '18 at 20:38
• That's great that you got the code working. The easiest way to code something like this is to sum over all the $2^n$ binary values of length $n$, using the binary values as representations of whether you are using an upper or lower bound for an argument in the CDF. – Ben Jul 22 '18 at 22:17

As suggested by @Ben, this problem boils down to an application of the inclusion-exclusion principle. For notational simplicity, let's replace $x_i^{lower}$ with $0$, and $x_i^{upper}$ with $1$, so that for example the vector $\vec{x}=(x_1^{lower},x_2^{upper},x_3^{upper},x_4^{lower})$ becomes $\vec{x}=(0,1,1,0)$.

Then the general formula for $\Phi^{upper}_{lower}(\vec{x}_{lower},\vec{x}_{upper})$ with $n$-dimensional $\vec{x}$ is:

$$\Phi^{upper}_{lower}\big((0,\dots,0),\ (1,\dots,1)\big) \ = \ \Phi^{upper}(1,\dots,1) \ - \sum_{x_1+\cdots+x_n=n-1} \Phi^{upper}(x_1,\dots,x_n) \ + \cdots \\+ \ (-1)^{n-1} \cdot \sum_{x_1+\cdots+x_n=1} \Phi^{upper}(x_1,\dots,x_n) \ + \ (-1)^n \cdot\Phi^{upper}(0,\dots,0)$$

• Thanks for clearing things up! But I'm still a bit confused. First of all, if we only have two dimensions, $\vec{x}=(x_1,x_2)$. So I believe you meant to say that $\vec{x}_{lower}=(0,0)$ and $\vec{x}_{upper}=(1,1)$. Or more generally, for an n-dimensional case, $\vec{x}_{lower}=(0,...,0)$ and $\vec{x}_{upper}=(1,...,1)$. But besides that, I'm still getting confused with the limits of the summations. What do you mean when you use $x_1+...+x_n=n-1$ as the limits/bounds of your summation? Is it summing all possible combinations where $x_1+...+x_n=n-1$ is true? – Felipe D. Jul 22 '18 at 9:23
• Furthermore, what are the terms $x_1, ..., x_n$ in your example? Are they ones or zeros? This seems a bit ambiguous. Thanks again for the help! – Felipe D. Jul 22 '18 at 9:27
• @FelipeD. My reference to $\vec{x}=(x_1^{lower},x_2^{upper},x_3^{upper},x_4^{lower})$, which I assume is what your confusion is about in your first comment, was just meant to show an example of the notation that I'm using, where it just so happens that $n=4$. You're correct that using this notation $\vec{x}_{lower}=(0,\dots,0)$ and $\vec{x}_{upper}=(1,\dots,1)$. You are also correct in your understanding of the summation notation, that it's referring to all combinations of $x_i$ such that the equality holds. This is a commonly-used notational convention. – jon_simon Jul 22 '18 at 18:33
• Regarding your second comment, yes, all of the $x_i$ are either $0$ (denoting $x_i^{lower}$) or $1$ (denoting $x_i^{upper}$) – jon_simon Jul 22 '18 at 18:36