Distribution of the magnitude of samples from the multivariate gaussian distribution

Let $\mathbf{X}\sim \mathcal{N}(0,I)$ be an $n-$dimensional standard Gaussian vector. We know that \begin{align} ||\mathbf{X}||^2&=\sum_{i=1}^nX_i^2\sim\chi^2_n\\ ||\mathbf{X}||&=\sqrt{\sum_{i=1}^nX_i^2}\sim\chi_n\\ \end{align}

Suppose the components are correlated, i.e., we consider $\mathbf{Y}\sim \mathcal{N}(0,\Sigma)$ where $\Sigma$ is positive definite. We know that $||\mathbf{Y}||^2$ has the generalized chi-squared distribution. Is there also a corresponding generalized chi distribution for $||\mathbf{Y}||$? If so, is there some asymptotic result on the mean of such a distribution for $n\to\infty$? This is related to

Expectation of square root of sum of independent squared uniform random variables

Central Limit Theorem for square roots of sums of i.i.d. random variables

But here the components of the random vector are neither independent or identically distributed.

• The method I used in the (first) linked post do not really depend on independence, it can be used directly. You just need to find the mgf of $\| \mathbf{ Y} \|^2$ – kjetil b halvorsen Dec 11 '17 at 17:16
• @kjetilbhalvorsen good point, but unfortunately cfs and mgfs weren't part of my (only) academic course on Probability, so I'm not too familiar with them. I know that the mgf of the sum of $n$ iid RVs with common mgf $f(t)$ is just $f(t)^n$ (because convolution becomes multiplication under a Fourier transform), but I know little else. Can you point me to some references which would help me calculate the mgf of correlated RVs? – DeltaIV Dec 12 '17 at 14:20
• You can start from my answer here: math.stackexchange.com/questions/442472/… or maybe cpaindex.naturalspublishing.com/files/published/…. Also, the book referenced in the first link do have the results. – kjetil b halvorsen Dec 12 '17 at 20:28
• Now I have written out the mgf as an answer here: stats.stackexchange.com/questions/262604/… – kjetil b halvorsen Dec 14 '17 at 20:48
• There is much more freedom here than in your former question, which only depended on $n$. Here there are $n$ eigenvalues, which might be all equal (easy case), or all different, or something in between. If you could indicate something about how $n$ tends to infinity (how do $\Sigma$ look like? – kjetil b halvorsen Dec 19 '17 at 21:36

This is not a complete answer so far, I will come back completing it when time. Using results and notation from What is the moment generating function of the generalized (multivariate) chi-square distribution?, we can write $$\| Y \|^2 = Y^T Y$$ which we can write as a sum of $n$ independent scaled chisquare random variables (here it is central chisquare since the expectation of $Y$ is zero) $Y^T Y = \sum_1^n \lambda_j U_j^2$ with the following definitions taken from the link above. Use the spectral theorem to write $\Sigma ^{1/2} A \Sigma^{1/2} = P^T \Lambda P$ with $P$ orthogonal and $\Lambda$ diagonal with diagonal elements $\lambda_j$, $U=P Y$ then have independent standard normal components. Then we find the mgf of $Y^T Y$ as $$M_{Y^T Y}(t) = \exp\left(-\frac12 \sum_1^n \log(1-2\lambda_j t) \right)$$ Then one can continue finding the expectation using the method from the link in comments. The expectation of $\sqrt{Y^T Y}$ is given by $$\DeclareMathOperator{\E}{\mathbb{E}} D^{1/2} M (0) = \Gamma(1/2)^{-1}\int_{-\infty}^0 |z|^{-1/2} M'(z) \; dz$$ where $M(t)$ is the mgf above and $'$ denotes differentiation. Using maple we can calculate

M := (t,lambda1,lambda2,lambda3) -> exp( -(1/2)*( log(1-2*t*lambda1)+log(1-2*t*lambda2)+log(1-2*t*lambda3) ) )

which is the case with $n=3$. As an example, if the covariance matrix $\Sigma$ is the equicorrelation matrix with offdiagonal elements $1/2$, that is, $$\Sigma=\begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & 1 & 1/2 \\ 1/2 & 1/2 & 1 \end{pmatrix}$$ then the eigenvalues $\lambda_j$ can be calculated as $2, 1/2, 1/2$. Then we can calculate

int( diff(M(z,2,1/2,1/2),z)/sqrt(abs(z)),z=-infinity..0 )/GAMMA(1/2)

obtaining the result. $$1/3\,{\frac {\sqrt {3}{\rm arctanh} \left(\sqrt {3}/2\right)+6}{ \sqrt {\pi}}}$$ It is probably better to just use the option numeric=true, we will not find some simple general symbolic answer. An example with other eigenvalues gives another symbolic form of the result:

int( diff(M(z,5,1,1/2),z)/sqrt(abs(z)),z=-infinity..0 )/GAMMA(1/2)

$$-1/15\,{\frac {\sqrt {5} \left( 10\,i\sqrt {5}{\it EllipticF} \left( 3 \,i,1/3 \right) -9\,i\sqrt {5}{\it EllipticE} \left( 3\,i,1/3 \right) -30 \right) }{\sqrt {\pi}}}$$ (I now also suspect that at least one of this results returned by maple is in error, see https://www.mapleprimes.com/questions/223567-Error-In-A-Complicated-Integralcomplex?sq=223567)

Now I will turn to numerical methods, programmed in R. First, one idea is to use the saddlepoint approximation or one of the other methods from the answer here: Generic sum of Gamma random variables combined with the relation $f_{\text{chi}}(x)=2xf_{\text{chi}^2}(x^2)$, but first I will look into use of the R package CompQuadForm, which implements several methods for the cumulative distribution function, and then numerical differentiation with the package numDeriv, code follows:

library(numDeriv)

make_pchi2   <-  function(lambda, h=rep(1, length(lambda)), delta=rep(0, length(lambda)), sigma=0, ...)  {
function(q) {
vals <- rep(0, length(q))
warning_indicator <- 0
for (j in seq(along=q)) {
res <-  davies(q[j], lambda=lambda, h=h, delta=delta, sigma=sigma, ...)
if (res$ifault > 0) warning_indicator <- warning_indicator+1 vals[j] <- res$Qq
}
if (warning_indicator > 0) warning("warnings reported", warning_indicator)
vals  <-  1-vals  # converting to cumulative probabilities
return(vals)
}
}

pchi2  <-  make_pchi2(c(2, 0.5), c(1, 2))

# Construction a density function via numerical differentiation:

make_dchi  <-  function(pchi) {
function(x, method="simple")  {
side  <-  ifelse(x==0, +1, NA)
vals  <-  rep(0,  length(x))
for (j in seq(along=x)) {
vals[j]  <-  grad(pchi,  x[j]^2, method=method, side=side[j])
}
vals  <-  2*x*vals
return(vals)
}
}

dchi  <-  make_dchi(pchi2)

plot( function(x) dchi(x, method="Richardson"),  from=0,  to=7, main="density of generalized chi distribution", xlab="x", ylab="density", col="blue")

which gives this result: and we can test with:

integrate(function(x) dchi(x,  method="Richardson"), lower=0,  upper=Inf)
0.9999565 with absolute error < 0.00011

Back to (approximation) for the expected value. I will try to get some approximation for the moment generating function of $Y^T Y$ for large $n$. Start with the expression for the mgf $$M_{Y^T Y}(t) = \exp\left(-\frac12 \sum_1^n \log(1-2\lambda_j t) \right)$$ found above. We concentrate on its logarithm, the cumulant generating function (cgf). Assume that when $n$ increases, then the eigenvalues will approximatly come from some distribution with density $g(\lambda)$ on $(0,\infty)$ (in the numerical examples we will just use the uniform distribution on $(0,1)$). Then the sum above divided into $n$ will be a Riemann sum for the integral $$I(t) = \int_0^\infty g(\lambda) \log(1-2\lambda t)\; d\lambda$$ for the uniform distribution example we will find $$I(t) = \frac{2t\ln(1-2t)-\ln(1-2t)-2t}{2t}$$ Then the mgf will be (approximately) $M_n(t)=e^{-\frac{n}2 I(t)}$ and the expected value of $\| Y \|$ is $$\Gamma(1/2)^{-1} \int_{-\infty}^0 |z|^{-1/2} M_n'(z) \; dz$$ This seems to be a difficult integral, some work seems to show that the integrand has a vertical asymptote at zero, for instance maple refuses to do a numerical integration (but uses abnormally long time). So it might be that a better strategy is to use numerical integration using the density approximations above. Nevertheless, this approach seems to indicate that the dependence on $n$ should not be very dependent on the eigenvalue distribution.

Going for numerical integration in R, that also shows to be difficult. First some code:

E  <-  function(lambda, h=rep(1, length(lambda)), lower, upper, ...) {
pchi2  <-  make_pchi2(lambda, h, acc=0.001, lim=50000)
dchi  <-   make_dchi(pchi2)
val  <-   integrate(function(x) x*dchi(x), lower=lower,   upper=upper, subdivisions=10000L, stop.on.error=FALSE)$value return(val) } Some experimentation with this code, using the integration limits lower=-Inf, upper=0, shows that for a large number of eigenvalues, the result is 0! The reason is that the integrand is zero "almost everywhere", and the integration routine misses the dump. But plotting the integrand first, and choosing a sensible interval, we can get reasonable results. For the example I am using uniformly distributed eigenvalues between 1 and 5, with varying$n$. So for$n=5000$the following code: E(seq(from=1, to=5, length.out=5000), lower=-130, upper=-110) works. This way we can compute the following results: n Expected value Approximation [1,] 5 3.656979 3.890758 [2,] 15 6.581486 6.738991 [3,] 25 8.559964 8.700000 [4,] 35 10.161615 10.293979 [5,] 45 11.544113 11.672275 [6,] 55 12.777101 12.904185 [7,] 65 13.901465 14.028328 [8,] 75 14.941275 15.068842 [9,] 85 15.915722 16.042007 [10,] 95 16.830085 16.959422 [11,] 105 17.697715 17.829694 [12,] 500 38.702620 38.907583 [13,] 1000 54.749000 55.023631 [14,] 5000 122.449200 123.036580 Where the last column is calculated as$1.74 \sqrt{n}$, and is found by regression analysis. We show this as a plot: The fit is quite good, indicating that a reasonable guess is that the expectation grows with the square root of$n$, with a prefactor which probably depends on the eigenvalue distribution. • Nice! I look forward to reading the rest of the answer. Thanks a lot! The complete answer will definitely deserve a bounty :-) – DeltaIV Dec 15 '17 at 10:37 • Thanks! In the meantime, the mgf can be used to approximate the density function of$Y^T Y$via the saddflepoint approximation, see stats.stackexchange.com/questions/72479/… (since chisquare is a special case of gamma). If that idea can be extended to work also for its squareroot, I do not know. – kjetil b halvorsen Dec 15 '17 at 16:35 • Could we not simplify your moment generating function into the following? $$M_{Y^2}(t)=\prod_1^n (1-2\lambda_jt)^{-\frac{1}{2}}$$ – Martijn Weterings Dec 15 '17 at 21:41 • Surely, but I am not sure that is a simplification---for instance, I had better results with maple using the exponential form. The exponential form is also often better for numerical computation – kjetil b halvorsen Dec 15 '17 at 22:03 • @kjetilbhalvorsen ok, if I understand correctly a generic answer is not possible, it depends on the variance-covariance matrix. That's fine. Given that we cannot find an explicit expression, do you think it could be possible at least to get an estimate for the growth of$||\mathbf{Y}||$with$n$? If I could at least say that, even with correlated components, the magnitude of the multivariate Gaussian vector goes to$\infty$as fast as, say,$\sqrt{n}$, or maybe$\Gamma(n+1/2)/\Gamma(n)\$, that would be a nice result. – DeltaIV Dec 18 '17 at 11:59