mnorm: An R Package for Calculation and Differentiation of Conditional Multivariate Normal Densities and Probabilities

mnorm: An R Package for Calculation and Differentiation of Conditional Multivariate Normal Densities and Probabilities

\noindentKeywords: multivariate normal distribution, analytical gradients, econometrics

JEL Codes: C63, C87

Introduction

The multivariate normal distribution (MVN) has various applications including maximum-likelihood estimation and Bayesian inference. There are many R packages dealing with MVN. The most popular one is mvtnorm mvtnorm. It provides tools for calculation of multivariate normal probabilities and densities. For many practical applications, these tools are sufficient. However, sometimes it is necessary to get a boost in terms of speed, accuracy, or functionality. To increase calculation speed in the low-dimensional case, the pbivnorm package may be used, which provides fast calculation of bivariate normal probabilities. The mvnfast package allows for fast parallel (multi-core) calculation of MVN densities. To calculate conditional or truncated MVN probabilities, one may use the condMVNorm and TruncatedNormal packages, respectively. The \textbf{tlrmvnmvt} package tlrmvnmvt implements methods suitable for accurate calculation of MVN probabilities in the high-dimensional case. To analyse MVN data, one may use the \textbf{norm} and \textbf{MVN} MVN packages. The \textbf{copula} package copula allows one to simulate from the Gaussian copula, which extends MVN to the case of non-normal margins.

The mnorm package aims to extend the MVN tools available in R mainly by providing functions for differentiation (calculation of gradients) of MVN probabilities and densities with respect to various parameters including covariances and integration limits. Particularly, it may be useful for maximization of likelihood functions for which gradient information may substantially speed up a numeric optimization routine. For example, the switchSelection package uses mnorm for estimation of multivariate and multinomial sample selection and endogenous switching econometric models.

Also, the package allows one to calculate and differentiate conditional MVN probabilities. Particularly, it is possible to calculate gradients with respect to conditioned (given) components of the distribution (elements of the MVN vector). Moreover, a user may provide an alternative (non-normal) specification of marginal distributions (i.e., deal with the Gaussian copula) and differentiate probabilities with respect to the parameters of these distributions. To facilitate speed, the package is written in Rcpp Rcpp and intensively utilizes RcppArmadillo RcppArmadillo. It also uses Open Multi-Processing (OMP) for parallel computing.

The article proceeds as follows. In Section (ref), we provide theoretical information on the MVN distribution. In Section (ref), we demonstrate the application of \detokenize{dmnorm} and \detokenize{pmnorm} functions allowing one to calculate and differentiate (conditional) MVN densities and probabilities, respectively. In Section (ref), we compare the efficiency (speed) of our analytical gradients to the lower bound time needed to compute numeric gradients. Moreover, we compare the calculation speed of densities and probabilities to that of other popular packages. In Section (ref), we provide the derivation of the differentiation formulas implemented in the package. In Section (ref), we derive some simple, useful bounds on the MVN log-probabilities and their derivatives, which may accompany the usage of the package. Particularly, we contribute to the literature on the bounds for inverse Mills ratios From2020 by providing a useful bound on the derivatives of log-probabilities.

Multivariate normal distribution

Consider an $m$-variate random row vector which has the multivariate normal distribution:

Parameters of this distribution are the row vector of expectations and the covariance matrix:

The density function of the multivariate normal vector $X$ is as follows:

where $x = (x_{1},...,x_{m})$ is a row vector and $\text{det}\left(\Sigma\right)$ is the determinant of the covariance matrix $\Sigma$.

Probabilities of the multivariate normal distribution are as follows:

where $x^{(l)}_{j}, x^{(u)}_{j}\in \mathbb{R}\cup\{-\infty,\infty\}$ are the elements of the row vectors $x^{(l)}$, $x^{(u)}$, representing lower and upper integration limits, respectively.

Suppose that some of the components of $X$ are conditioned. Denote by $I_{g}$ a vector of indices of the conditioned components. The following vector has a conditional multivariate normal distribution:

where the vector of conditional expectations and the conditional covariance matrix are as follows:

Denote by $\Sigma_{(A, B)}$ a submatrix of $\Sigma$, where $A,B\in\{d, g\}$. This matrix is constructed from $\Sigma$ by the intersection of rows and columns with indices $I_{A}$ and $I_{B}$, respectively. Particularly, $\Sigma_{d,g}^{T}=\Sigma_{g,d}$.

Random vectors $X_{d}$ and $X_{g}$ are subvectors of $X$, with elements with indices $I_{d}$ and $I_{g}$, respectively. Similarly, $\mu_{d}$ and $\mu_{g}$ are subvectors of $\mu$. The density function of the conditional multivariate normal distribution at point $x^{(d)}$ is denoted by $f\left(x^{(d)};\mu,\Sigma | x^{(g)},I_{g}\right)$, where the row vector $x^{(g)}$ contains elements at which the elements of $X_{g}$ are conditioned.

The numbers of conditioned and unconditioned elements of $X$ are as follows:

Consider $i\in I_{d}$ and $j\in I_{g}$. Also, denote by $I_{g,k}$ and $I_{d,k}$ the $k$-th elements of the vectors $I_{g}$ and $I_{d}$, respectively. Let

where $q_{d}(i)$ is the index of the $i$-th element of $X$ in $I_{d}$, and $q_{g}(j)$ is the index of the $j$-th element of $X$ in $I_{g}$; hence $x_{i}=x^{(d)}_{q_{d}(i)}$ and $x_{j}=x^{(g)}_{q_{g}(j)}$. Finally, denote by $q_{g}^{-1}(\cdot)$ and $q_{d}^{-1}(\cdot)$ the inverse functions of $q_{g}(\cdot)$ and $q_{d}(\cdot)$, respectively. For example, let $I_{d} = \{1,3,5\}$ and $I_{g} = \{2, 4\}$. Then $q_{d}(3)=2$ and $q_{g}(2) = 1$. Moreover, $q_{d}^{-1}(2)=3$ and $q_{g}^{-1}(1) = 2$.

The notation $M(a, b; m_{1}, m_{2})$ corresponds to an $\left(m_{1}\times m_{2}\right)$ matrix whose elements are $0$, except for the $(a, b)$-th element, which is equal to $1$. Similarly, $M(a, b; m_{1})$ is a square matrix $\left(m_{1}\times m_{1}\right)$ of zeros, except that the elements with indices $(a, b)$ and $(b, a)$ are equal to $1$.

Denote by $f_{X_{i}}\left(.\right)$ the density function of $X_{i}$. Similarly, denote by $f_{X_{i},X_{j}}\left(.\right)$ the joint density function of $X_{i}$ and $X_{j}$. The notation $x^{(l)}_{-(i, j)}$ corresponds to the vector $x^{(l)}$ without the $i$-th and $j$-th elements. The notations $x^{(u)}_{-(i, j)}$ and $X_{-(i, j)}$ are defined in the same way.

The cumulative distribution function of the multivariate normal distribution is as follows:

Examples

Calculation and differentiation of densities

The function \detokenize{dmnorm} calculates MVN densities $f(x;\mu,\Sigma)$. Its interface is as follows:

The main arguments of this function are a numeric vector (or matrix) \detokenize{x}, a numeric vector \detokenize{mean}, and a symmetric positive-definite numeric matrix \detokenize{sigma}, representing $x$, $\mu$, and $\Sigma$, respectively. The function returns a list which contains information on the density value (named \detokenize{den}). If \detokenize{grad_x = TRUE}, then the output list also contains the gradient with respect to $x$ named \detokenize{grad_x}. Similarly, if \texttt{\detokenize{grad_sigma = TRUE}}, then the output list additionally contains the gradient with respect to $\Sigma$ named \texttt{\detokenize{grad_sigma}}.

To demonstrate an application of this function, we start with a simple example and gradually move to more advanced cases. Consider a $5$-dimensional MVN vector $X$ with the following parameters:

To calculate $f(x;\mu,\Sigma)$ at the point $x = (-1, 0, 1, 2, 3)$, type:

Note that the function returns a list whose element named \detokenize{den} contains the density value:

The function is vectorized with respect to \detokenize{x}. If \detokenize{x} is a matrix, then the density values will be calculated for each row of \detokenize{x}. For example, to additionally calculate the density at the point $y = (2.5, 2, 1.5, 1, 0.5)$, add the following code:

Therefore, \detokenize{d2$den[1]}} and \texttt{\detokenize{d2$den[2]} are the calculated values of $f(x;\mu,\Sigma)$ and $f(y;\mu,\Sigma)$, respectively.

To calculate the conditional densities, indices of conditioned elements $I_{g}$ should be provided through the \detokenize{given_ind} argument. Then, elements of \detokenize{x} with corresponding indices will be treated as given values $x^{(g)}$ and other elements will represent $x^{(d)}$. For example, suppose that $X_{1}$ and $X_{3}$ are conditioned at values $-1$ and $2$, respectively, so $x^{(g)}=(-1, 2)$. Then, to calculate the density function of:

at the point $x^{(d)} = (0, 0.5, 1)$, namely:

add the following code:

Note that vectorization with respect to \detokenize{x} is still applicable for the conditional densities. At the same time, the indices of the conditioned elements \detokenize{given_ind} remain identical for all the observations. For example, to additionally calculate the density of:

at the point $y^{(d)} = (2, 1, 0.5)$, add the code:

All aforementioned densities may be differentiated with respect to various parameters. For example, to calculate the gradient of $f(x;\mu,\Sigma)$ with respect to $x$, set \detokenize{grad_x = TRUE} and get the element named \detokenize{grad_x} from the output list:

Therefore, \detokenize{d5$grad_x}} represents $\nabla f(x;\mu,\Sigma)$. For example, \texttt{\detokenize{d5$grad_x[3]} is $\frac{\partial f(x;\mu,\Sigma)}{\partial x_{3}}$.

To simultaneously calculate gradients at the points $x$ and $y$ (Jacobian), one may use vectorization:

Note that \detokenize{grad_x[i, ]} is the gradient at the point \detokenize{x[i, ]}. Particularly, since the first row of \detokenize{x_mat} is $x$ and the second row is $y$, \detokenize{d6$grad_x[1, 3]}} is $\frac{\partial f(x;\mu,\Sigma)}{\partial x_{3}}$ and \texttt{\detokenize{d6$grad_x[2, 5]} is $\frac{\partial f(y;\mu,\Sigma)}{\partial y_{5}}$.

To differentiate the density function with respect to the elements of the covariance matrix, set \detokenize{grad_sigma = TRUE} and get the element named \detokenize{grad_sigma} from the output list:

Note that the output value \detokenize{grad_sigma} is a $3$-dimensional array such that \detokenize{grad_sigma[, , i]} is the gradient with respect to $\Sigma$ at the point corresponding to \detokenize{x[i, ]}. For example, \detokenize{d7$grad_sigma[2, 1, 1]}} is $\frac{\partial f(x;\mu,\Sigma)}{\partial \Sigma_{2,1}}$ and \texttt{\detokenize{d7$grad_sigma[3, 5, 2]} is $\frac{\partial f(y;\mu,\Sigma)}{\partial \Sigma_{3,5}}$.

Finally, let's differentiate both aforementioned conditional densities:

Note that the output value \detokenize{grad_x} contains derivatives with respect to $x^{(d)}$ and $x^{(g)}$. To distinguish between these two types of derivatives, it is convenient to think of \detokenize{grad_x} as a gradient with respect to the input argument \detokenize{x}. For example, \detokenize{d8$grad_x[1, 3]}} is $\frac{\partial f\left(x^{(d)};\mu,\Sigma | x^{(g)},I_{g}\right)}{\partial x^{(g)}_{2}}$ since \texttt{\detokenize{z_mat[1, 3]}} represents the value of $x^{(g)}_{2}$. Similarly, \texttt{\detokenize{d8$grad_x[2, 5]} is $\frac{\partial f\left(y^{(d)};\mu,\Sigma | y^{(g)},I_{g}\right)}{\partial y^{(d)}_{3}}$ because \detokenize{z_mat[2, 5]} is $y^{(d)}_{3}$. Interpretation of derivatives with respect to $\Sigma$ is identical to the aforementioned case with unconditional densities.

To calculate and differentiate the log-densities, set \detokenize{log = TRUE}. The number of cores to be used for a parallel computation may be set via the \detokenize{n_cores} argument. If \detokenize{validation = FALSE}, then validation of the input arguments will be disabled, which may result in a slight boost in terms of speed. These arguments work similarly for \detokenize{pmnorm}.

Calculation and differentiation of probabilities

The function \detokenize{pmnorm} calculates MVN probabilities $\text{P}\left(x^{(l)}\leq X\leq x^{(u)}\right)$. Its interface is as follows:

To calculate the bivariate probabilities (by default), the function applies the method of Drezner and Wesolowsky, described in Genz2004. In contrast to the classical implementation of this method, \detokenize{pmnorm} applies Gauss-Legendre quadrature with $30$ sample points to approximate the integral (1) of Genz2004. The classical implementation of this method uses up to 20 points but requires some additional transformations of (1). During preliminary testing, it has been found that the approach with $30$ points provides similar accuracy, being slightly faster because of better vectorization and parallel computing capabilities.

To calculate the trivariate ($3$-dimensional) probabilities (by default), the function uses the Drezner method, following formula (14) of Genz2004. Similarly to the bivariate case, 30 points are used in the Gauss-Legendre quadrature.

For the $m$-variate probabilities, where $m > 3$, the function applies the Geweke-Hajivassiliou-Keane (GHK) algorithm described in section 4.2 of Genz2009. The implementation of GHK is based on the deterministic Richtmyer sequence with \detokenize{n_sim} draws and uses variable reordering suggested in Section 4.1.3 of Genz2009. The ordering algorithm may be determined via the \detokenize{ordering} argument. Available options are \detokenize{NO}, \detokenize{mean} (default), and \detokenize{variance}. To switch from the Richtmyer sequence to another alternative, the \detokenize{random_sequence} parameter of the \texttt{\detokenize{control}} input argument may be used. A user may supply an $m$-column matrix of uniform (between 0 and 1) numbers generated via some algorithm. For example, the scrambled Halton sequence may be generated via the \texttt{\detokenize{halton}} function of the \texttt{\detokenize{mnorm}} package. The bivariate and trivariate probabilities may also be calculated via the GHK algorithm (not recommended in terms of speed and accuracy) by setting \texttt{\detokenize{method = "GHK"}}.

If \detokenize{method = "Gassmann"}, then the function applies the method of Gassmann2003 for the calculation of the $m>3$-dimensional normal probabilities. It uses the matrix representation (5) of Gassmann2003 and $30$ points of the Gauss-Legendre quadrature.

The main arguments are \detokenize{lower}, \detokenize{upper}, \detokenize{mean}, and \detokenize{sigma} representing $x^{(l)}$, $x^{(u)}$, $\mu$, and $\Sigma$, respectively. Note that \detokenize{pmnorm} is optimized to perform much faster (especially in the $2$- and $3$-dimensional cases) when some integration limits are infinite. The function returns a list which contains information on the probability value (named \detokenize{prob}). The list may also contain the derivatives named after the gradient-related arguments \texttt{\detokenize{grad_lower}}, \texttt{\detokenize{grad_upper}}, \texttt{\detokenize{grad_sigma}}, \texttt{\detokenize{grad_given}}, and \texttt{\detokenize{grad_marginal}}. For example, if \texttt{\detokenize{grad_upper = TRUE}}, then the output list will contain the element named \texttt{\detokenize{grad_upper}} which stores the information on the derivatives with respect to the upper integration limits. Let's clarify the application of the \texttt{\detokenize{pmnorm}} function and the structure of its output by investigating a series of pedagogical examples.

Consider a $5$-dimensional MVN vector from the previous section. Let's calculate the following probability:

To calculate this probability, run the following code:

The function is vectorized with respect to the \detokenize{lower} and \detokenize{upper} arguments. If these arguments are matrices, then the limits should be provided in the rows. Let's additionally calculate the probability with the following limits:

To simultaneously calculate both aforementioned probabilities, add the code:

Therefore, \detokenize{p2$prob[1]}} and \texttt{\detokenize{p2$prob[2]} are the calculated values of the MVN probabilities under the limits $x^{(l)}$, $x^{(u)}$ and $y^{(l)}$, $y^{(u)}$, respectively.

Let's calculate the conditional probabilities:

and:

where:

The conditional values $x^{(g)}$ and $y^{(g)}$ should be provided as the rows of the vectorized argument \detokenize{given_x}. The indices of the conditioned components $I_{g}$ should be supplied to the \detokenize{given_ind} argument. Similarly to the previous example, $x^{(l)}, x^{(u)}$ and $y^{(l)}, y^{(u)}$ should be the rows of the \detokenize{lower} and \detokenize{upper} arguments, respectively. Therefore, to calculate these conditional probabilities, add the following code:

The arguments \detokenize{grad_lower}, \detokenize{grad_upper}, \detokenize{grad_sigma}, and \detokenize{grad_given} allow one to differentiate the aforementioned conditional probabilities with respect to $x^{(l)}$ and $y^{(l)}$, $x^{(u)}$ and $y^{(u)}$, $\Sigma$, and $x^{(g)}$ and $y^{(g)}$, respectively:

Note that \detokenize{grad_lower[i, ]}, \detokenize{grad_upper[i, ]}, and \detokenize{grad_given[i, ]} are the gradients of the $i$-th probability with respect to the lower integration limits, the upper integration limits, and the conditioned values, respectively. For example, \detokenize{grad_lower[1, 2]} is $\frac{\partial P_{1}}{\partial x^{(l)}_{2}}$, \detokenize{grad_upper[2, 3]} is $\frac{\partial P_{2}}{\partial y^{(u)}_{3}}$, and \detokenize{grad_given[2, 1]} is $\frac{\partial P_{2}}{\partial y^{(g)}_{1}}$. The value of \texttt{\detokenize{grad_sigma[j, k, i]}} represents the derivative of the $i$-th probability with respect to $\Sigma_{j, k}$. Particularly, \texttt{\detokenize{grad_sigma[4, 5, 1]}} is $\frac{\partial P_{1}}{\partial \Sigma_{4, 5}}$.

Finally, let's substitute the marginal distributions using the \detokenize{marginal} argument. It should be a list such that \detokenize{names(marginal)[i]} is the name of the $i$-th distribution and \detokenize{marginal[[i]]} represents its parameters. Possible names (distributions) are \detokenize{"PGN"} (or \detokenize{"hpa"}), \detokenize{"student"} (or \texttt{\detokenize{"t"}}), \texttt{\detokenize{"logistic"}}, and \texttt{\detokenize{"normal"}}. Normal and logistic distributions should take \texttt{\detokenize{NULL}} or an empty vector \texttt{\detokenize{numeric()}} as the parameter. For the Student distribution, it should be a single value representing the degrees of freedom. To incorporate the PGN distribution proposed by PGN, we use the \texttt{\detokenize{phpa0}} function from the \textbf{hpa} package, allowing one to calculate the probabilities of the standardized PGN distribution as described in Trifonov. Therefore, the values from the corresponding elements of \texttt{\detokenize{marginal}} are transferred to the argument \texttt{\detokenize{pc}} of the \texttt{\detokenize{phpa0}} function.

Suppose that $X_{1}$ and $X_{2}$ have the PGN distribution (with some parameters), $X_{3}$ has the Student distribution with $5$ degrees of freedom, $X_{4}$ has the standard logistic distribution, and $X_{5}$ has the normal distribution. The distribution of each $X_{i}$ is automatically transformed to have mean $\mu_{i}$ and variance $\Sigma_{i,i}$. To calculate the probabilities given these marginal distributions (connected by the Gaussian copula), run the following code:

By setting \detokenize{grad_marginal = TRUE}, we additionally calculate the derivatives of the probabilities with respect to the parameters of the marginal distributions. Since normal and logistic distributions have no parameters, \detokenize{p5$grad_marginal[[4]]}} and \texttt{\detokenize{p5$grad_marginal[[5]]} are empty. The value of \detokenize{p5$grad_marginal[[3]][i]}} represents the derivative of the $i$-th probability with respect to the degrees of freedom of the Student distribution. Finally, \texttt{\detokenize{p5$grad_marginal[[1]][i, ]} and \detokenize{p5$grad_marginal[[2]][i, ]}} contain the gradient of the $i$-th probability with respect to the parameters of the PGN distributions associated with $X_{1}$ and $X_{2}$, respectively

Efficiency analysis

Derivatives

To motivate the application of the analytical derivatives, we compare their computation speed (in seconds) to that of numeric derivatives. The numeric gradient requires at least one more call to the function than the number of its elements. For example, to get the numeric gradient of $f(x)$ with respect to $x$ in the $m=5$-dimensional case, the numeric gradient requires at least $m+1=6$ calls to $f(x)=f(x_{1},...,x_{5})$. We will use this fact to calculate a lower bound of the time needed to compute the numeric gradient. Note that accurate numeric differentiation techniques may require $2$ or more calls to the function per parameter. Also, this lower bound ignores the time needed to combine the calls to the function into a gradient. Nevertheless, the following results suggest that in most cases the calculation time of the analytical derivatives of the MVN densities and probabilities notably outperforms this lower bound.

Consider the same MVN distribution setup as in the previous sections. First, let's compare the calculation time of the densities, the conditional densities, and their derivatives on $10^6$ simulations.

The computation time for $f(x)$ and $f(x^{(d)}; x^{(g)})$, denoted by \detokenize{den_uncond} and \detokenize{den_cond}, respectively (similarly for the following derivatives of these densities), is as follows:

The computation time of the gradients with respect to $x$ (unconditional) and $x^{(d)}$, $x^{(g)}$ (conditional), where prefix \detokenize{num} stands for the computed value of the lower bound of time needed to calculate the numeric derivatives (similarly for further examples), is as follows:

The computation time of the gradients with respect to all the parameters, i.e., $x$, $x^{(d)}$, $x^{(g)}$, and $\Sigma$, is as follows:

The results suggest that computation time for the unconditional and conditional densities is just slightly greater than that of their derivatives with respect to $x$, $x^{(d)}$, and $x^{(g)}$. As a result, analytical derivatives with respect to these parameters are much faster than numeric derivatives. In contrast, the time needed to calculate the analytical derivatives of the (conditional) densities with respect to all the parameters (including elements of $\Sigma$) is just mildly smaller than the minimum time needed to compute the numeric derivatives of the (conditional) densities. It may be associated with the fact that the analytical derivatives with respect to $\Sigma$ are rather complicated, especially in the conditional case.

Let's use $10^3$ simulations to calculate the efficiency of the analytical gradients for the unconditional probabilities.

The results indicate that the analytical gradients are much faster than the numeric gradients for $\Sigma$ and just slightly faster for $x^{(u)}$. It may be explained by the fact that the derivatives of the probabilities with respect to $x^{(u)}$ require the calculation of the $(m-1)$-dimensional integrals, while the derivatives with respect to $\Sigma$ require the calculation of the $(m-2)$-dimensional integrals. Notice that in the $5$-dimensional case, the $4$-dimensional integrals are still calculated via the GHK algorithm, which is much slower than the Drezner method used for the $3$-dimensional probabilities. This explains why the extra time needed to calculate the derivatives with respect to $\Sigma$ is so small.

Let's reproduce the experiment for the conditional probabilities using $10^5$ simulations.

For the conditional probabilities, the difference in the calculation speed between analytical and numeric derivatives is more substantial. This is because computation of the $5$-dimensional conditional MVN probabilities requires the $m_{d}=3$-dimensional integration (since $2$ out of $5$ components are conditioned). Therefore, the derivatives are calculated extremely fast via the $1$-dimensional and $2$-dimensional integration methods. Finally, let's consider the $10$-dimensional case when all the derivatives will be calculated using the probabilities calculated via the GHK method. Let's use $10^3$ simulations, a zero mean vector and a constant correlation matrix with the correlation coefficients equal to $0.5$.

In line with the aforementioned expectations, the results demonstrate a substantial decline in the relative performance of the analytical derivatives. Nevertheless, their calculation time is still slightly less than the lower bound for the numeric derivatives.

Finally, notice that the derivatives with respect to $x^{(g)}$ and $\Sigma_{i,j}$ (such that $i\in I_{g}$ or $j\in I_{g}$) are just functions of the derivatives with respect to $x^{(u)}$ and $\Sigma_{k, l}$ (such that $k,l\in I_{d}$). Hence, extra time needed to calculate these derivatives is very low. Consequently, the aforementioned decline should be much lower for the conditional probabilities. Indeed, suppose that $I_{g}=\{1,2,3,4\}$. Then both the absolute and the relative (to the numeric gradient) time to calculate the analytical derivatives with respect to $x^{(u)}$, $x^{(g)}$ and $\Sigma$ decrease substantially.

Densities

We compare the efficiency of MVN density computation by mnorm with that of the mvnfast and mvtnorm packages, which calculate these densities via the functions \detokenize{dmvn} and \detokenize{dmvnorm}, respectively. To do so, we simulate $10$ million observations from the $5$-dimensional MVN distribution. Then we compute the MVN density value at each of the simulated points and calculate computation time (in seconds) for each function. Since both \detokenize{dmnorm} and \texttt{\detokenize{dmvn}} support parallel computing, we also compare these functions using $8$ cores.

The results suggest that both \detokenize{dmnorm} and \detokenize{dmvn} notably outperform \detokenize{dmvnorm}. Also, under the same number of cores, \detokenize{dmnorm} is faster than \detokenize{dmvn}.

Let's replicate the example for the $100$-dimensional case with a zero mean vector and a constant correlation matrix whose correlation coefficients equal $0.5$. We consider $10^5$ simulations.

The results suggest that, in the high-dimensional case, all the non-parallel functions provide similar computation speed. However, the parallel implementations of \detokenize{dmnorm} and \detokenize{dmvn} substantially outperform \detokenize{dmvnorm}.

Bivariate and trivariate probabilities

First, we compare the speed of the \detokenize{pmnorm} and \detokenize{pbivnorm} functions using the standard bivariate normal distribution with correlation coefficients equal to $0.2$, $0.5$, $0.8$, and $0.95$.

The results suggest that, independent of the correlation coefficient value, \detokenize{pmnorm} is notably faster than \detokenize{pbivnorm}. To check that the accuracy of both functions is similar, we also calculate the maximum absolute and relative differences between the probabilities computed via \detokenize{pmnorm} and \detokenize{pbivnorm}.

The accuracy of both functions is very similar since the maximum absolute and relative differences in the calculated probabilities are extremely small.

To our knowledge, the fastest way to compute the trivariate probabilities in R is by using the TVPACK algorithm from the \detokenize{mvtnorm} package. Therefore, we compare our algorithm to this function on $20000$ simulations from the standard trivariate normal distribution with the correlation coefficients $0.8$, $0.7$, and $0.6$. Since \detokenize{mvtnorm} is not vectorized, we will use a loop to calculate the corresponding probabilities.

The results suggest that \detokenize{pmnorm} dramatically outperforms \detokenize{pmvnorm}. We have additionally checked that the difference in performance remains very similar even for large values of the \detokenize{abseps} input argument of \detokenize{TVPACK} (omitted for brevity). To ensure that the results are not driven by substantial accuracy losses, we have also checked that the absolute and relative differences in the calculated probabilities are small.

High-dimensional probabilities

To test the efficiency of \detokenize{pmnorm} in the high-dimensional case, we have replicated Tables 1-3 of tlrmvnmvt (using their code from the supplementary materials) for $16$, $64$, $128$, and $512$ dimensions. Our modified version of this code may be provided upon request. All tables have a constant correlation matrix. For Table 1 and Table 2, the correlation coefficients equal $0.5$. For Table 3 they equal $0.8$. All lower integration limits are $-\infty$. For Table 2 and Table 3, the upper integration limits are $-1$. For Table 1 they equal $0$. Following tlrmvnmvt, we use the default settings for all the functions. For \detokenize{pmnorm}, we set \detokenize{n_sim = 20000} since with this value its calculation time is similar to that of other functions. Columns of these tables represent the number of dimensions. Rows correspond to the functions or packages used to calculate the probabilities. The row labeled TruncatedNormal stands for the \detokenize{pmvnorm} function from the TruncatedNormal package. The rows labeled tlrmvnmvt (GenzBretz) and tlrmvnmvt (TLRQMC) represent the function \detokenize{pmvn} from \textbf{tlrmvnmvt} with the \texttt{\detokenize{algorithm}} argument set to \texttt{\detokenize{GenzBretz}} and \texttt{\detokenize{TLRQMC}}, respectively.

The results are similar (identical in terms of accuracy but slightly differ in terms of speed since we ran the code using another computer) to those reported in tlrmvnmvt and suggest that the speed and accuracy of \detokenize{pmnorm} (from mnorm) are comparable to those of the other functions. Therefore, it is reasonable to consider \detokenize{pmnorm} for calculation of the high-dimensional probabilities.

Formulas for the derivatives

In this section, we derive the formulas used in the package to calculate the derivatives of the MVN probabilities and densities. To simplify the derivations, suppose that $\mu$ is a vector of zeros. To apply the formulas from this section to the case when $\mu$ is not a vector of zeros, simply replace $x$ with $x-\mu$. Moreover, for brevity, $f(x;\mu,\Sigma)$, $F(x;\mu,\Sigma)$, and $\overline{F}(x^{(l)}, x^{(u)};\mu,\Sigma)$ are denoted by $f(x)$, $F(x)$, and $\overline{F}(x^{(l)}, x^{(u)})$, respectively.

Useful formulas from the linear algebra

Let $c$, $x$, and $A$ be a constant, a row vector, and a matrix, respectively. A matrix $U=U(c)$ is a function of the constant $c$. We will frequently use the following well-known formulas:

Moreover, we use the following lemma.

\noindentProof of Lemma (ref).

Density function

In this section, we provide the derivatives of the log-density function of $X$ with respect to the argument $x$ and the elements of the covariance matrix $\Sigma$. The notation $\Sigma_{ij}^{-1}$ corresponds to the element in the $i$-th row and the $j$-th column of the inverse matrix $\Sigma^{-1}$.

\noindentProof of Lemma (ref).

Conditional expectation

In this section, we provide the derivatives of the conditional expectation $\mu_{c}$ with respect to the fixed (given) values $x^{(g)}$ and the covariance matrix $\Sigma$.

\noindentProof of Lemma (ref).

Conditional covariance

In this section, we provide the derivatives of the elements of the conditional covariance matrix $\Sigma_{c}$ with respect to the elements of the unconditional covariance matrix $\Sigma$.

\noindentProof of Lemma (ref).

Conditional density function

In this section, we provide the derivatives of the conditional log-density of $X$ with respect to the argument $x^{(d)}$, the given values $x^{(g)}$, and the elements of the covariance matrix $\Sigma$. For simplicity, suppose that $\mu_{d}$ and $\mu_{g}$ are vectors of zeros. To relax this assumption, simply replace $x_{d}$ and $x_{g}$ with $x_{d}-\mu_{d}$ and $x_{g}-\mu_{g}$, respectively. Finally, for brevity, denote the conditional density function $f\left(x^{(d)};\mu,\Sigma | x^{(g)},I_{g}\right)$ as $f(x^{(d)}\mid x^{(g)})$.

\noindentProof of Lemma (ref).

Cumulative distribution function

In this section, we provide the derivatives of the cumulative distribution function of $X$ with respect to the argument $x$ and the elements of the covariance matrix $\Sigma$.

\noindentProof of Lemma (ref).

Probabilities

In this section, we provide the derivatives of the probability of the event that $X$ takes values between $x^{(l)}$ and $x^{(u)}$, with respect to these arguments and the elements of the covariance matrix $\Sigma$.

\noindentProof of Lemma (ref).

Conditional probabilities

In this section, we provide the derivatives of the conditional probabilities with respect to the integration limits $x^{(l)}$, $x^{(u)}$, and the elements of the covariance matrix $\Sigma$.

Specifically, we redefine $x^{(l)}$ and $x^{(u)}$ as the integration limits associated only with the components indexed by $I_{d}$. Moreover, for brevity, we introduce the following notations:

Despite $\tilde{\text{P}}=\text{P}$, sometimes it is more convenient to differentiate $\tilde{\text{P}}$ rather than $\text{P}$.

\noindentProof of Lemma (ref).

Bounds

In this section, we consider some useful bounds on the probabilities of the multivariate normal distribution and their derivatives. In particular, these inequalities are useful for proving consistency and asymptotic normality of maximum-likelihood estimators in econometric models such as the multivariate probit model. Therefore, these bounds will be useful for the users of the package.

Probabilities

In this section, we provide a lower bound on the logarithm of the cumulative distribution function of the multivariate normal distribution.

\noindentProof of Lemma (ref).

Notice that the lower bound of Lemma (ref) is less strict than the one obtained by Hashorva2003. However, the derived inequality may be useful in some practical scenarios due to its simplicity.

Derivatives of log-probabilities

In this section, we provide an upper bound on the derivatives of the log-probabilities of the multivariate normal distribution.

\noindentProof of Lemma (ref).

\noindentProof of Lemma (ref).

Conclusion

We have proposed the mnorm package which allows one to calculate and differentiate the conditional MVN densities and probabilities. The ability to calculate the gradients with respect to various parameters of MVN densities and probabilities is the main contribution of the package. Simulated data analysis suggests that analytical gradients notably outperform numeric gradients in terms of speed. The difference is especially great for probabilities when $m_{d}\leq 5$ because the derivatives are calculated via fast bivariate and trivariate computation methods. Moreover, the analytical derivatives with respect to the parameters associated with the given (conditioned) parts of the conditional distribution are very fast, since they are just functions of the derivatives with respect to other parameters. Therefore, we hope that our package will substantially ease the implementation of maximum-likelihood estimators involving calculation of MVN densities and probabilities. The package may be especially useful for econometric models involving the conditional MVN probabilities, for example, multivariate and multinomial sample selection and endogenous switching models {abort, {Dolgikh2023A}, {Dolgikh2023B}, {Dolgikh2024}, {Dolgikh2026gender}, {Dolgikh2026}}.

Also, according to the results of simulated data analysis, \detokenize{pmnorm} notably outperforms other popular alternatives in terms of calculation speed of the bivariate and trivariate MVN probabilities. Moreover, \detokenize{dmnorm} is slightly faster than the alternative functions for calculation of the MVN densities.

Finally, we have provided convenient bounds on the derivatives of the log-probabilities of the MVN distribution.