Separation-of-Variables Methods

Genz and Bretz (1999) continued the separation of variables (SV) transformation that results in equation (7) with additional transformations (similar to those used by Genz, 1992, for MVN problems), to produce methods based on a complete transformation of the original integration to the unit hypercube $[0,1]^m$. The first step is to notice that the lower triangular structure of $C$ allows a separation of the integration limits in equation (7), so that

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath$\Sigma$}, \nu ) = ... ...}^{(1)}}{(1+\frac{u_m^2}{m+\nu-1})^{\frac{m+\nu}{2}}}d{\bf u};$$ (11)

with ${\tilde a}_i=\sqrt{\frac{\nu+i-1}{\nu+\sum_{j=1}^{i-1}y_j^2}} (a_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i},$ and ${\tilde b}_i=\sqrt{\frac{\nu+i-1}{\nu+\sum_{j=1}^{i-1}y_j^2}} (b_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i}$ .

Next, let $u_i=t_{\nu+i-1}^{-1}(z_i)$, and then

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath $\Sigma$},\nu ) = ... ...}_m(t_\nu^{-1}(z_1),...,t_{\nu+m-2}^{-1}(z_{m-1})))} d{\bf z}.$$

Finally, let $z_i=d_i+(e_i-d_i)w_i,$ so $dz_i=(e_i-d_i)dw_i,$ with

$$\begin{eqnarray*} d_i(w_1,..., w_{i-1})&=&t_{\nu+i-1} ({\tilde a}_i(t_\nu^{-1}(z... ...i(t_\nu^{-1}(z_1(w_1)),...,t_{\nu+i-2}^{-1}(z_{i-1}(w_{i-1})))). \end{eqnarray*}$$

Then

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath $\Sigma$},\nu ) = ... ...,w_{m-1})-d_{m}(w_1,\ldots,w_{m-1})) \int\limits_0^1 d{\bf w}.$$

The innermost integral has value equal to one, so the number of integration variables has been reduced to $m-1$, and standard multidimensional numerical integration methods can be used for the transformed ${\bf T}$. A simple Monte-Carlo algorithm for the MVT problem uses

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath$\Sigma$},\nu ) \ap... ...}{N} \sum_{k=1}^N \prod_{i=1}^m(e_i({\bf w}_k)-d_i({\bf w}_k))$$ (12)

with $w_{i,k} \sim Uniform(0,1)$.

Schervish (1984) originally suggested that the computation of MVN probabilities should be easier for numerical integration methods if the variables are reordered (and appropriate rows and columns of $\mbox{\boldmath$\Sigma$}$ are permuted) so that the innermost integrals have the larger integration intervals. This sorting heuristic often has the effect that the innermost integrals have expected value closer to one, thereby reducing the overall variation in the integrand. Gibson, Glasbey and Elston (1992), who independently developed a Monte-Carlo method similar to the one developed by Genz (1992) for MVN probabilities, suggested an improved prioritization of the variables. With this technique, the variables are sorted so that the innermost integrals have the largest expected integration intervals. This is more complicated than just sorting the integration limits and permuting the respective rows and columns of $\mbox{\boldmath$\Sigma$}$, because the Cholesky factor $C$ must be computed dynamically during the sorting of the variables. This method uses ${\bf a}$, ${\bf b}$ and $\mbox{\boldmath$\Sigma$}$ in the sorting process, and it should therefore further increases the likelihood that the innermost integrals have values close to one and improve the convergence of the numerical integration methods.

The Gibson, Glasbey and Elston method can be generalized to MVT problems, which use the SV transformations, in the following manner, using the MVT definition given by equation (11). The first (outermost) integration variable is chosen by selecting a variable $i$ where $\min_{1\leq i \leq m}\{ t_\nu(b_i/\sqrt{\sigma_{i,i}})-t_\nu(a_i/\sqrt{\sigma_{i,i}})\}$ is achieved. The limits and rows and columns of $\mbox{\boldmath$\Sigma$}$ for variables $1$ and $i$ are interchanged. Then the first column of the Cholesky decomposition $C$ of $\mbox{\boldmath$\Sigma$}$ is computed using $c_{1,1} = \sqrt{\sigma_{1,1}}$ and $c_{i,1} = \sigma_{i,1}/c_{1,1}$, for $i = 2, \ldots, m$, and we set

$$y_1 = u_1 = \frac{K_\nu^{(1)}\int\limits_{a_1}^{b_1} s(1+\frac{s^2}{\nu})^{-\frac{\nu+1}{2}}ds}{t_\nu(b_1)-t_\nu(a_1)}.$$

Given this (expected) value for $y_1$, the second integration variable is chosen by selecting a variable $i$ where

$$\min\limits_{2\leq i \leq m}\Big\{ t_{\nu+1}\Big(\sqrt{\frac... ...frac{a_i-c_{i,1}y_1}{\sqrt{\sigma_{i,i}-c_{i,1}^2}}\Big)\Big\}$$

is achieved. The integration limits, rows and columns of $\mbox{\boldmath$\Sigma$}$, and rows of $C$ for variables $2$ and $i$ are interchanged. Then the second column of $C$ is computed using $c_{2,2} = \sqrt{\sigma_{2,2}-c_{2,1}^2}$ and $c_{i,2} = (\sigma_{i,2}-c_{2,1}c_{i,1})/c_{2,2}$, for $i = 3, \ldots, m$, and we compute the expected value for $u_2$ using

$$u_2 = \frac{K_{\nu+1}^{(1)}\int\limits_{{\tilde a}_2}^{{\til... ...u+2}{2}}ds}{ t_{\nu+1}({\tilde b}_2)-t_{\nu+1}({\tilde a}_2)},$$

and we set $y_2 = u_1\sqrt{\frac{\nu+y_1^2}{\nu+1}}$. At stage $j$, given the expected values for $y_1, y_2, \ldots, y_{j-1}$, the $j$th integration variable is chosen by selecting a variable $i$, where

$$\min\limits_{j\leq i \leq m}\Big\{ t_{\nu+j-1}\Big(\sqrt{\fr... ...k} {\sqrt{\sigma_{i,i}-\sum_{k=1}^{j-1}c_{i,k}^2}}\Big) \Big\}$$

is achieved. The integration limits, rows and columns of $\mbox{\boldmath$\Sigma$}$, and rows of $C$ for variables $j$ and $i$ are interchanged. Then the $j$th column of $C$ is computed using $c_{j,j} = \sqrt{\sigma_{j,j}-\sum_{k=1}^{j-1}c_{j,k}^2}$ and $c_{i,j} = (\sigma_{i,j}-\sum_{k=1}^{j-1}c_{j,k}c_{i,k})/c_{j,j}$, for $i = j+1, \ldots, m$, and we let

$$u_j = \frac{K_{\nu+j-1}^{(1)}\int\limits_{{\tilde a}_j}^{{\t... ...{2}}ds}{ t_{\nu+j-1}({\tilde b}_j)-t_{\nu+j-1}({\tilde a}_j)}.$$

and set $y_j = u_j\sqrt{\frac{\nu+\sum_{i=1}^{j-1}y_i^2}{\nu+j-1}}$. The complete $m-1$ stage process has overall cost $O(m^3)$, which is not significant compared to the rest of the computation cost for the methods discussed here, and is therefore a relatively cheap preconditioning step that can be used with the algorithms. This sorting technique was included in the implementations for some of the algorithms used for the tests reported later in this paper.

A second set of MVT SV methods are based on the combination of equation (2) with MVN SV methods using

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath$\Sigma$}, \nu ) = ... ...\hat{b}_m(s,y_1,...,y_{m-1})} e^{-\frac{y_m^2}{2}}d{\bf y}ds,$$ (13)

with $\hat{a}_i(s,y_1....,y_{i-1})=\frac{s}{\sqrt{\nu}} (a_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i},$ and $\hat{b}_i(s,y_1....,y_{i-1})=\frac{s}{\sqrt{\nu}} (b_i-\sum_{j=1}^{i-1}c_{i,j}y_j)/c_{i,i}$. Then

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath $\Sigma$},\nu ) = ... ...s,{\bf w})-\hat{d}_{m}(s,{\bf w}))\int\limits_0^1 d{\bf w}ds,$$

with $\hat{z}_i=\hat{d}_i+(\hat{e}_i-\hat{d}_i)w_i$, and

$$\begin{eqnarray*} \hat{d}_i(s,{\bf w})&=& \Phi(\hat{a}_m(s,\Phi^{-1}(\hat{z}_1(w... ...hi^{-1}(\hat{z}_1(w_1)),...,\Phi^{-1}(\hat{z}_{i-1}(w_{i-1})))). \end{eqnarray*}$$

The last ${\bf w}$ component integral has value one, so the ${\bf T}$ integral is determined as an $m$-dimensional integral. A simple Monte-Carlo algorithm for this formulation of the MVT problem uses

$${\bf T}( {\bf a}, {\bf b},\mbox{\boldmath$\Sigma$},\nu ) \ap... ...rod_{i=1}^m(\hat{e}_{i}(s_k,{\bf w})-\hat{d}_{i}(s_k,{\bf w}))$$ (14)

with $w_{i,k} \sim Uniform(0,1)$, for $i=1, 2, \ldots, m-1$, and $s_k\sim \chi_\nu$. An inverse $\chi$ distribution function can be used to generate the required $s_k$'s from $Uniform(0,1)$ random numbers, via the formula $s_k=\chi^{-1}(w_{m,k})$, so each term in the Monte-Carlo sum for ${\bf T}$ requires $m Uniform(0,1)$ random numbers.

The Gibson, Glasbey and Elston technique can also be used as a preconditioning step for this method. In order to simplify our implementation of this preconditioning, we use $\sqrt{\nu}$ for the expected value for $s$, instead of the theoretical $\sqrt{2}\Gamma((\nu+1)/2)/\Gamma(\nu/2)$ (which approaches $\sqrt{\nu}$ for large $\nu$). This value for $s$ cancels the $\sqrt{\nu}$'s in the integration limits $\hat{{\bf a}}$ and $\hat{{\bf b}}$, so the preconditioning step for this method used the original Gibson, Glasbey and Elston technique with input ${\bf a}$, ${\bf b}$, and $\mbox{\boldmath$\Sigma$}$ (and $\mbox{\boldmath$\delta$}$ for NCMVT problems) and completes the preconditioning step as if the problem were an MVN problem.