# Probabilities of (conditional) multivariate normal distribution This function calculates and differentiates probabilities of (conditional) multivariate normal distribution. ```r pmnorm( lower, upper, given_x = numeric(), mean = numeric(), sigma = matrix(), given_ind = numeric(), n_sim = 1000L, method = "default", ordering = "mean", log = FALSE, grad_lower = FALSE, grad_upper = FALSE, grad_sigma = FALSE, grad_given = FALSE, is_validation = TRUE, control = NULL, n_cores = 1L, marginal = NULL, grad_marginal = FALSE, grad_marginal_prob = FALSE ) ``` ## Arguments - `lower`: numeric vector representing lower integration limits. If `lower` is a matrix then each row determines new limits. Negative infinite values are allowed while positive infinite values are prohibited. - `upper`: numeric vector representing upper integration limits. If `upper` is a matrix then each row determines new limits. Positive infinite values are allowed while negative infinite values are prohibited. - `given_x`: numeric vector which `i`-th element corresponds to the given value of the `given_ind[i]`-th element (component) of multivariate normal vector. If `given_x` is numeric matrix then it's rows are such vectors of given values. - `mean`: numeric vector representing expectation of multivariate normal vector (distribution). - `sigma`: positively defined numeric matrix representing covariance matrix of multivariate normal vector (distribution). - `given_ind`: numeric vector representing indexes of multivariate normal vector which are conditioned at values given by `given_x` argument. - `n_sim`: positive integer representing the number of draws from Richtmyer sequence in GHK algorithm. More draws provide more accurate results by the cost of additional computational burden. Alternative types of sequences could be provided via `random_sequence` argument. - `method`: string representing the method to be used to calculate multivariate normal probabilities. Possible options are `"GHK"` and `"Gassmann"` recommended for high dimensional (more than 5) and low dimensional probabilities correspondingly. Additional option `"default"` selects optimal method depending on the number of dimensions. See 'Details' for additional information. - `ordering`: string representing the method to be used to order the integrals. See Details section below. - `log`: logical; if `TRUE` then probabilities (or densities) p are given as log(p) and derivatives will be given respect to log(p). - `grad_lower`: logical; if `TRUE` then the vector of partial derivatives of the probability will be calculated respect to each element of `lower`. If `lower` is a matrix then gradients will be estimated for each row of `lower`. - `grad_upper`: logical; if `TRUE` then the vector of partial derivatives of the probability will be calculated respect to each element of `upper`. If `upper` is a matrix then gradients will be estimated for each row of `upper`. - `grad_sigma`: logical; if `TRUE` then the vector of partial derivatives (gradient) of the probability will be calculated respect to each element of `sigma`. If `lower` and `upper` are matrices then gradients will be estimated for each row of these matrices. - `grad_given`: logical; if `TRUE` then the vector of partial derivatives of the density function will be calculated respect to each element of `given_x`. If `given_x` is a matrix then gradients will be estimated for each row of `given_x`. - `is_validation`: logical value indicating whether input arguments should be validated. Set it to `FALSE` to get performance boost (default value is `TRUE`). - `control`: a list of control parameters. See Details. - `n_cores`: positive integer representing the number of CPU cores used for parallel computing. Currently it is not recommended to set `n_cores > 1` if vectorized arguments include less then 100000 elements. - `marginal`: list such that `marginal[[i]]` represents parameters of marginal distribution of the `i`-th component of the random vector and `names(marginal)[i]` is a name of this distribution. If `names(marginal)[i] = "logistic"` or `names(marginal)[i] = "normal"` then `marginal[[i]]` should be an empty vector or `NULL`. If `names(marginal)[i] = "student"` then `marginal[[i]]` should contain a single element representing degrees of freedom. If `names(marginal)[i] = "PGN"` or `names(marginal)[i] = "hpa"` then `marginal[[i]]` should be a numeric vector which elements correspond to `pc` argument of `phpa0`. - `grad_marginal`: logical; if `TRUE` then the vector of partial derivatives (gradient) of probability will be calculated respect to each parameter of marginal distributions i.e. respect to each element of `marginal`. The gradient respect to the parameters of the `i`-th marginal distribution will be stored in the `grad_marginal[[i]]` output matrix which `j`-th column stores derivatives respect to `marginal[[i]][j]`. - `grad_marginal_prob`: logical; if `TRUE` then the vector of partial derivatives (gradient) of probability will be calculated respect to each cumulative marginal probability of marginal distributions. ## Returns This function returns an object of class "mnorm_pmnorm". An object of class "mnorm_pmnorm" is a list containing the following components: * `prob` - probability that multivariate normal random variable will be between `lower` and `upper` bounds. * `grad_lower` - gradient of probability respect to `lower` if `grad_lower` or `grad_sigma` input argument is set to `TRUE`. * `grad_upper` - gradient of probability respect to `upper` if `grad_upper` or `grad_sigma` input argument is set to `TRUE`. * `grad_sigma` - gradient respect to the elements of `sigma` if `grad_sigma` input argument is set to `TRUE`. * `grad_given` - gradient respect to the elements of `given_x` if `grad_given` input argument is set to `TRUE`. * `grad_marginal` - gradient respect to the elements of `marginal_par` if `grad_marginal` input argument is set to `TRUE`. Currently only derivatives respect to the parameters of `"PGN"` distribution are available. If `log` is `TRUE` then `prob` is a log-probability so output `grad_lower`, `grad_upper`, `grad_sigma` and `grad_given` are calculated respect to the log-probability. Output `grad_lower` and `grad_upper` are Jacobian matrices which rows are gradients of the probabilities calculated for each row of `lower` and `upper` correspondingly. Similarly `grad_given` is a Jacobian matrix respect to `given_x`. Output `grad_sigma` is a 3D array such that `grad_sigma[i, j, k]` is a partial derivative of the probability function respect to the `sigma[i, j]` estimated for `k`-th observation. Output `grad_marginal` is a list such that `grad_marginal[[i]]` is a Jacobian matrice which rows are gradients of the probabilities calculated for each row of `lower` and `upper` correspondingly respect to the elements of `marginal_par[[i]]`. ## Details Consider notations from the Details sections of `cmnorm` and `dmnorm`. The function calculates probabilities of the form: $$ P\left(x^{(l)}\leq X_{I_{d}}\leqx^{(u)}|X_{I_{g}}=x^{(g)}\right) $$ where $x^{(l)}$ and $x^{(u)}$ are lower and upper integration limits respectively i.e. `lower` and `upper` correspondingly. Also $x^{(g)}$ represents `given_x`. Note that `lower` and `upper` should be matrices of the same size. Also `given_x` should have the same number of rows as `lower` and `upper`. To calculate bivariate probabilities the function applies the method of Drezner and Wesolowsky described in A. Genz (2004). In contrast to the classical implementation of this method the function applies Gauss-Legendre quadrature with 30 sample points to approximate integral (1) of A. Genz (2004). Classical implementations of this method use up to 20 points but requires some additional transformations of (1). During preliminary testing it has been found that approach with 30 points provides similar accuracy being slightly faster because of better vectorization capabilities. To calculate trivariate probabilities the function uses Drezner method following formula (14) of A. Genz (2004). Similarly to bivariate case 30 points are used in Gauss-Legendre quadrature. The function may apply the method of Gassmann (2003) for estimation of $m>3$ dimensional normal probabilities. It uses matrix $5$ representation of Gassmann (2003) and 30 points of Gauss-Legendre quadrature. For $m$-variate probabilities, where $m > 1$, the function may apply GHK algorithm described in section 4.2 of A. Genz and F. Bretz (2009). The implementation of GHK is based on deterministic Halton sequence with `n_sim` draws and uses variable reordering suggested in section 4.1.3 of A. Genz and F. Bretz (2009). The ordering algorithm may be determined via `ordering` argument. Available options are `"NO"`, `"mean"` (default), and `"variance"`. Univariate probabilities are always calculated via standard approach so in this case `method` will not affect the output. If `method = "Gassmann"` then the function uses fast (aforementioned) algorithms for bivariate and trivariate probabilities or the method of Gassmann for $m>3$ dimensional probabilities. If `method = "GHK"` then GHK algorithm is used. If `method = "default"` then `"Gassmann"` is used for bivariate and trivariate probabilities while `"GHK"` is used for $m>3$ dimensional probabilities. During future updates `"Gassmann"` may become a default method for calculation of the $4-5$ dimensional probabilities. We are going to provide alternative estimation algorithms during future updates. These methods will be available via `method` argument. The function is optimized to perform much faster when all upper integration limits `upper` are finite while all lower integration limits `lower` are negative infinite. The derivatives could be also calculated much faster when some integration limits are infinite. For simplicity of notations further let's consider unconditioned probabilities. Derivatives respect to conditioned components are similar to those mentioned in Details section of `dmnorm`. We also provide formulas for $m \geq 3$. But the function may calculate derivatives for $m \leq 2$ using some simplifications of the formulas mentioned below. If `grad_upper` is `TRUE` then function additionally estimates the gradient respect to `upper`: $$ \frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial x^{(u)}_{i}}=P\left(x^{(l)}_{(-i)}\leq X_{(-i)}\leq x^{(u)}_{(-i)}|X_{i}=x^{(u)}_{i}\right)f_{X_{i}}\left(x^{(u)}_{i};\mu_{i},\Sigma_{i,i}\right) $$ If `grad_lower` is `TRUE` then function additionally estimates the gradient respect to `lower`: $$ \frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial x^{(l)}_{i}}=-P\left(x^{(l)}_{(-i)}\leq X_{(-i)}\leq x^{(u)}_{(-i)}|X_{i}=x^{(l)}_{i}\right)f_{X_{i}}\left(x^{(l)}_{i};\mu_{i},\Sigma_{i,i}\right) $$ If `grad_sigma` is `TRUE` then function additionally estimates the gradient respect to `sigma`. For $i\ne j$ the function calculates derivatives respect to the covariances: $$ \frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial \Sigma_{i, j}}= $$ $$ =P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}|X_{i}=x^{(u)}_{i}, X_{j}=x^{(u)}_{j}\right)f_{X_{i}, X_{j}}\left(x^{(u)}_{i}, x^{(u)}_{j};\mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) - $$ $$ -P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}|X_{i}=x^{(l)}_{i}, X_{j}=x^{(u)}_{j}\right)f_{X_{i}, X_{j}}\left(x^{(l)}_{i}, x^{(u)}_{j};\mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) - $$ $$ -P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}|X_{i}=x^{(u)}_{i}, X_{j}=x^{(l)}_{j}\right)f_{X_{i}, X_{j}}\left(x^{(u)}_{i}, x^{(l)}_{j};\mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) + $$ $$ +P\left(x^{(l)}_{(-(i, j))}\leq X_{-(i, j)}\leq x^{(u)}_{(-(i, j))}|X_{i}=x^{(l)}_{i}, X_{j}=x^{(l)}_{j}\right)f_{X_{i}, X_{j}}\left(x^{(l)}_{i}, x^{(l)}_{j};\mu_{(i,j)},\Sigma_{(i, j),(i, j)}\right) $$ Note that if some of integration limits are infinite then some elements of this equation converge to zero which highly simplifies the calculations. Derivatives respect to variances are calculated using derivatives respect to covariances and integration limits: $$ \frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial \Sigma_{i, i}}= $$ $$ -\frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial x^{(l)}_{i}} \frac{x^{(l)}_{i}}{2\Sigma_{i, i}}-\frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial x^{(u)}_{i}} \frac{x^{(u)}_{i}}{2\Sigma_{i, i}}- $$ $$ -\sum\limits_{j\ne i}\frac{\partial P\left(x^{(l)}\leq X\leqx^{(u)}\right)}{\partial \Sigma_{i, j}}\frac{\Sigma_{i, j}}{2\Sigma_{i, i}} $$ If `grad_given` is `TRUE` then function additionally estimates the gradient respect to `given_x` using formulas similar to those described in Details section of `dmnorm`. More details on abovementioned differentiation formulas could be found in the appendix of E. Kossova and B. Potanin (2018). If `marginal` is not empty then Gaussian copula will be used instead of a classical multivariate normal distribution. Without loss of generality let's assume that $\mu$ is a vector of zeros and introduce some additional notations: $$ q_{i}^{(u)} = \Phi^{-1}\left(P_{i}\left(\frac{x_{i}^{(u)}}{\sigma_{i}}\right)\right) $$ $$ q_{i}^{(l)} = \Phi^{-1}\left(P_{i}\left(\frac{x_{i}^{(l)}}{\sigma_{i}}\right)\right) $$ where $\Phi(.)^{-1}$ is a quantile function of a standard normal distribution and $P_{i}$ is a cumulative distribution function of the standartized (to zero mean and unit variance) marginal distribution which name and parameters are defined by `names(marginal)[i]` and `marginal[[i]]` correspondingly. For example if `marginal[i] = "logistic"` then: $$ P_{i}(t) = \frac{1}{1+e^{-\pi t / \sqrt{3}}} $$ Let's denote by $X^{*}$ random vector which is distributed with Gaussian (its covariance matrix is $\Sigma$) copula and marginals defined by `marginal`. Then probabilities for these random vector are calculated as follows: $$ P\left(x^{(l)}\leq X^{*}\leqx^{(u)}\right) =P\left(\sigma q^{(l)}\leq X\leq\sigma q^{(u)}\right) = P_{0}\left(\sigma q^{(l)}, \sigma q^{(u)}\right) $$ where $q^{(l)} = (q_{1}^{(l)},...,q_{m}^{(l)})$, $q^{(u)} = (q_{1}^{(u)},...,q_{m}^{(u)})$ and $\sigma = (\sqrt{\Sigma_{1, 1}},...,\sqrt{\Sigma_{m, m}})$. Therefore probabilities of $X^{*}$ may be calculated using probabilities of multivariate normal random vector $X$ (with the same covariance matrix) by substituting lower and upper integration limits $x^{(l)}$ and $x^{(u)}$ with $\sigma q^{(l)}$ and $\sigma q^{(u)}$ correspondingly. Therefore differentiation formulas are similar to those mentioned above and will be automatically adjusted if `marginal` is not empty (including conditional probabilities). Argument `control` is a list with the following input parameters: * `random_sequence` -- numeric matrix of uniform pseudo random numbers (like Halton sequence). The number of columns should be equal to the number of dimensions of multivariate random vector. If omitted than this matrix will be generated automatically using `n_sim` number of simulations. ## Examples ```r # Consider multivariate normal vector: # X = (X1, X2, X3, X4, X5) ~ N(mean, sigma) # Prepare multivariate normal vector parameters # expected value mean <- c(-2, -1, 0, 1, 2) n_dim <- length(mean) # correlation matrix cor <- c( 1, 0.1, 0.2, 0.3, 0.4, 0.1, 1, -0.1, -0.2, -0.3, 0.2, -0.1, 1, 0.3, 0.2, 0.3, -0.2, 0.3, 1, -0.05, 0.4, -0.3, 0.2, -0.05, 1) cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE) # covariance matrix sd_mat <- diag(c(1, 1.5, 2, 2.5, 3)) sigma <- sd_mat %*% cor %*% t(sd_mat) # Estimate probability: # P(-3 < X1 < 1, -2.5 < X2 < 1.5, -2 < X3 < 2, -1.5 < X4 < 2.5, -1 < X5 < 3) lower <- c(-3, -2.5, -2, -1.5, -1) upper <- c(1, 1.5, 2, 2.5, 3) p.list <- pmnorm(lower = lower, upper = upper, mean = mean, sigma = sigma) p <- p.list$prob print(p) # Additionally estimate a probability lower.1 <- c(-Inf, 0, -Inf, 1, -Inf) upper.1 <- c(Inf, Inf, 3, 4, 5) lower.mat <- rbind(lower, lower.1) upper.mat <- rbind(upper, upper.1) p.list.1 <- pmnorm(lower = lower.mat, upper = upper.mat, mean = mean, sigma = sigma) p.1 <- p.list.1$prob print(p.1) # Estimate the probabilities # P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4) lower.2 <- c(-1, -3, -5) upper.2 <- c(1, 3, 5) given_ind <- c(2, 4) given_x <- c(-2, 4) p.list.2 <- pmnorm(lower = lower.2, upper = upper.2, mean = mean, sigma = sigma, given_ind = given_ind, given_x = given_x) p.2 <- p.list.2$prob print(p.2) # Additionally estimate the probability # P(-Inf < X1 < 1, -3 < X3 < Inf, -Inf < X5 < Inf | X2 = 4, X4 = -2) lower.3 <- c(-Inf, -3, -Inf) upper.3 <- c(1, Inf, Inf) given_x.1 <- c(-2, 4) lower.mat.2 <- rbind(lower.2, lower.3) upper.mat.2 <- rbind(upper.2, upper.3) given_x.mat <- rbind(given_x, given_x.1) p.list.3 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2, mean = mean, sigma = sigma, given_ind = given_ind, given_x = given_x.mat) p.3 <- p.list.3$prob print(p.3) # Estimate the gradient of previous probabilities respect various arguments p.list.4 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2, mean = mean, sigma = sigma, given_ind = given_ind, given_x = given_x.mat, grad_lower = TRUE, grad_upper = TRUE, grad_sigma = TRUE, grad_given = TRUE) p.4 <- p.list.4$prob print(p.4) # Gradient respect to 'lower' grad_lower <- p.list.4$grad_lower # for the first probability print(grad_lower[1, ]) # for the second probability print(grad_lower[2, ]) # Gradient respect to 'upper' grad_upper <- p.list.4$grad_upper # for the first probability print(grad_upper[1, ]) # for the second probability print(grad_upper[2, ]) # Gradient respect to 'given_x' grad_given <- p.list.4$grad_given # for the first probability print(grad_given[1, ]) # for the second probability print(grad_given[2, ]) # Gradient respect to 'sigma' grad_given <- p.list.4$grad_given # for the first probability print(grad_given[1, ]) # for the second probability print(grad_given[2, ]) # Compare analytical gradients from the previous example with # their numeric (forward difference) analogues for the first probability n_dependent <- length(lower.2) n_given <- length(given_x) n_dim <- n_dependent + n_given delta <- 1e-6 grad_lower.num <- rep(NA, n_dependent) grad_upper.num <- rep(NA, n_dependent) grad_given.num <- rep(NA, n_given) grad_sigma.num <- matrix(NA, nrow = n_dim, ncol = n_dim) for (i in 1:n_dependent) { # respect to lower lower.delta <- lower.2 lower.delta[i] <- lower.2[i] + delta p.list.delta <- pmnorm(lower = lower.delta, upper = upper.2, given_x = given_x, mean = mean, sigma = sigma, given_ind = given_ind) grad_lower.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta # respect to upper upper.delta <- upper.2 upper.delta[i] <- upper.2[i] + delta p.list.delta <- pmnorm(lower = lower.2, upper = upper.delta, given_x = given_x, mean = mean, sigma = sigma, given_ind = given_ind) grad_upper.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta } for (i in 1:n_given) { # respect to lower given_x.delta <- given_x given_x.delta[i] <- given_x[i] + delta p.list.delta <- pmnorm(lower = lower.2, upper = upper.2, given_x = given_x.delta, mean = mean, sigma = sigma, given_ind = given_ind) grad_given.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta } for (i in 1:n_dim) { for(j in 1:n_dim) { # respect to sigma sigma.delta <- sigma sigma.delta[i, j] <- sigma[i, j] + delta sigma.delta[j, i] <- sigma[j, i] + delta p.list.delta <- pmnorm(lower = lower.2, upper = upper.2, given_x = given_x, mean = mean, sigma = sigma.delta, given_ind = given_ind) grad_sigma.num[i, j] <- (p.list.delta$prob - p.list.4$prob[1]) / delta } } # Comparison of gradients respect to lower integration limits h.lower <- cbind(analytical = p.list.4$grad_lower[1, ], numeric = grad_lower.num) print(h.lower) # Comparison of gradients respect to upper integration limits h.upper <- cbind(analytical = p.list.4$grad_upper[1, ], numeric = grad_upper.num) print(h.upper) # Comparison of gradients respect to given values h.given <- cbind(analytical = p.list.4$grad_given[1, ], numeric = grad_given.num) print(h.given) # Comparison of gradients respect to the covariance matrix h.sigma <- list(analytical = p.list.4$grad_sigma[, , 1], numeric = grad_sigma.num) print(h.sigma) # Let's again estimate probability # P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4) # But assume that standardized (to zero mean and unit variance): # 1) X1 and X2 have standardized PGN distribution with coefficients vectors # pc1 = (0.5, -0.2, 0.1) and pc2 = (0.2, 0.05) correspondingly. # 2) X3 has standardized student distribution with 5 degrees of freedom # 3) X4 has standardized logistic distribution # 4) X5 has standard normal distribution marginal <- list(PGN = c(0.5, -0.2, 0.1), hpa = c(0.2, 0.05), student = 5, logistic = numeric(), normal = NULL) p.list.5 <- pmnorm(lower = lower.2, upper = upper.2, mean = mean, sigma = sigma, given_ind = given_ind, given_x = given_x, grad_lower = TRUE, grad_upper = TRUE, grad_sigma = TRUE, grad_given = TRUE, marginal = marginal, grad_marginal = TRUE) # Lets investigate derivatives respect to parameters # of marginal distributions # respect to pc1 of X1 p.list.5$grad_marginal[[1]] # respect to pc2 of X2 p.list.5$grad_marginal[[2]] # derivative respect to degrees of freedom of X5 is # currently unavailable and will be set to 0 p.list.5$grad_marginal[[3]] ``` ## References Genz, A. (2004), Numerical computation of rectangular bivariate and trivariate normal and t-probabilities, Statistics and Computing, 14, 251-260. Genz, A. and Bretz, F. (2009), Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Vol. 195. Springer-Verlag, Heidelberg. E. Kossova, B. Potanin (2018). Heckman method and switching regression model multivariate generalization. Applied Econometrics, vol. 50, pages 114-143. H. I. Gassmann (2003). Multivariate Normal Probabilities: Implementing an Old Idea of Plackett's. Journal of Computational and Graphical Statistics, vol. 12 (3), pages 731-752.

pmnorm function