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!
 A: 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.
A: 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)$$
