# The Cumulative Distribution Unction

We have previously seen[1][2] how we can generalise normally distributed random variables to multiple dimensions by defining vectors with elements that are linear functions of independent standard normally distributed random variables, having means of zero and standard deviations of one, with
$\mathbf{Z}^\prime = \mathbf{L} \times \mathbf{Z} + \mathbf{\mu}$
where $$\mathbf{L}$$ is a constant matrix, $$\mathbf{Z}$$ is a vector whose elements are the independent standard normally distributed random variables and $$\mathbf{\mu}$$ is a constant vector.
So far we have derived and implemented the probability density function and the characteristic function of the multivariate normal distribution that governs such random vectors but have yet to do the same for its cumulative distribution function since it's a rather more difficult task and thus requires a dedicated treatment, which we shall have in this post.

Recall that the elements of $$\mathbf{Z}^\prime$$ are themselves normally distributed but, unlike the elements of $$\mathbf{Z}$$, are linearly dependent with a covariance matrix of
$\mathbf{\Sigma} = \mathbf{L} \times \mathbf{L}^\mathrm{T}$
where $$\mathbf{L}^\mathrm{T}$$ is the transpose of $$\mathbf{L}$$ in which its rows and columns are swapped, satisfying
$\Sigma_{ij} = \mathrm{E}\left[\left(Z^\prime_i - \mu_i\right) \times \left(Z^\prime_j - \mu_j\right)\right]$
where $$\mathrm{E}[X]$$ stands for the expected value of $$X$$. Compare this to a pair of independent normally distributed random variables $$Z_i$$ and $$Z_j$$, having means of $$\mu_i$$ and $$\mu_j$$ and standard deviations of $$\sigma_i$$ and $$\sigma_j$$ respectively, for which
\begin{align*} \mathrm{E}\left[\left(Z_i - \mu_i\right) \times \left(Z_i - \mu_i\right)\right] &= \sigma_i^2\\ \mathrm{E}\left[\left(Z_j - \mu_j\right) \times \left(Z_j - \mu_j\right)\right] &= \sigma_j^2\\ \mathrm{E}\left[\left(Z_i - \mu_i\right) \times \left(Z_j - \mu_j\right)\right] &= 0 \end{align*}
Note that we indicate that a random variable $$\mathbf{Z}$$ is multivariate normally distributed with a mean of $$\mathbf{\mu}$$ and a covariance matrix of $$\mathbf{\Sigma}$$ with
$\mathbf{Z} \sim \mathbf{N}(\mathbf{\mu}, \mathbf{\Sigma})$
where $$\sim$$ means drawn from and $$\mathbf{N}(\mathbf{\mu}, \mathbf{\Sigma})$$ represents the distribution, much as we do for a univariate normal distribution with a mean of $$\mu$$ and a standard deviation of $$\sigma$$ with
$Z \sim N(\mu, \sigma)$
Furthermore, the CDF of a multivariate distribution yields the probability that every element of an observation of a random variable $$\mathbf{X}$$ governed by that distribution is less than or equal to the corresponding element of its argument
$P_\mathbf{X}(\mathbf{x}) = \Pr\left(X_0 \leqslant x_0 \wedge X_1 \leqslant x_1 \wedge X_2 \leqslant x_2 \dots\right)$
where $$P_\mathbf{X}$$ is the CDF and $$\wedge$$ means and, generalising the definition for univariate distributions
$P_X(x) = \Pr\left(X \leqslant x\right)$

### The Independent Case

Now if the covariance matrix is diagonal, where the only non-zero elements are those having equal row and column indices, then the elements of the random variable will have the expectations given above for $$Z_i$$ and $$Z_j$$ and will consequently be independent of each other. Since the probability of independent events happening simultaneously equals the product of the probabilities of their each happening alone we have
$\Pr\left(Z_i \leqslant x_i \wedge Z_j \leqslant x_j\right) = \Pr\left(Z_i \leqslant x_i\right) \times \Pr\left(Z_j \leqslant x_j\right)$
and similarly, if we denote the multivariate normal CDF as $$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}$$ and the univariate normal CDF as $$\Phi_{\mu,\sigma}$$, for diagonal $$\mathbf{\Sigma}$$ we also have
$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) = \prod_i \Phi_{\mu_i,\sigma_i}(x_i)$
where $$\prod$$ is the product sign and $$\sigma_i$$ is the square root of $$\Sigma_{ii}$$.

Noting that we can express a univariate normal CDF with a mean of $$\mu$$ and a standard deviation of $$\sigma$$ in terms of the standard univariate normal CDF with
$\Phi_{\mu,\sigma}(x) = \Phi\left(\frac{x-\mu}{\sigma}\right)$
the uncorrCDF function given in listing 1 implements this formula for the uncorrelated CDF with the arguments mu, root and ucdf being an array of the values $$\mu_i$$, an array of the reciprocals of the values $$\sigma_i$$ and a standard univariate normal CDF respectively.

Listing 1: The uncorrCDF Helper Function
function uncorrCDF(mu, root, ucdf, x) {
var n = mu.length;
var c = 1;
var i, rt;

if(ak.type(x)!==ak.VECTOR_T) {
throw new Error('invalid argument in ak.multiNormalCDF');
}

if(x.dims()!==n) {
throw new Error('dimension mismatch in ak.multiNormalCDF');
}

for(i=0;i<n;++i) c *= ucdf((x.at(i)-mu[i])*root[i]);
return c;
}


Whilst this uncorrelated CDF has been trivial to derive and implement, the correlated case is unfortunately a much more difficult proposition.

### The Correlated Bivariate Case

We shall begin with the correlated two dimensional, or bivariate, normal CDF which we can define in terms of the bivariate normal PDF $$\phi_{\mathbf{\mu},\mathbf{\Sigma}}$$ with the double integral
$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) = \int_{-\infty}^{x_0} \int_{-\infty}^{x_1} \phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{z}) \; \mathrm{d}z_1 \; \mathrm{d}z_0$
much as we can define the univariate CDF as the integral of the univariate PDF $$\phi_{\mu,\sigma}$$
$\Phi_{\mu,\sigma}(x) = \int_{-\infty}^{x} \phi_{\mu,\sigma}(z) \; \mathrm{d}z$
Now, given this definition, it's rather tempting to try using a nested numerical integrator[3], such as our ak.nestedIntegral of our ak.rombergIntegral, on the bivariate PDF. One thing standing in our way is the fact that the lower bounds of integration equal minus infinity and so we won't be able to divide up the intervals of integration into sub-intervals of equal width.
Fortunately we can get around the problem with integration by substitution in which we change the variable over which we're integrating to something more manageable, specifically the CDFs of the elements of $$\mathbf{x}$$ when taken independently, $$\Phi_{\mu_0,\sigma_0}(x_0)$$ and $$\Phi_{\mu_1,\sigma_1}(x_1)$$ where $$\sigma_0$$ and $$\sigma_1$$ are the square roots of $$\Sigma_{00}$$ and $$\Sigma_{11}$$ respectively, yielding
$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) = \int_0^{\Phi_{\mu_0,\sigma_0}(x_0)} \int_0^{\Phi_{\mu_1,\sigma_1}(x_1)} \phi_{\mathbf{\mu},\mathbf{\Sigma}}\left(\begin{pmatrix}z_0\\z_1\end{pmatrix}\right) \; \frac{\mathrm{d}z_1}{\mathrm{d}u_1} \mathrm{d}u_1 \; \frac{\mathrm{d}z_0}{\mathrm{d}u_0}\mathrm{d}u_0$
where
\begin{align*} u_0 &= \Phi_{\mu_0,\sigma_0}\left(z_0\right)\\ u_1 &= \Phi_{\mu_1,\sigma_1}\left(z_1\right) \end{align*}
This means that
\begin{align*} z_0 &= \Phi^{-1}_{\mu_0,\sigma_0}\left(u_0\right)\\ z_1 &= \Phi^{-1}_{\mu_1,\sigma_1}\left(u_1\right) \end{align*}
and
\begin{align*} \frac{\mathrm{d}z_0}{\mathrm{d}u_0} &= \frac{1}{\frac{\mathrm{d}u_0}{\mathrm{d}z_0}} = \frac{1}{\phi_{\mu_0,\sigma_0}(z_0)} = \frac{1}{\phi_{\mu_0,\sigma_0}\left(\Phi^{-1}_{\mu_0,\sigma_0}\left(u_0\right)\right)}\\ \frac{\mathrm{d}z_1}{\mathrm{d}u_1} &= \frac{1}{\frac{\mathrm{d}u_1}{\mathrm{d}z_1}} = \frac{1}{\phi_{\mu_1,\sigma_1}(z_1)} = \frac{1}{\phi_{\mu_1,\sigma_1}\left(\Phi^{-1}_{\mu_1,\sigma_1}\left(u_1\right)\right)} \end{align*}
and so
$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) = \int_0^{\Phi_{\mu_0,\sigma_0}(x_0)} \int_0^{\Phi_{\mu_1,\sigma_1}(x_1)} \frac{\phi_{\mathbf{\mu},\mathbf{\Sigma}}\left(\begin{pmatrix}\Phi^{-1}_{\mu_0,\sigma_0}(u_0)\\\Phi^{-1}_{\mu_1,\sigma_1}(u_1)\end{pmatrix}\right)}{\phi_{\mu_0,\sigma_0}\left(\Phi^{-1}_{\mu_0,\sigma_0}(u_0)\right) \times \phi_{\mu_1,\sigma_1}\left(\Phi^{-1}_{\mu_1,\sigma_1}(u_1)\right)} \; \mathrm{d}u_1 \; \mathrm{d}u_0$
Now we have to be a little bit careful when both the numerator and the denominator of the integrand are equal to zero, which will happen when either $$u_0$$ or $$u_1$$ are equal to either zero or one, since it will evaluate to NaN rather than zero as the original integrand would have. To handle this we can set the integrand to zero whenever its denominator equals zero, as demonstrated by program 1 which compares the result of this CDF integral to the proportion of a sample of the appropriate bivariate normal random variable that have elements that are all less than or equal to its argument.

Program 1: Integrating The Bivariate PDF

That the results of the integral are roughly equal to those of the sample is reassuring, but unfortunately it's not very efficient; if you try decreasing the numerical integrator's convergence threshold eps to improve the accuracy, say by dividing it by ten, you'll find that it will fail to converge.

Fortunately, we can do a lot better.

### Flying The Nest

As unlikely as it might seem, it is possible to approximate a correlated bivariate normal CDF with a one dimensional numerical integral. In particular a correlated standard bivariate normal CDF, for which both elements have a mean of zero and a standard deviation of one, having a correlation of $$\rho$$ can be expressed as
$\Phi_\rho\left(x_0, x_1\right) = \int_0^{\Phi\left(x_1\right)} \Phi\left(\frac{x_0-\rho\Phi^{-1}\left(u_1\right)}{\sqrt{1-\rho^2}}\right) \; \mathrm{d}u_1$
where $$\Phi$$ is the standard univariate normal CDF having a mean of zero and a standard deviation of one and $$\Phi^{-1}$$ is its inverse, as proven in derivation 1.

Derivation 1: The Correlated Standard Bivariate Normal CDF
For a correlation $$\rho$$ the correlated standard bivariate normal covariance matrix equals
$\Sigma = \begin{pmatrix} 1 & \rho\\ \rho & 1 \end{pmatrix}$
which has a determinant of
$|\Sigma| = \begin{vmatrix} 1 & \rho\\ \rho & 1 \end{vmatrix} = 1 - \rho^2$
and an inverse of
$\Sigma^{-1} = \frac{1}{|\Sigma|} \begin{pmatrix} 1 & -\rho\\ -\rho & 1 \end{pmatrix}$
and so the correlated standard bivariate PDF is given by
\begin{align*} \phi_\rho\left(z_0, z_1\right) &= \frac{1}{2\pi\sqrt{1-\rho^2}} e^{-\tfrac12 \tfrac{z_0^2 - 2 \rho z_0z_1 + z_1^2}{1-\rho^2}} = \frac{1}{2\pi\sqrt{1-\rho^2}} e^{-\tfrac12 \tfrac{\left(z_0 - \rho z_1\right)^2 - \rho^2 z_1^2 + z_1^2}{1-\rho^2}}\\ &= \frac{1}{2\pi\sqrt{1-\rho^2}} e^{-\tfrac12 \tfrac{\left(z_0 - \rho z_1\right)^2 + \left(1 - \rho^2\right) z_1^2}{1-\rho^2}} = \frac{1}{2\pi\sqrt{1-\rho^2}} e^{-\tfrac12 \tfrac{\left(z_0 - \rho z_1\right)^2}{1-\rho^2}} e^{-\tfrac12 z_1^2}\\ &= \frac{1}{\sqrt{2\pi}}e^{-\tfrac12 z_1^2} \times \frac{1}{\sqrt{2\pi}\sqrt{1-\rho^2}} e^{-\tfrac12 \tfrac{\left(z_0 - \rho z_1\right)^2}{1-\rho^2}} = \phi\left(z_1\right) \times \phi_{\rho z_1,\sqrt{1-\rho^2}}\left(z_0\right) \end{align*}
where $$\phi$$ is the standard univariate normal PDF.
The correlated standard bivariate CDF is consequently given by
\begin{align*} \Phi_\rho\left(x_0, x_1\right) &= \int_{-\infty}^{x_1} \int_{-\infty}^{x_0} \phi_\rho\left(z_0, z_1\right) \mathrm{d}z_0 \; \mathrm{d}z_1 = \int_{-\infty}^{x_1} \int_{-\infty}^{x_0} \phi\left(z_1\right) \times \phi_{\rho z_1,\sqrt{1-\rho^2}}\left(z_0\right) \mathrm{d}z_0 \; \mathrm{d}z_1\\ &= \int_{-\infty}^{x_1} \phi\left(z_1\right) \times \left(\int_{-\infty}^{x_0} \phi_{\rho z_1,\sqrt{1-\rho^2}}\left(z_0\right) \mathrm{d}z_0 \right) \; \mathrm{d}z_1\\ &= \int_{-\infty}^{x_1} \phi\left(z_1\right) \times \Phi_{\rho z_1,\sqrt{1-\rho^2}}\left(x_0\right) = \int_{-\infty}^{x_1} \phi\left(z_1\right) \times \Phi\left(\frac{x_0 - \rho z_1}{\sqrt{1-\rho^2}}\right) \mathrm{d}z_1 \end{align*}
Finally, by substituting $$u_1 = \Phi\left(z_1\right)$$ we have
\begin{align*} \Phi_\rho\left(x_0, x_1\right) &= \int_0^{\Phi\left(x_1\right)} \phi\left(\Phi^{-1}\left(u_1\right)\right) \times \Phi\left(\frac{x_0-\rho\Phi^{-1}\left(u_1\right)}{\sqrt{1-\rho^2}}\right) \; \frac{\mathrm{d}u_1}{\phi\left(\Phi^{-1}\left(u_1\right)\right)}\\ &= \int_0^{\Phi\left(x_1\right)} \Phi\left(\frac{x_0-\rho\Phi^{-1}\left(u_1\right)}{\sqrt{1-\rho^2}}\right) \; \mathrm{d}u_1 \end{align*}

We can trivially generalise this to the non-standard case by noting that, for a two dimensional covariance matrix $$\mathbf{\Sigma}$$
$\rho = \frac{\Sigma_{0,1}}{\sqrt{\Sigma_{0,0}} \times \sqrt{\Sigma_{1,1}}}$
and so
$\Phi_{\mathbf{\mu},\mathbf{\Sigma}}\left(\mathbf{x}\right) = \Phi_\rho\left(\frac{x_0-\mu_0}{\sqrt{\Sigma_{0,0}}}, \frac{x_1-\mu_1}{\sqrt{\Sigma_{1,1}}}\right)$
since this transforms the elements of $$\mathbf{x}$$ into standard normal random variables whilst preserving their correlation.

Program 2 replaces the two dimensional integral of the PDF in program 1 with this one dimensional integral to demonstrate that it works.

Program 2: Integrating The Univariate CDF

Well that's certainly a lot more efficient! In fact, you can divide the convergence threshold by ten thousand and it will still reliably converge.

Whilst this demonstrates the principal of reducing the dimension of the integral, it is not the most numerically stable way of doing so. Instead we shall turn to Drezner and Wesolowsky's formula[4]
$\Phi_\rho\left(x_0, x_1\right) = \Phi\left(x_0\right) \times \Phi\left(x_1\right) + \frac{1}{2\pi} \int_0^\rho \frac{1}{\sqrt{1-r^2}} e^{-\tfrac{x_0^2 - 2 x_0 x_1 r + x_1^2}{2\left(1-r^2\right)}}\mathrm{d}r$
as transformed by Genz[5] by substituting $$\sin(\theta)$$ for $$r$$ to yield
$\Phi_\rho\left(x_0, x_1\right) = \Phi\left(x_0\right) \times \Phi\left(x_1\right) + \frac{1}{2\pi} \int_0^{\sin^{-1}(\rho)} e^{-\tfrac{x_0^2 - 2 x_0 x_1 \sin(\theta) + x_1^2}{2\cos^2(\theta)}}\mathrm{d}\theta$
as shown by derivation 2.

Derivation 2: Genz's Substitution
Firstly, we have
\begin{align*} r &= \sin(\theta)\\ \theta &= \sin^{-1}(r)\\ \frac{\mathrm{d}r}{\mathrm{d}\theta} &= \cos(\theta) \end{align*}
Next, we have
\begin{align*} &\cos^2(\theta) + \sin^2(\theta) = 1\\ &\cos^2(\theta) = 1 - \sin^2(\theta) = 1 - r^2 \end{align*}
Changing the variable of integration to $$\theta$$ therefore yields
\begin{align*} \int_0^\rho \frac{1}{\sqrt{1-r^2}} e^{-\tfrac{x_0^2 - 2 x_0 x_1 r + x_1^2}{2\left(1-r^2\right)}}\mathrm{d}r &=\int_{\sin^{-1}(0)}^{\sin^{-1}(\rho)} \frac{1}{\sqrt{\cos^2(\theta)}} e^{-\tfrac{x_0^2 - 2 x_0 x_1 \sin(\theta) + x_1^2}{2\cos^2(\theta)}} \frac{\mathrm{d}r}{\mathrm{d}\theta}\;\mathrm{d}\theta\\ &=\int_0^{\sin^{-1}(\rho)} \frac{1}{\cos(\theta)} e^{-\tfrac{x_0^2 - 2 x_0 x_1 \sin(\theta) + x_1^2}{2\cos^2(\theta)}} \cos(\theta)\;\mathrm{d}\theta\\ &=\int_0^{\sin^{-1}(\rho)} e^{-\tfrac{x_0^2 - 2 x_0 x_1 \sin(\theta) + x_1^2}{2\cos^2(\theta)}} \mathrm{d}\theta \end{align*}

The corrStd2dCDF helper function given in listing 2 implements this formula for the correlated standard bivariate normal distribution, with the arguments r, x0, x1, c0 and c1 representing the correlation, the first and second elements of the argument and their standard univariate normal CDF results respectively, and the threshold and steps arguments controlling the behaviour of the ak.rombergIntegral integrator.

Listing 2: The corrStd2dCDF Helper Function
function corrStd2dCDF(r, x0, x1, c0, c1, threshold, steps) {
var c = 0;
var p, s, f;

if(c0!==0 && c0!==1 && c1!==0 && c1!==1) {
p = x0*x1;
s = (x0*x0+x1*x1)/2;

f = function(t) {
var st = Math.sin(t);
var ct = Math.cos(t);
return Math.exp(-(s-p*st)/(ct*ct));
};

c = ak.rombergIntegral(f, threshold, steps)(0, Math.asin(r))/(2*ak.PI);
}

return Math.max(0, Math.min(1, c0*c1 + c));
}


Now it's trivially the case that if the value of the CDF for either x0 or x1 is zero then their bivariate CDF must also equal zero. Furthermore, if either equals one then the value of the bivariate CDF must equal the value of the univariate CDF of the other and we handle both cases by simply not including the integral in the result.
Furthermore, just in case numerical error pushes the result outside of its valid range of zero to one, we clamp it to that range in the return statement.

The corr2dCDF helper function given in listing 3 generalises this to the non-standard bivariate CDF with exactly the same transformation of the arguments that we used in the double integral of the PDF.

Listing 3: The corr2dCDF Helper Function
function corr2dCDF(mu, s0, s1, r, x, cdf, threshold, steps) {
var z0, z1, c0, c1;

if(ak.type(x)!==ak.VECTOR_T) {
throw new Error('invalid argument in ak.multiNormalCDF');
}

if(x.dims()!==2) {
throw new Error('dimension mismatch in ak.multiNormalCDF');
}

x0 = (x.at(0)-mu[0])/s0;
x1 = (x.at(1)-mu[1])/s1;
if(isNaN(x0) || isNaN(x1)) return ak.NaN;

c0 = cdf(x0);
c1 = cdf(x1);
return corrStd2dCDF(r, x0, x1, c0, c1, threshold, steps);
}


Here, mu is an ak.vector representing the mean, s0 and s1 are the standard deviations of the elements of the argument, r is the correlation, x is the ak.vector argument, cdf is a standard univariate normal CDF and threshold and steps are the control parameters for the numerical integrator.
Now this is where we check that the argument is, in fact, a two dimensional vector and handle NaNs before forwarding to the corrStd2dCDF function to do the actual work.

Finally, the makeCorr2dCDF function given in listing 4 computes the standard deviations s0 and s1 and the correlation r from the ak.matrix covariance matrix sigma before returning a function that forwards to corr2dCDF with its argument x.

Listing 4: The makeCorr2dCDF Helper Function
function makeCorr2dCDF(mu, sigma, threshold, steps) {
var cdf = ak.normalCDF();
var s0 = Math.sqrt(sigma[0][0]);
var s1 = Math.sqrt(sigma[1][1]);
var r = sigma[0][1]/(s0*s1);

return function(x) {
return corr2dCDF(mu, s0, s1, r, x, cdf, threshold, steps);
};
}


### The Correlated Trivariate Case

We can pull a similar trick with the three dimensional correlated standard trivariate normal distribution, reducing the integral from three to two dimensions with Gupta's[6] formula
$\Phi_{\rho_{01}, \rho_{02}, \rho_{12}}\left(x_0, x_1, x_2\right) = \int_{-\infty}^{x_0} \Phi_\rho \left(\frac{x_1 - \rho_{01}z}{\sqrt{1-\rho_{01}^2}}, \frac{x_2 - \rho_{02}z}{\sqrt{1-\rho_{02}^2}}\right) \phi(z) \; \mathrm{d}z$
where $$\rho_{01}$$, $$\rho_{02}$$ and $$\rho_{12}$$ are the correlations between the first and second, the first and third and the second and third elements of the random variable respectively and
$\rho = \frac{\rho_{12} - \rho_{01}\rho_{02}}{\sqrt{\left(1-\rho_{01}^2\right)\left(1-\rho_{02}^2\right)}}$
which follows from much the same argument as we made in derivation 1.
As was the case there, we can't numerically integrate from minus infinity and so we must again change the variable of integration to the CDF with $$u = \Phi(z)$$, yielding
$\Phi_{\rho_{01}, \rho_{02}, \rho_{12}}\left(x_0, x_1, x_2\right) = \int_0^{\Phi\left(x_0\right)} \Phi_\rho \left(\frac{x_1 - \rho_{01}\Phi^{-1}\left(u\right)}{\sqrt{1-\rho_{01}^2}}, \frac{x_2 - \rho_{02}\Phi^{-1}\left(u\right)}{\sqrt{1-\rho_{02}^2}}\right) \; \mathrm{d}u$
This integral is implemented by the corrStd3dCDF helper function given in listing 5.

Listing 5: The corrStd3dCDF Helper Function
function corrStd3dCDF(r01, r02, rt01, rt02, rho, x0, x1, x2,
cdf, inv, threshold, steps) {
var f = function(u) {
var z = inv(u);
var z0 = (x1-r01*z)*rt01;
var z1 = (x2-r02*z)*rt02;

if(isNaN(z0) || isNaN(z1)) return ak.NaN;
return corrStd2dCDF(rho, z0, z1, cdf(z0), cdf(z1), threshold, steps);
};
var c = ak.rombergIntegral(f, threshold, steps)(0, cdf(x0));
return Math.max(0, Math.min(1, c));
}


Here r01 and r02 represent $$\rho_{01}$$ and $$\rho_{02}$$ respectively, rt01 and rt02 represent the reciprocals of the square roots of $$1-\rho_{01}^2$$ and $$1-\rho_{02}^2$$ respectively, rho represents $$\rho$$ and x0, x1 and x2 the three arguments $$x_0$$, $$x_1$$ and $$x_2$$. The cdf and inv arguments are the standard univariate normal CDF and its inverse and threshold and steps are the control parameters for the numerical integrator.
Note that we use the corrStd2dCDF helper function directly to compute the correlated standard bivariate normal CDF in the integrand f.

The corr3dCDF helper function given in listing 6 converts the vector argument of a general trivariate normal CDF to the standard arguments expected by corrStd3dCDF.

Listing 6: The corr3dCDF Helper Function
function corr3dCDF(mu, s0, s1, s2, r01, r02, r12, rt01, rt02, rho,
x, cdf, inv, threshold, steps) {
var x0, x1, x2, c0, c1, c2;

if(ak.type(x)!==ak.VECTOR_T) {
throw new Error('invalid argument in ak.multiNormalCDF');
}

if(x.dims()!==3) {
throw new Error('dimension mismatch in ak.multiNormalCDF');
}

x0 = (x.at(0)-mu[0])/s0;
x1 = (x.at(1)-mu[1])/s1;
x2 = (x.at(2)-mu[2])/s2;
if(isNaN(x0) || isNaN(x1) || isNaN(x2)) return ak.NaN;

c0 = cdf(x0);
c1 = cdf(x1);
c2 = cdf(x2);
if(c0===0 || c1===0 || c2===0) return 0;
if(c0===1) return corrStd2dCDF(r12, x1, x2, c1, c2, threshold, steps);
if(c1===1) return corrStd2dCDF(r02, x0, x2, c0, c2, threshold, steps);
if(c2===1) return corrStd2dCDF(r01, x0, x1, c0, c1, threshold, steps);

return corrStd3dCDF(r01, r02, rt01, rt02, rho,
x0, x1, x2, cdf, inv, threshold, steps);
}


Here we also have the vector mean passed in mu, the standard deviations in s0, s1 and s2, the correlation between the second and third arguments in r12 and the vector argument in x.
After checking that x is a three dimensional vector we transform its elements into their standard equivalents x0, x1 and x2 and store the values of the CDF for them in c0, c1 and c2 respectively.
Next, if any of those values equal zero the trivariate CDF must also be zero and if any of them equal one then it must equal the bivariate CDF for the other two and so we handle these corner cases before forwarding to our standard trivariate helper function.

Finally, the makeCorr3dCDF given in listing 7 creates the standard univariate CDF and its inverse, converts the covariance matrix sigma to the standard deviations and correlations expected by corr3dCDF and computes its rt01, rt02 and rho arguments.

Listing 7: The makeCorr3dCDF Helper Function
function makeCorr3dCDF(mu, sigma, threshold, steps) {
var cdf = ak.normalCDF();
var inv = ak.normalInvCDF();
var s0 = Math.sqrt(sigma[0][0]);
var s1 = Math.sqrt(sigma[1][1]);
var s2 = Math.sqrt(sigma[2][2]);
var r01 = sigma[0][1]/(s0*s1);
var r02 = sigma[0][2]/(s0*s2);
var r12 = sigma[1][2]/(s1*s2);
var rt01 = 1/Math.sqrt(1-r01*r01);
var rt02 = 1/Math.sqrt(1-r02*r02);
var rho = (r12-r01*r02)*rt01*rt02;

if(ak.nativeType(threshold)===ak.UNDEFINED_T) {
threshold = Math.pow(ak.EPSILON, 0.375);
}
else {
threshold = Math.sqrt(Math.abs(threshold));
}

if(ak.nativeType(steps)===ak.UNDEFINED_T) steps = 12;
else steps = ak.ceil(Math.abs(steps)*2/3);

return function(x) {
return corr3dCDF(mu, s0, s1, s2, r01, r02, r12, rt01, rt02, rho,
x, cdf, inv, threshold, steps);
};
}


Now if the threshold and steps numerical integrator parameters are undefined they are set to the square root and two thirds of their defaults for univariate integration respectively to reflect the fact that nested integration is inherently less accurate. Similarly, if they are defined, they are replaced by the square root and two thirds of their given values respectively.
This done, we simply return a function that forwards to the corr3dCDF helper function with its argument x in addition to those that we have already initialised.

### The Correlated Multivariate Case

For more than two dimensions, nested numerical integration starts to get rather impractical and so we shall instead turn to our quasi random integrator[7] ak.quasiRandomIntegral. As before, we can get away with integrating over one fewer dimension than the multivariate normal distribution has itself, this time with the formula for those having a mean of zero described by Genz and Bretz[8]
\begin{align*} x^\prime_i\left(y_0,\dots,y_{i-1}\right) &= \frac{x_i - \sum_{j=0}^{i-1} l_{ij} y_j}{l_{ii}}\\ e_0 &= \Phi\left(\frac{x_0}{l_{00}}\right)\\ e_i\left(u_0,\dots,u_{i-1}\right) &= \Phi\left(x^\prime_i\left(\Phi^{-1}\left(e_0 u_0\right),\dots,\Phi^{-1}\left(e_{i-1} u_{i-1}\right)\right)\right)\\ \Phi_{\mathbf{0},\mathbf{\Sigma}}(x) &= e_0 \int_0^1 e_1\left(u_0\right) \int_0^1 e_2\left(u_0, u_1\right) \dots \int_0^1 e_{n-1}\left(u_0, u_1, \dots, u_{n-2}\right) \; \mathrm{d}u_{n-2} \dots \mathrm{d}u_1 \mathrm{d}u_0 \end{align*}
where $$\sum$$ is the summation sign and $$l_{ij}$$ are the elements of the lower triangular root $$\mathbf{L}$$ of the $$n$$ dimensional covariance matrix $$\mathbf{\Sigma}$$, as shown in excruciating detail by derivation 3.

Derivation 3: The Multivariate Correlated CDF
Firstly, we have
\begin{align*} \Phi_{\mathbf{0},\mathbf{\Sigma}} \left(\mathbf{x}\right) &= \int_{-\mathbf{\infty}}^{\mathbf{x}} \phi_{\mathbf{0},\mathbf{\Sigma}}\left(\mathbf{z}\right) \; \mathrm{d}\mathbf{z} = \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{\Sigma}|^\frac12} e^{-\frac12 \mathbf{z} \mathbf{\Sigma}^{-1} \mathbf{z}} \; \mathrm{d}\mathbf{z} = \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \mathbf{z} \left(\mathbf{L} \mathbf{L}^\mathrm{T}\right)^{-1}\mathbf{z}} \; \mathrm{d}\mathbf{z}\\ &= \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \mathbf{z} \left(\mathbf{L}^{\mathrm{T}^{-1}} \mathbf{L}^{-1}\right)\mathbf{z}} \; \mathrm{d}\mathbf{z} = \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \mathbf{z} \left(\mathbf{L}^{\mathrm{T}^{-1}} \mathbf{L}^{-1}\right)\mathbf{z}} \; \mathrm{d}\mathbf{z}\\ &= \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \left(\mathbf{z} \mathbf{L}^{\mathrm{T}^{-1}}\right) \left(\mathbf{L}^{-1}\mathbf{z}\right)} \; \mathrm{d}\mathbf{z} = \int_{-\mathbf{\infty}}^{\mathbf{x}} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \left(\mathbf{L}^{-1}\mathbf{z}\right) \left(\mathbf{L}^{-1}\mathbf{z}\right)} \; \mathrm{d}\mathbf{z} \end{align*}
Next, let us define $$\mathbf{z} = \mathbf{L} \mathbf{z}^\prime$$ so that
$z_i = \sum_{j=0}^i l_{ij} z^\prime_j \implies \frac{\partial z_i}{\partial z^\prime_j} = l_{ij} \implies \left|\frac{\partial\mathbf{z}}{\partial\mathbf{z}^\prime}\right| = |\mathbf{L}|$
where $$\implies$$ means implies and $$\left|\frac{\partial\mathbf{z}}{\partial\mathbf{z}^\prime}\right|$$ is the determinant of the Jacobian matrix of partial derivatives of $$\mathbf{z}$$ with respect to the elements of $$\mathbf{z}^\prime$$, and
$z_i \leqslant x_i \implies \sum_{j=0}^i l_{ij} z_j^\prime \leqslant x_i \implies l_{ii} z^\prime_i \leqslant x_i - \sum_{j=0}^{i-1} l_{ij} z^\prime_j \implies z^\prime_i \leqslant \frac{x_i - \sum_{j=0}^{i-1} l_{ij} z^\prime_j}{l_{ii}} = x^\prime_i$
Using integration by substitution to change the variable of integration to $$\mathbf{z}^\prime$$ consequently yields $$\mathrm{d}\mathbf{z} = |\mathbf{L}| \mathrm{d}\mathbf{z}^\prime$$ and so
\begin{align*} \Phi_{\mathbf{0},\mathbf{\Sigma}} \left(\mathbf{x}\right) &= \int_{-\mathbf{\infty}}^{\mathbf{x}^\prime} \phi_{\mathbf{0},\mathbf{\Sigma}}\left(\mathbf{L} \mathbf{z}^\prime\right) \; |\mathbf{L}|\mathrm{d}\mathbf{z}^\prime = \int_{-\mathbf{\infty}}^{\mathbf{x}^\prime} \frac{1}{(2\pi)^{\frac{n}2} |\mathbf{L}|} e^{-\frac12 \left(\mathbf{L}^{-1}\mathbf{L}\mathbf{z}^\prime\right) \left(\mathbf{L}^{-1}\mathbf{L}\mathbf{z}^\prime\right)} \; |\mathbf{L}|\mathrm{d}\mathbf{z}^\prime\\ &= \int_{-\mathbf{\infty}}^{\mathbf{x}^\prime} \frac{1}{(2\pi)^{\frac{n}2}} e^{-\frac12 \mathbf{z}^\prime \times \mathbf{z}^\prime} \; \mathrm{d}\mathbf{z}^\prime = \int_{-\mathbf{\infty}}^{\mathbf{x}^\prime} \phi_{\mathbf{0},\mathbf{I}}(\mathbf{z}^\prime) \; \mathrm{d}\mathbf{z}^\prime = \int_{-\mathbf{\infty}}^{\mathbf{x}^\prime} \prod_{i=0}^{n-1}\phi(z^\prime_i) \; \mathrm{d}\mathbf{z}^\prime\\ &= \int_{-\infty}^{x^\prime_0} \phi\left(z^\prime_0\right) \int_{-\infty}^{x^\prime_1} \phi\left(z^\prime_1\right) \dots \int_{-\infty}^{x^\prime_{n-1}} \phi\left(z^\prime_{n-1}\right) \; \mathrm{d}z^\prime_{n-1} \dots \mathrm{d}z^\prime_1 \mathrm{d}z^\prime_0 \end{align*}
Now if we define $$u^\prime_i = \Phi\left(z^\prime_i\right)$$ we have $$z^\prime_i = \Phi^{-1}\left(u^\prime_i\right)$$ and
$\frac{\mathrm{d}u^\prime_i}{\mathrm{d}z^\prime_i} = \phi\left(z^\prime_i\right) \implies \mathrm{d}z^\prime_i = \frac{\mathrm{d}u^\prime_i}{\phi\left(z^\prime_i\right)}$
and substituting the $$u^\prime_i$$ for the $$z^\prime_i$$ in the integrals yields
\begin{align*} \Phi_{\mathbf{0},\mathbf{\Sigma}} \left(\mathbf{x}\right) &= \int_{\Phi(-\infty)}^{\Phi\left(x^\prime_0\right)} \phi\left(z^\prime_0\right) \int_{\Phi(-\infty)}^{\Phi\left(x^\prime_1\right)} \phi\left(z^\prime_1\right) \dots \int_{\Phi(-\infty)}^{\Phi\left(x^\prime_{n-1}\right)} \phi\left(z^\prime_{n-1}\right) \; \frac{\mathrm{d}u^\prime_{n-1}}{\phi\left(z^\prime_{n-1}\right)} \dots \frac{\mathrm{d}u^\prime_1}{\phi\left(z^\prime_1\right)} \frac{\mathrm{d}u^\prime_0}{\phi\left(z^\prime_0\right)}\\ &= \int_0^{\Phi\left(x^\prime_0\right)} \int_0^{\Phi\left(x^\prime_1\right)} \dots \int_0^{\Phi\left(x^\prime_{n-1}\right)} \mathrm{d}u^\prime_{n-1} \dots \mathrm{d}u^\prime_1 \mathrm{d}u^\prime_0 \end{align*}
Further defining $$u^\prime_i = \Phi\left(x^\prime_i\right) u_i$$ yields $$z^\prime_i = \Phi^{-1}\left(\Phi\left(x^\prime_i\right) u_i\right)$$ and
$u_i = \frac{u^\prime_i}{\Phi\left(x^\prime_i\right)} \implies \frac{\mathrm{d}u_i}{\mathrm{d}u^\prime_i} = \frac{1}{\Phi\left(x^\prime_i\right)} \implies \mathrm{d}u^\prime_i = \Phi\left(x^\prime_i\right)\mathrm{d}u_i$
and substituting the $$u_i$$ for the $$u^\prime_i$$ gives us
\begin{align*} \Phi_{\mathbf{0},\mathbf{\Sigma}} \left(\mathbf{x}\right) &= \int_0^1 \int_0^1 \dots \int_0^1 \Phi\left(x^\prime_{n-1}\right) \dots \Phi\left(x^\prime_1\right) \Phi\left(x^\prime_0\right) \mathrm{d}u_{n-1} \dots \mathrm{d}u_1 \mathrm{d}u_0 \end{align*}
Since each $$x^\prime_i$$ only depends upon $$z^\prime_j$$, and hence $$u_j$$, if $$j$$ is strictly less than $$i$$, we have
\begin{align*} \Phi_{\mathbf{0},\mathbf{\Sigma}} \left(\mathbf{x}\right) &= \Phi\left(x^\prime_0\right) \int_0^1 \Phi\left(x^\prime_1\right) \int_0^1 \Phi\left(x^\prime_2\right) \dots \int_0^1 \Phi\left(x^\prime_{n-1}\right) \int_0^1 \mathrm{d}u_{n-1} \mathrm{d}u_{n-2} \dots \mathrm{d}u_1 \mathrm{d}u_0\\ &= \Phi\left(x^\prime_0\right) \int_0^1 \Phi\left(x^\prime_1\right) \int_0^1 \Phi\left(x^\prime_2\right) \dots \int_0^1 \Phi\left(x^\prime_{n-2}\right) \mathrm{d}u_{n-2} \dots \mathrm{d}u_1 \mathrm{d}u_0 \end{align*}
Finally, we define $$e_i = \Phi\left(x^\prime_i\right)$$ so that
$e_i u_i = \Phi\left(x^\prime_i\right) u_i = u^\prime_i = \Phi\left(z^\prime_i\right) \implies z^\prime_i = \Phi^{-1}\left(e_i u_i\right) \implies x^\prime_i = \frac{x_i - \sum_{j=0}^{i-1} l_{ij} \Phi^{-1}\left(e_i u_i\right)}{l_{ii}}$
from which the formula trivially follows.

Since we shall be approximating this with a quasi random sample we don't care about the way the integrals nest and so we shall just need the integrand
$f = \prod_{i=0}^{n-1} e_i$
as implemented in listing 8.

Listing 8: The corrNdIntegrand Helper Function
function corrNdIntegrand(root, x, cdf, inv, e0, y) {
var n = x.length;

return function(u) {
var em = e0;
var fm = e0;
var m, ly, i;

for(m=1;m<n;++m) {
y[m-1] = inv(u.at(m-1)*em);

ly = 0;
for(i=0;i<m;++i) ly += root[m][i]*y[i];

em = cdf((x[m]-ly)*root[m][m]);
fm *= em;
}
return fm;
};
}


Here, root is an array of arrays of the elements of $$\mathbf{L}$$ with those on the leading diagonal, having row indices equal to their column indices, being their reciprocals. Since $$x_0$$ doesn't depend upon any of the variables over which we're integrating we pass it in as an argument to save recalculating it each time we evaluate the integrand. Finally, x is an array of the elements of the argument of the CDF, cdf is the standard normal CDF, inv is its inverse and y is a buffer for the $$y_i$$ values.

Now, we'll have a bit of a problem if any of the elements of x are infinite since the sum ly might end up subtracting infinity from infinity to yield NaN. It is trivially the case that if any of them are negative infinity then the CDF must equal zero but if any of them are positive infinity then we must be a little more careful. Specifically, we must reduce the integral to one over the finite elements, much as we did for the two and three dimensional cases, as done by the corrNdMarginalCDF helper function given in listing 9.

Listing 9: The corrNdMarginalCDF Helper Function
function corrNdMarginalCDF(sigma, devs, x, cdf, threshold, steps, qrnd, prnd) {
var n = x.length;
var c = 0;
var f = [];
var i, j, k, s, y, mcdf;

for(i=0;i<n;++i) if(cdf(x[i]*devs[i])<1) f[c++] = i;
if(c===0) return 1;
if(c===1) return cdf(x[f[0]]*devs[f[0]]);

s = new Array(c);
y = new Array(c);

for(i=0;i<c;++i) {
k = f[i];
s[i] = new Array(c);
y[i] = x[k];
for(j=0;j<c;++j) s[i][j] = sigma[k][f[j]];
}
mcdf = ak.multiNormalCDF(ak.matrix(s), threshold, steps, qrnd, prnd);
return mcdf(ak.vector(y));
}


The sigma argument is an array of arrays of the elements of the covariance matrix, devs is an array of the reciprocals of the standard deviations of each element of the CDF's argument x. The threshold, steps, qrnd and prnd arguments are the convergence threshold, the maximum number of steps, the quasi random and pseudo random number generators required for our quasi random integrator.
We first count the number of elements for which the CDF is less than one. If there aren't any then we return one and if there is just one then we return its value for the univariate CDF.
Otherwise, we create a new array of arrays s that we populate with the elements of sigma whose row and column indices are those of the finite elements of x, which we keep track of in the array f. Similarly, the array y is filled with those finite elements and we then simply return the value of the CDF for them.

The corrNdCDF helper function given in listing 10 glues these functions together to implement the multivariate normal CDF.

Listing 10: The corrNdCDF Helper Function
function corrNdCDF(mu, sigma, root, devs, x, zero, one,
cdf, inv, y, threshold, steps, qrnd, prnd) {
var n = mu.length;
var i, c, inf, e0, f;

if(ak.type(x)!==ak.VECTOR_T) {
throw new Error('invalid argument in ak.multiNormalCDF');
}

if(x.dims()!==n) {
throw new Error('dimension mismatch in ak.multiNormalCDF');
}

x = x.toArray();
for(i=0;i<n;++i) x[i] -= mu[i];

if(x.some(isNaN)) return ak.NaN;

inf = 0;
for(i=0;i<n && inf!==-1;++i) {
c = cdf(x[i]*devs[i]);
if(c===0 || c===1) inf = 2*c-1;
}

if(inf===-1) {
return 0;
}

if(inf===+1) {
return corrNdMarginalCDF(sigma, devs, x, cdf, threshold, steps, qrnd, prnd);
}

e0 = cdf(x[0]*root[0][0]);
f = corrNdIntegrand(root, x, cdf, inv, e0, y);
return ak.quasiRandomIntegral(f, threshold, steps, qrnd, prnd)(zero, one);
}


This first checks that the CDF's argument x is a vector of the correct number of dimensions before subtracting the mean to centre its elements on zero, returning NaN if any of them are NaN. It then loops over the elements testing whether the standard univariate normal cdf returns zero or one, indicating that they are negative or positive infinity respectively. If any are the former then it returns zero and if any are the latter then it forwards to the corrNdMarginalCDF function to remove them from the integral. Finally, it uses ak.quasiRandomIntegral to integrate the function returned by corrNdIntegrand from zero to one.

Finally, the makeCorrNdCDF helper function given in listing 11 creates a function that takes the CDF's argument x and forwards to the corrNdCDF function to compute the value of the CDF, having created the ak.vector objects zero and one for the bounds of integration, the univariate cdf and its inverse inv, the reciprocals of the standard deviations devs and having reciprocated the diagonal elements of root.

Listing 11: The makeCorrNdCDF Helper Function
function makeCorrNdCDF(mu, sigma, root, threshold, steps, qrnd, prnd) {
var n = mu.length;
var zero = ak.vector(n-1, 0);
var one = ak.vector(n-1, 1);
var cdf = ak.normalCDF();
var inv = ak.normalInvCDF();
var devs = new Array(n);
var y = new Array(n-1);
var i;

for(i=0;i<n;++i) {
root[i][i] = 1/root[i][i];
devs[i] = 1/Math.sqrt(sigma[i][i]);
}

return function(x) {
return corrNdCDF(mu, sigma, root, devs, x, zero, one,
cdf, inv, y, threshold, steps, qrnd, prnd);
};
}


### Putting It All Together

We're now ready to implement the general purpose multivariate normal CDF with the ak.multiNormalCDF function given in listing 12.

Listing 12: ak.multiNormalCDF
ak.multiNormalCDF = function() {
var state = {mu: [0], sigma: [1], root: [1]};
var arg0  = arguments[0];
var n, i, cdf, f;

constructors[ak.type(arg0)](state, arg0, arguments);

n = state.sigma.length;

if(n===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) {
cdf = ak.normalCDF();
for(i=0;i<n;++i) state.root[i] = 1/state.root[i];

f = function(x){return uncorrCDF(state.mu, state.root, cdf, x);};
f.sigma = function(){return ak.matrix('diagonal', state.sigma);};
}
else {
if(n===2) {
f = makeCorr2dCDF(state.mu, state.sigma, state.threshold, state.steps);
}
else if(n===3) {
f = makeCorr3dCDF(state.mu, state.sigma, state.threshold, state.steps);
}
else {
f = makeCorrNdCDF(state.mu, state.sigma, state.root,
state.threshold, state.steps, state.qrnd, state.prnd);
}

f.sigma = function(){return ak.matrix(state.sigma);};
}
f.mu = function(){return ak.vector(state.mu);};

return Object.freeze(f);
};


This first checks for an uncorrelated CDF by seeing if the elements of sigma are numbers and, if not, forwards to the two, three or general multivariate CDF helper functions depending upon the length of state.sigma array of arrays of the elements of the covariance matrix.

Now we need two more constructors to allow the user to specify the parameters of the quasi random integrator, as given in listing 13.

Listing 13: The Multivariate Integration Constructors
constructors[ak.VECTOR_T][ak.MATRIX_T][ak.NUMBER_T] =
function(state, threshold, args) {
state.threshold = threshold;
state.steps = args[3];
state.qrnd = args[4];
state.prnd = args[5];
};

constructors[ak.MATRIX_T][ak.NUMBER_T] =
function(state, threshold, args) {
state.threshold = threshold;
state.steps = args[2];
state.qrnd = args[3];
state.prnd = args[4];
};


Note that, if these are undefined, ak.quasiRandomIntegral will provide defaults for them.

To demonstrate that ak.multiNormalCDF works correctly we shall need to exercise these four different algorithms, which we can do by comparing its values for an uncorrelated CDF and for correlated two, thee and four dimensional CDFs to sample based approximations, just as we did in programs 1 and 2.
This is done in programs 3, 4, 5 and 6 respectively with each using one hundred thousand samples to approximate the CDF rather than the ten thousand that we used in programs 1 and 2.

Program 3: An Uncorrelated CDF

Program 4: A Two Dimensional Correlated CDF

Program 5: A Three Dimensional Correlated CDF

Program 6: A Four Dimensional Correlated CDF

These programs confirm that the CDFs created by ak.multiNormalCDF are indeed consistent with the random samples, giving us some confidence that we have correctly implemented the four cases.
Derivation 4: The Symmetry Of The CDF
Firstly, for an $$n$$ dimensional normally distributed random variable we have
\begin{align*} \Phi_{-\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{-x}) &= \int_{-\mathbf{\infty}}^{-\mathbf{x}\,-\,-\mathbf{\mu}} \phi_{\mathbf{0},\mathbf{\Sigma}}(\mathbf{z}) \; \mathrm{d}\mathbf{z}\\ &= (-1)^n \int_{-\mathbf{x}\,-\,-\mathbf{\mu}}^{-\mathbf{\infty}} \phi_{\mathbf{0},\mathbf{\Sigma}}(\mathbf{z}) \; \mathrm{d}\mathbf{z} \end{align*}
Next, if we define $$\mathbf{z}^\prime = -\mathbf{z}$$ we have
$\mathrm{d}\mathbf{z} = \left|\frac{\partial \mathbf{z}}{\partial \mathbf{z}^\prime}\right| \mathrm{d}\mathbf{z}^\prime = \left|-\frac{\partial \mathbf{z}}{\partial \mathbf{z}}\right| \mathrm{d}\mathbf{z}^\prime = \left|-\mathbf{I}\right| \mathrm{d}\mathbf{z}^\prime = (-1)^n \mathrm{d}\mathbf{z}^\prime$
where $$\mathbf{I}$$ is the identity matrix, and so substituting $$\mathbf{z}^\prime$$ for $$\mathbf{z}$$ in the integral yields
\begin{align*} \Phi_{-\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{-x}) &= (-1)^n \int_{-(-\mathbf{x}\,-\,-\mathbf{\mu})}^{-(-\mathbf{\infty})} \phi_{\mathbf{0},\mathbf{\Sigma}}(\mathbf{z}^\prime) \; (-1)^n \mathrm{d}\mathbf{z}^\prime\\ &= \int_{\mathbf{x}-\mathbf{\mu}}^{\mathbf{\infty}} \phi_{\mathbf{0},\mathbf{\Sigma}}(\mathbf{z}^\prime) \; \mathrm{d}\mathbf{z}^\prime = \bar{\Phi}_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) \end{align*}

### The Complementary CDF

Finally, we must implement the complementary CDF $$\bar{\Phi}$$ which measures the probability that a multivariate normal random variable will yield a vector $$\mathbf{X}$$ whose elements are all strictly greater than those of a given vector $$\mathbf{x}$$
\begin{align*} \mathbf{X} &\sim \mathbf{N}(\mathbf{\mu},\mathbf{\Sigma})\\ \bar{\Phi}_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) &= \Pr\left(X_0 > x_0 \wedge X_1 > x_1 \wedge \dots \right)\\ &=\int_{\mathbf{x}}^{\mathbf{\infty}} \phi_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{z}) \; \mathrm{d}\mathbf{z} \end{align*}
We can do this by exploiting the symmetry of the multivariate normal distribution about its mean; specifically, that
$\bar{\Phi}_{\mathbf{\mu},\mathbf{\Sigma}}(\mathbf{x}) = \Phi_{\mathbf{-\mu},\mathbf{\Sigma}}(\mathbf{-x})$
as shown by derivation 4 and used in the implementation of ak.multiNormalCompCDF given in listing 14.

Listing 14: ak.multiNormalCompCDF
ak.multiNormalCompCDF = function() {
var state = {mu: [0], sigma: [1], root: [1]};
var arg0  = arguments[0];
var nmu = [];
var n, i, cdf, f, g;

constructors[ak.type(arg0)](state, arg0, arguments);

n = state.sigma.length;
nmu.length = n;
for(i=0;i<n;++i) nmu[i] = -state.mu[i];

if(n===0 || ak.nativeType(state.sigma[0])===ak.NUMBER_T) {
cdf = ak.normalCDF();
for(i=0;i<n;++i) state.root[i] = 1/state.root[i];

f = function(x){return uncorrCDF(nmu, state.root, cdf, ak.neg(x));};
f.sigma = function(){return ak.matrix('diagonal', state.sigma);};
}
else {
if(n===2) {
g = makeCorr2dCDF(nmu, state.sigma, state.threshold, state.steps);
}
else if(n===3) {
g = makeCorr3dCDF(nmu, state.sigma, state.threshold, state.steps);
}
else {
g = makeCorrNdCDF(nmu, state.sigma, state.root,
state.threshold, state.steps, state.qrnd, state.prnd);
}

f = function(x){return g(ak.neg(x));};
f.sigma = function(){return ak.matrix(state.sigma);};
}
f.mu = function(){return ak.vector(state.mu);};

return Object.freeze(f);
};


Program 7 demonstrates that this too is consistent with the proportion of a random sample of multivariate normally distributed vectors having all of their elements greater than those of its argument.

Program 7: A Correlated Complementary CDF

And with that we're done, except for me to remind you that all of these functions can be found in MultiNormalDistribution.js.

### References

[1] Every Which Way Is Normal, www.thusspakeak.com, 2017.

[2] Multiple Multiply Normal Functions, www.thusspakeak.com, 2017.

[3] The Nesting Instinct, www.thusspakeak.com, 2016.

[4] Drezner, Z. and Wesolowsky, G., On the computation of the bivariate normal integral, Journal of Statistical Computation and Simulation, Vol. 35, No. 1-2, 1990.

[5] Genz, A., Numerical computation of rectangular bivariate and trivariate normal and t probabilities, Statistics and Computing, Vol. 14, No. 3, 2004.

[6] Gupta, S., Probability Integrals of Multivariate Normal and Multivariate t, The Annals of Mathematical Statistics, Vol. 34, No. 3, 1963.

[7] Integration, Quasi Mode, www.thusspakeak.com, 2017.

[8] Genz, A. and Bretz, F., Computation of Multivariate Normal and t Probabilities, Springer, 2009.