Transformations

A sequence of three transformations will be used to transform the original integral into an integral over a unit hyper-cube. This sequence begins with a Cholesky decomposition transformation $\mbox{\boldmath$\theta$}= C{\bf y}$, where $CC^t$ is the Cholesky decomposition of the covariance matrix $\Sigma$. Now $\mbox{\boldmath$\theta$}^t\Sigma^{-1}\mbox{\boldmath$\theta$}= {\bf y}^tC^tC^{-t}C^{-1}C{\bf y} = {\bf y}^t{\bf y}$, and $d\mbox{\boldmath$\theta$}= \vert C\vert d{\bf y} = \vert\Sigma\vert^{\frac{1}{2}}d{\bf y}$. Since ${\bf a} \leq \mbox{\boldmath$\theta$}= C{\bf y} \leq {\bf b}$ implies $(a_i - \sum_{j=1}^{i-1}c_{ij}y_j)/c_{ii} \leq y_i \leq (b_i - \sum_{j=1}^{i-1}c_{ij}y_j)/c_{ii}$ for $i = 1, 2, ..., m$, we have

$$F({\bf a, b}) = \frac{1}{\sqrt{(2\pi)^m}} \int_{a'_1}^{b'_1... ...1})}^{b'_m(y_1, ..., y_{m-1})} e^{-\frac{y_m^2}{2}} d{\bf y},$$

with $a'_i(y_1, ..., y_{i-1}) = (a_i - \sum_{j=1}^{i-1}c_{ij}y_j)/c_{ii}$ and $b'_i(y_1, ..., y_{i-1}) = (b_i - \sum_{j=1}^{i-1}c_{ij}y_j)/c_{ii}$.

Now each of the $y_i$'s can be transformed separately using $y_i = \Phi^{-1}(z_i)$, where

$$\Phi(y) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^y e^{-\frac{1}{2}\theta^2}d\theta$$

is the standard univariate normal distribution function. After these transformations, $F({\bf a, b})$ becomes

$$F({\bf a, b}) = \int_{d_1}^{e_1} \int_{d_2(z_1)}^{e_2(z_1)} ... ...t_{d_m(z_1, ..., z_{m-1})}^{e_m(z_1, ..., z_{m-1})} d{\bf z},$$

with $d_i(z_1, ..., z_{i-1}) = \Phi((a_i - \sum_{j=1}^{i-1}c_{ij}\Phi^{-1}(z_j))/c_{ii})$ and $e_i(z_1, ..., z_{i-1}) = \Phi((b_i - \sum_{j=1}^{i-1}c_{ij}\Phi^{-1}(z_j))/c_{ii})$.

The integrand in this form is much simpler than the original integrand. The integration region is more complicated, however, and cannot be handled directly with standard numerical multiple integration algorithms. A solution to this problem is to put the integral into a constant limit form using $z_i = d_i + w_i(e_i-d_i)$. After this final set of transformations,

$$F({\bf a, b}) = (e_1-d_1) \int_0^1 (e_2-d_2) \ ... \int_0^1 (e_{m}-d_{m}) \int_0^1 d{\bf w},$$

with $d_i = \Phi((a_i-\sum_{j=1}^{i-1}c_{ij} \Phi^{-1}(d_j+w_j(e_j-d_j)))/c_{ii})$ and $e_i = \Phi((b_i-\sum_{j=1}^{i-1}c_{ij} \Phi^{-1}(d_j+w_j(e_j-d_j)))/c_{ii})$.

The innermost integral over $w_m$ can be done explicitly because $d_m$ and $e_m$ have no dependence on $w_m$, so the complete sequence of transformations has reduced the number of integration variables by one.

The sequence of transformations described here might appear to have made the final integrand more complicated. However, the overall transformation has forced a priority ordering on the integration variables. The $w_1$ variable is the most important one, because all of the integrand factors $e_i-d_i$ for $i > 1$, depend on it. The $w_2$ variable is the next most important one, and so on. This priority ordering actually makes the problem more suitable for the type of subregion adaptive algorithm used for some of the tests described in Section 5. This type of adaptive algorithm subdivides along one coordinate axis at a time and therefore works most efficiently when the integrand has a few priority variables that define the directions of most variation in the integrand.