# Subtraction of two squared multivariate normal variables

Suppose we have a vector of random variables $X \sim N(\mu_0, \Sigma_0)$. We are interested in the distribution of $Y$ defined as follows:

$Y = (X - \mu_0)^T\Sigma_0^{-1}(X - \mu_0) - (X - \mu_1)^T\Sigma_1^{-1}(X - \mu_1)$

Note that the distribution of $(X - \mu_0)^T\Sigma_0^{-1}(X - \mu_0)$ is Chi-squared distribution. The question is that what is the distribution of $Y$.

Thank you.

• Since the two terms both involve the same $X$, they certainly won't be independent in general, so that's not something we can "note." Precisely what assumptions are you making about the relationship between $\Sigma_0$ and $\Sigma_1$?
– whuber
Commented Sep 17, 2015 at 15:05
• You are right. We can assume that $\Sigma_0$ and $\Sigma_1$ are independent (the same for $\mu_0$ and $\mu_1$). Commented Sep 17, 2015 at 15:11
• What does "independence" of $\Sigma_0$ and $\Sigma_1$ mean? Are these random matrices in some sense? And if so, how does one guarantee that they are invertible and that at least $\Sigma_0$ is positive semidefinite? Commented Sep 17, 2015 at 15:34
• By "independence", I just mean that there is no relationships between $\Sigma_0$ and $\Sigma_1$. $\Sigma_0$ and $\Sigma_1$ are invertible and positive semidefinite. Commented Sep 17, 2015 at 16:07

In case you wanted $Y$ to be defined as follows:

$$X_i \sim N(\mu_i, \Sigma_i) \; \text{for} \;i={0,1}$$ so that $$Y = (X_0-\mu_0)'\Sigma_0^{-1}(X_0-\mu_0) \; - \; (X_1-\mu_1)'\Sigma_1^{-1}(X_1-\mu_1)$$ where the subtrahend and minuend both are $\sim \chi^2(n)$. I therefore assume that both $X_i$ follow a $n$-variate normale distribution and that $X_0$ and $X_1$ are independent. Then the difference $Y$ follows the variance-gamma distribution. There was already a similar question to this problem in the math community on stackexchange:

Transformation of Difference of Random Variables

In R you can verfiy the result you find following the above link via the following code chunk:

library(MASS)
library(VarianceGamma)

sigma0 <- diag(4)
sigma1 <- diag(4)*2 + 1

X0 <- mvtnorm::rmvnorm(50000, mean=rep(0, 4), sigma0)
X1 <- mvtnorm::rmvnorm(50000, mean=rep(0, 4), sigma1)

X <- cbind(X0, X1)

Y <- apply(X, 1, function(z) { t(z[1:4])%*%solve(sigma0)%*%z[1:4] - t(z[5:8])%*%solve(sigma1)%*%z[5:8] } )

truehist(Y, nbins = 50, ylim=c(0,0.2))
lines(x=seq(-20,20,0.01)-4, y=dvg(x=seq(-20,20,0.01), vgC=4, sigma=4, theta=0, nu=0.5))


• The problem is that in my case, $X_0 = X_1$. Commented Sep 19, 2015 at 13:35
• Let's assume $\Sigma_0 = \Sigma_1$. By this assumption, I think that the problem becomes easier Commented Sep 22, 2015 at 8:10