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:


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.

  • $\begingroup$ R demands everything reside in memory. Likely that it is memory issues that are causing the problem. $\endgroup$
    – meh
    Oct 26, 2017 at 13:00
  • $\begingroup$ What is the full code you are currently using? $\endgroup$ Oct 26, 2017 at 13:02
  • $\begingroup$ Hi @Greenparker. I included the R-code I am using. I will appreciate any suggestion. $\endgroup$ Oct 26, 2017 at 13:30
  • $\begingroup$ Thanks, @aginensky. I think so. However, I tried to find (with no success) any algorithm to make more efficient computations. $\endgroup$ Oct 26, 2017 at 13:31
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 26, 2017 at 14:01

1 Answer 1


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$.

  • $\begingroup$ Thanks, @AdamO, for sharing your perspective of the problem. $\endgroup$ Oct 26, 2017 at 16:44
  • $\begingroup$ @Student1981 I admit i haven't done much research in the area. However, I feel confident that it is an open problem to find more efficient multivariate normal implementations, that is if the manuscript I linked is to be trusted. Not a perspective, at least not in the subjective sense. $\endgroup$
    – AdamO
    Oct 26, 2017 at 21:01
  • $\begingroup$ One might even go further and say that the whole concept of "efficiency" here is about the rate at which the computation becomes longer as $n$ increases (i.e., we can't even say that it is less efficient as $n$ increases, so much as we can saay it takes longer as $n$ increases). $\endgroup$
    – Ben
    Apr 20, 2020 at 5:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.