Multivariate gaussian integral over positive reals The multivariate gaussian integral over the whole $\mathbf{R}^n$ has closed form solution
$$P = \int_{\mathbf{x} \in \mathbf{R}^n} \exp \left(-\frac12 \mathbf{x}^T \mathbf{A} \mathbf{x}\right)\,d\mathbf{x} = \sqrt{\frac{(2\pi)^n}{\det \mathbf{A}}}$$
where $\mathbf{A}$ is a symmetric positive-definite covariance matrix.
However, I need to solve the integral for positive reals $\{\mathbf{x} \in \mathbf{R}^n :\, \mathbf{x}_i \geq 0\ \forall i\}$ only and in at least 6 dimensions:
$$P = \int_{\{\mathbf{x} \in \mathbf{R}^n :\, \mathbf{x}_i \geq 0\ \forall i\}} \exp \left(-\frac12 \mathbf{x}^T \mathbf{A} \mathbf{x}\right)\,d\mathbf{x}$$
For diagonal $\mathbf{A}$ with zero covariance, a solution has been published.
For non-diagonal covariance, my approach so far is to apply affine coordinate transforms to rotate and rescale the gaussian ellipsoid into the unit sphere (see here).
In two dimensions, the solution to the integral then reduces to comparing the area enclosed by the transformed positive coordinate axes (blue) to the area of the unit circle:

In three dimensions, the solution is given by the ration of the surface area of an enclosed spherical polygon to the surface area of the unit sphere.
In four dimensions, this approach becomes quite complicated, and I don't know how to use the usual spherical excess formulas for higher dimensions.
Any ideas or alternative approaches? Is there a multivariate error function? Any treatment on the multivariate half normal distribution?

Addition (2018-12-03):
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:


*

*Completing the square in one variable, say $x$, leaves me with
$$\int_{\mathbb{R}_+^2} \mathrm{d}y\mathrm{d}z \exp\left(-\frac{1}{2} \frac{\mathrm{det}\,A_3}{\mathrm{det}\,A_2}z^2\right) \exp\left(-\frac{1}{2} \frac{\mathrm{det}\, A_2}{a}(y-m z)^2\right) \left[1 - \mathrm{erf}\left(\frac{a_{12}y+a_{13}z}{\sqrt{2a}}\right) \right] $$
where $A_2=\begin{pmatrix} a & a_{12}\\ & b\end{pmatrix}$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_{12}y+a_{13}z$ or something similar.

*Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral. 
Any help / hint would be greatly appreciated! Thank you in advance.
 A: The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:


*

*DR Cox, N Wermuth, Biometrika, 2002

*JH Macke, P Berens et al., Neural Computation, 2009

A: Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
A: For Przemo's example with $n=5$, define a higher arcsine:
$$\arcsin(a,b,c) = \int_0^b \frac { \arcsin\left( a \beta c / v / w \right) } { \sqrt{1-\beta^2} } d\beta$$
$$\text{where}\; v = \sqrt{ 1-a^2-\beta^2 },\; w = \sqrt{ 1-\beta^2-c^2 },$$
and the orthoscheme volume:
$$V( a, b, c, d ) = \frac {\pi^2} {60}
  + \frac {\pi} {30} \left( \arcsin( a ) + \arcsin( b ) + \arcsin( c ) + \arcsin( d ) \right)
  + \frac {1} {15} \left(
    \arcsin( a ) \arcsin( c ) + \arcsin( a ) \arcsin( d ) + \arcsin( b ) \arcsin( d )
    + \arcsin( a, b, c ) + \arcsin( b, c, d ) \right)
.$$
Then the (normalized) orthant probability for A is $\frac{15S}{8\pi^2}$ where:
$$ S = 8 \cdot V \left( abc, -\sqrt{\frac{1-a}{2}}, -\frac{1+a-2ab}{2(1-ab)}, -\frac{1-a}{2(1-ab)} \right) \\
+ 16 \cdot V \left( abc, -\sqrt{\frac{1-ab}{2}}, -\frac{1}{2}\sqrt{\frac{1-a}{1-ab}}, -\frac{1}{2}\sqrt{\frac{1+a-2ab}{1-ab}} \right)
. $$
E.g., for $a = 0.3, b = 0.2, c = 0.1$, we get $\approx 0.049972036321177378010658463826149025346$.
A: There is a solution in the form of a Taylor series.
See this Wikipedia page which refers to "Measuring Solid Angles Beyond Dimension Three" by Jason M. Ribando (DOI: 10.1007/s00454-006-1253-4).
For that Taylor series you get (Lemma 2.1 from Ribando) for a matrix $\mathbf{A} = (\alpha_{ij})$
$$\exp \left(-\frac12 \mathbf{x}^T \mathbf{A} \mathbf{x}\right) = \sum_{\mathbf{a} \in N^{{n}\choose{2}}} \left[ \frac{(-2)^{\sum_{i<j} a_{ij}}}{\prod_{i<j} a_{ij}!} \prod_{l=1}^n x_{l}^{\sum_{m \neq l} a_{lm}} \prod_{i} e^{-x_{i}^2} \right] \boldsymbol{\alpha}^{\mathbf{a}}$$
Where the subscript $i<j$ means that we sum over one half of the off-diagonal elements, and the subscript $m \neq l$ means that we sum over all cases in this half of the off-diagonal elements where either the first or the second subscript is equal to $l$ (This was not intuitive to me at first but say $l=3$ then the set $a_{31}, a_{32}, a_{34}, a_{35}, \dots a_{3n}$ is equal to $a_{13}, a_{23}, a_{34}, a_{35}, \dots a_{3n}$ because $a_{ij} = a_{ji}$ due to the symmetry).
If we integrate this then we can bring the integral inside the sum and express integrals as Gamma functions $\int_0^\infty x_i^k e^{-x_i^2} dx = \frac{\Gamma((k+1)/2) }{2}$
$$\begin{array}{} \int_{\mathbf{x} \in \mathbb{R}^n_{\geq 0}}\exp \left(-\frac12 \mathbf{x}^T \mathbf{A} \mathbf{x}\right) d\mathbf{x} &=& \int_{\mathbf{x} \in \mathbb{R}^n_{\geq 0}} \sum_{\mathbf{a} \in N^{{n}\choose{2}}} \left[ \frac{(-2)^{\sum_{i<j} a_{ij}}}{\prod_{i<j} a_{ij}!} \prod_{l=1}^n x_{l}^{\sum_{m \neq l} a_{lm}} \prod_{i} e^{-x_{i}^2} \right] \boldsymbol{\alpha}^{\mathbf{a}} d\mathbf{x} \\
 &=&  \sum_{\mathbf{a} \in N^{{n}\choose{2}}} \left[ \frac{(-2)^{\sum_{i<j} a_{ij}}}{\prod_{i<j} a_{ij}!} \prod_{i} \int_{x_i \in \mathbb{R}_{\geq 0}} x_{i}^{\sum_{m \neq l} a_{lm}} e^{-x_{i}^2} dx \right] \boldsymbol{\alpha}^{\mathbf{a}}\\
&=&  \sum_{\mathbf{a} \in N^{{n}\choose{2}}} \left[ \frac{(-2)^{\sum_{i<j} a_{ij}}}{\prod_{i<j} a_{ij}!} \prod_{i} \Gamma\left( \frac{1+\sum_{m \neq l} a_{lm}}{2} \right)  \right] \boldsymbol{\alpha}^{\mathbf{a}}\\
\end{array}$$
The R-code below demonstrates this approach and compares it with a Monte Carlo approach to estimate the value. The Monte Carlo approach considers a multivariate normal distribution which has a density of the form ${\det(2 \pi \boldsymbol{\Sigma})}^{-0.5} \exp \left(-\frac12 \mathbf{x}^T \boldsymbol{\Sigma}^{-1} \mathbf{x}\right)$ .
library(MASS)
set.seed(1)

###
### simulations to compute orthant probability
### 
P_mc = function(sig, n = 10^3) {
  dim = length(sig[1,])
  mu = rep(0,dim)
  X = mvrnorm(n, mu, sig)   ### sample from multivariate normal
  Y = (X > 0)               
  return(mean(rowSums(Y) == dim))  ### return the mean number of cases when all X>0
}

###
### taylor series to compute orthant probability
###

P_comp = function(sig, order = 3) {
  k_t = order + 1
  dim = length(sig[1,])  
  
  invsig = solve(sig) ### use inverse of Sig to relate with covariance marix of Gaussian used in P_mc
  
  ### rescale vectors such that diagonal is equal to 1
  ### this does not change the orthant
  f = diag(invsig) %*% t(diag(invsig))
  invsig = invsig / sqrt(f)
  
  ### create a vector with elements from the matrix with condition on subscripts i<j
  alpha = invsig[lower.tri(invsig)]

  ### create a matrix with each row containing the combination of powers for the multivariate taylor series
  l = length(alpha)
  a_base = c(1:k_t) - 1
  a = a_base
  for (k in 2:l) {
    a = cbind(rep(1, k_t) %x% a, rep(a_base, each = k_t^(k-1)))
  }
  n_rows = length(a[,1])

  ### compute the Taylor series terms
  sum_term = sapply(1:n_rows, FUN = function(j) {
    aj = a[j,]
    
    ### making a matrix with the power terms we use this to compute term_i
    Maj = matrix(rep(1,dim^2),dim)
    Maj[lower.tri(Maj) == 1] = aj
    Maj = t(Maj)
    Maj[lower.tri(Maj) == 1] = aj 
    
    prod1 = prod(factorial(aj))
    term_i = sapply(1:dim, FUN = function(i) {gamma( 0.5 * (sum(Maj[i,-i])+1) )})
    prod2 = prod(term_i)
    prod3 = prod(alpha^aj)
    pow = (-2)^sum(aj)
    pow/prod1*prod2*prod3 
  })
  
  value = sum(sum_term) * det(invsig)^0.5 / (4*pi)^(dim/2)
  
  return(value)
}

### some matrix to try out
set.seed(1)
rho = 0.3
sig = matrix(c(1,rho,rho,rho,
               rho,1,rho,rho,
               rho,rho,1,rho,
               rho,rho,rho,1),4)

P_mc(sig = sig, n  =10^6)
### 0.140251
P_comp(sig, order = 4)
### 0.1400075


sig = matrix(c(1,  -0.3, 0.1, 0.2,
              -0.3, 1,   0.2, 0.1,
               0.1, 0.2, 1,   0.4,
               0.2, 0.1, 0.4, 1), 4)
P_mc(sig = sig, n  =10^6)
### 0.088968
P_comp(sig, order = 4)
### 0.08863386

