Deák's Methods

These methods use a transformation to a spherical coordinate system. First, let ${\bf x}= C{\bf y}$, where $\Sigma = CC^t$ is the Cholesky decomposition of $\Sigma$. Then

$${\bf x}^t \Sigma^{-1} {\bf x}= {\bf y}^tC^{t}C^{-t}C^{-1}C{\bf y}= {\bf y}^t{\bf y}$$

so

$$P({\bf b}) = (2\pi)^{-\frac{m}{2}} \int_{-\infty < Cy \leq b} e^{-\frac{1}{2}{\bf y}^t{\bf y}} d{\bf y}.$$

Next let ${\bf y}= r{\bf z}$ with $\vert\vert{\bf z}\vert\vert \equiv \sqrt{{\bf z}^t{\bf z}} = 1$, so

$$\begin{eqnarray*}P({\bf b}) & = & (2\pi)^{-\frac{m}{2}} \int_{\vert\vert{\bf z}\... ...\int_{\vert\vert{\bf z}\vert\vert = 1}F(\rho({\bf z}))d{\bf z}, \end{eqnarray*}$$

where K_m is a normalization constant,

$$\rho({\bf z}) = \left\{ \begin{array}{cl} \begin{array}{c} \... ...{for all i} \\ \infty & \mbox{otherwise} \end{array}, \right.$$

with $w_i = \sum_{j=1}^ic_{i,j}z_j > 0$, and

$$F(\rho)= H_m\int_0^\rho r^{m-1}e^{-\frac{r^2}{2}}dr,$$

with H_m chosen so that $F(\infty) = 1$. $\rho({\bf z})$ is the distance, in the ${\bf z}$ direction, from the origin to the boundary of the integration region. The definition given here assumes b_i > 0 for all i, but this can been generalized (see Deák, 1986) to include negative b_i's.

Deák's methods all compute

$$P({\bf b}) \approx \frac{1}{N}\sum_{j=1}^N F(\rho({\bf z}_j))$$

using points $\{{\bf z}_j\}$ that are randomly chosen from the surface of m-sphere. His simplest method uses uniformly random $\{{\bf z}_j\}$ from the surface of m-sphere. Better methods improve the sampling by using antithetic variates. Let Z be an $m \times m$ random orthogonal matrix with columns $\{{\bf z}_j\}$ and let

$$S_n(Z) = \frac{1}{2^n{m \choose n}} \sum_{{\bf s}}\sum_{1 \le... ...ho\bigg(\frac{\sum_{l=1}^n s_l z_{j_l}}{\sqrt{n}}\bigg)\bigg),$$

where ${\bf s}= (s_1, s_2, ..., s_n) = (\pm 1, ..., \pm 1)$ and the outer sum is taken over the 2ⁿ possible sign combinations for the components of ${\bf s}$. The sample points used by S_n(Z) are very evenly spread over the surface of the unit m-sphere. Final estimates for $P({\bf b})$ are given by

$$\hat{P} = \frac{1}{M}\sum_{k=1}^M S_m(Z_k).$$

The standard error for the sample can be used to provide an error estimate.