# Efficiency of the command pmvnorm in R

Let $X_1,\,X_2,\ldots,\,X_n$ be $n$ independent random variables, where $X_i\sim\text{N}(\mu_i,\,\sigma_i^2)$, for all $i\in\{1,\ldots,\,n\}$. Consider the $n$-dimensional random vector
$$\boldsymbol{Y}_n=\left(X_1,\,X_1+X_2,\ldots,\,X_1+\cdots+X_n\right)^T$$ I learned that $\boldsymbol{Y}_n$ follows a $n$-dimensional normal distribution, with $$\mathbb{E}(\boldsymbol{Y}_n)=\left(\mu_1,\,\mu_1+\mu_2,\ldots,\,\mu_1+\cdots+\mu_n\right)^T$$ and the covariance matrix of $\boldsymbol{Y}_n$ is a $n\times n$ matrix $[a_{i,\,j}]$ such that $$a_{i,\,j}=\sum_{k=1}^{\min\{i,\,j\}}\sigma^2_k$$ When computing the joint cumulative distribution function of $\boldsymbol{Y}_n$ in ${\tt R}$-software, I realized that the command ${\tt pmvnorm}$ becomes less efficient as $n$ increases (e.g., when $n=3000$).

The following is the R-code I am using:

library(mvtnorm)

n <- 3000
Mean_X <- runif(n, min = -20, max = 20)
Var_X <- rexp(n, rate = 50)
CumSum_Var_X <- cumsum(Var_X)

Mean_Y <- cumsum(Mean_X)
Var_Y <- matrix(Var_X[1], nrow = 1, ncol = n)
for(k in 2:n) Var_Y <- rbind(Var_Y, c(CumSum_Var_X[1:k], rep(CumSum_Var_X[k], n - k)))

y <- runif(n, min = -20, max = 20)
# Joint cumulative distribution function of Yn, evaluated in y.
Pr <- pmvnorm(upper = y, mean = Mean_Y, sigma = Var_Y)


Question: I was wondering if you could tell me how to compute the joint cumulative distribution function of $\boldsymbol{Y}_n$ in ${\tt R}$-software in a more efficient way than ${\tt pmvnorm}$ for large values of $n$.

Thanks a lot for your help.

• R demands everything reside in memory. Likely that it is memory issues that are causing the problem. – meh Oct 26 '17 at 13:00
• What is the full code you are currently using? – Greenparker Oct 26 '17 at 13:02
• Hi @Greenparker. I included the R-code I am using. I will appreciate any suggestion. – Student1981 Oct 26 '17 at 13:30
• Thanks, @aginensky. I think so. However, I tried to find (with no success) any algorithm to make more efficient computations. – Student1981 Oct 26 '17 at 13:31
• You can view $(Y_m)_{m\le n}$ as a realisation of an inhomogeneous Markov chain (or time changed Brownian motion if you prefer). Use the Markov property to recursively calculate the path probability without computing the whole covariance matrix. – P.Windridge Oct 26 '17 at 14:01

In general, any R operation will become less efficient as $n$ increases. Especially problematic is that we can't immediately assume any convenient representation of these data, like independence, so the covariance matrix must be supplied explicitly in the form of an $n \times n$ matrix. Complicating matters further is that you are basically calling for the evaluation of an $n$-dimensional integral.
Inspecting the workhorse code for pmvnorm reveals a bit more insight to the problem. Much work has gone into developing better algorithms. See here from the package author. The approach used here is an adaptation of MCMC integration by using what he calls "quasi" random samples. It is an iterative, sampling based procedure that has an unproven complexity, but some intuition can be gained from Figure 3.
To improve computation, pre-process your data. If $\Sigma$ is block diagonal or approximately block diagonal, parse the problem into multiple smaller problems with smaller $n$.