41

Ok, so the Chi-Squared distribution with $n$ degrees of freedom is the sum of the squares of $n$ independent Gaussian random variables.

The trouble is, my Gaussian random variables are not independent. They do however all have zero mean and the same variance. Supposing I have a covariance matrix---which again is not a diagonal matrix because they aren't independent, but all the elements along the diagonal are equal to each other because they have the same variance, and in fact the covariance matrix is a symmetric Toeplitz matrix (and I'm not saying that this is important to the solution if there is one, but if it's a necessary property to get anywhere, by all means use that fact)---is there some way to decompose this sum of squares of these Gaussian random variables into perhaps a sum of chi-squared random variables and possibly Gaussian random variables? In other words, I can't directly just square them all and add them together and call it a chi squared distribution because a chi squared distribution is a sum of independent Gaussian squares, and they aren't independent.

I know how to find a linear transformation of the Gaussian random variables which are $n$ independent Gaussians, but that's no help because they aren't the things being squared, you see.

kjetil b halvorsen
  • 6,625
zortharg
  • 682

1 Answers1

54

Lets assume you have $X=(X_1, \dots, X_n)$ a random vector with multinormal distribution with expectation vector $\mu$ and covariance matrix $\Sigma$. We are interested in the quadratic form $Q(X)= X^T A X = \sum \sum a_{ij} X_i X_j$. Define $Y = \Sigma^{-1/2} X$ where we are assuming $\Sigma$ is invertible. Write also $Z=(Y-\Sigma^{-1/2} \mu)$, which will have expectation zero and variance matrix the identity.

Now $$ Q(X) = X^T A X= (Z+\Sigma^{-1/2} \mu)^T \Sigma^{1/2}A\Sigma^{1/2} (Z+\Sigma^{-1/2} \mu). $$ Use the spectral theorem now and write $\Sigma^{1/2}A \Sigma^{1/2} = P^T \Lambda P$ where $P$ is an orthogonal matrix (so that$P P^T=P^T P=I$) and $\Lambda$ is diagonal with positive diagonal elements $\lambda_1, \dotsc, \lambda_n$. Write $U = P Z$ so that $U$ is multivariate normal with identity covariance matrix and expectation zero.

We can compute $$ Q(X) = (Z+\Sigma^{-1/2} \mu)^T \Sigma^{1/2}A\Sigma^{1/2} (Z+\Sigma^{-1/2} \mu) \\ = (Z+\Sigma^{-1/2} \mu)^T P^T \Lambda P (Z+\Sigma^{-1/2} \mu) \\ = (PZ+P\Sigma^{-1/2} \mu)^T \Lambda (PZ+P\Sigma^{-1/2} \mu) \\ = (U+b)^T \Lambda (U+b) $$ where now $b = P \Sigma^{-1/2} \mu $. (There was a small typo in above defs of $U$ and $b$, now corrected.) So: $$ Q(X) = X^T A X = \sum_{j=1}^n \lambda_j (U_j+b_j)^2 $$ In your case, $\mu=0$ so $b=0$ so your quadratic form is a linear combination of independent chi-square variables, each with one degree of freedom. In the general case, we will get a linear combination of independent non-central chisquare variables.

If you want to work numerically with that distribution, there is an CRAN package (that is, package for R) implementing it, called CompQuadForm.

If you want (much) more detail, there is a book dedicated to the topic, Mathai & Provost: "Quadratic forms in random variables".

kjetil b halvorsen
  • 6,625
  • 2
    Wow, you are awesome. I had pretty much given up on anyone answering this. – zortharg Jul 15 '13 at 03:14
  • 1
    What it actually is that I'm trying to do, is find the variance of the sum of the square change in the value of a random walk after a time of t over a certain range of time. So for instance, the average square change after 3 time units if there are 7 samples is ((x(3)-x(0))^2+(x(4)-x(1))^2+(x(5)-x(2))^2+(x(6)-x(3))^2) / 4 - and if the random walk has a variance of 1 per time unit, then the expected value of this sum is 3 - but what is its variance? And I THOUGHT I answered it yesterday using induction, but the results aren't matching the results of my Monte Carlo simulations. Any ideas? – zortharg Jul 15 '13 at 03:26
  • First a fast question? How did you do the empirical estimation? As you present it, it seems like you did the different sums on the same random walk, but then they will not be independent replications, so you have no reason to expect they will converge on the expected value. – kjetil b halvorsen May 09 '14 at 12:35
  • 1
    How does $(U+b)^T P^T \Sigma^{1/2}A \Sigma^{1/2} P (U+b)$ simplify to $\Lambda$? If I solve that I end up with $P^T P^T \Lambda P P$ but how do I know that this is equal to $\Lambda$? – drahnr Feb 08 '15 at 10:15
  • 1
    http://math.stackexchange.com/questions/1138999/orthogonal-matrices-and-diagonal-matrices-multiplied – drahnr Feb 08 '15 at 23:07
  • 2
    This is only true for input distributions with equal std dev. – drahnr Feb 11 '15 at 11:44
  • @drahnr: The question by the OP specified equal variance, so that is assumed. – kjetil b halvorsen Feb 11 '15 at 18:09
  • @kjetilbhalvorsen do you know by chance a way to get the distribution for a case with different variances? – drahnr Feb 11 '15 at 23:32
  • 2
    @drahnr: Looking closer at this, there is a small typo I am correcting now, but the development is completely general and works for any invertible $\Sigma$. – kjetil b halvorsen Feb 13 '15 at 12:42
  • How do we know the spectral theorem holds for all $Σ^{.5}Σ^{.5}$? @kjetilbhalvorsen – swedishfished Jun 01 '20 at 01:10
  • @swedishfished: Because it is symmetric. $A$ is symmetric, $\Sigma^{1/2}$ is symmetric square root of symmetric matrix. – kjetil b halvorsen Jun 01 '20 at 14:39
  • 1
    @kjetilbhalvorsen Does the result holds if $\Sigma$ is singular? – Alphie Jan 11 '23 at 14:23