# The multivariate normal distribution and its calculation

Let a random vector $$X = (X_1, X_2,\ldots, X_p)^{\mathrm{T}}$$ has the multivariate normal distribution with mean vector $$\mu_1$$ and covariance matrix $$\Sigma_1 > 0$$, i.e. $$X \sim \mathcal{N}_p(\mu_1, \Sigma_1)$$, and has the density $$f_1(x)$$, $$x \in \mathbb{R}^p$$. Let a random vector $$Y = (Y_1, Y_2,\ldots, Y_p)^{\mathrm{T}}$$ has the multivariate normal distribution with mean vector $$\mu_2$$ and covariance matrix $$\Sigma_2 > 0$$, i.e. $$Y \sim \mathcal{N}_p(\mu_2, \Sigma_2)$$, and has the density $$f_2(y)$$, $$y \in \mathbb{R}^p$$.

I encountered $$D_p = \int_{\mathbb{R}^p}\int_{\mathbb{R}^p} (x - y)^{\mathrm{T}}(x - y) \, f_1(x)f_2(y) \, dy dx.$$

Mathematica 12 showed me: \begin{align*} D_1 & = \int_{-\infty}^\infty\int_{-\infty}^\infty (x - y)^2 f_1(x) f_2(y) \, dydx \\ & = (\mu_1 - \mu_2)^2 + (\sigma_1^2 + \sigma_2^2), \end{align*} where $$f_1(x) = \frac{1}{\sqrt{2\pi}\,\sigma_1} \exp\biggl[-\frac{1}{2}\left(\frac{x - \mu_1}{\sigma_1}\right)^2 \, \biggr], \quad f_2(y) = \frac{1}{\sqrt{2\pi}\,\sigma_2} \exp\biggl[-\frac{1}{2}\left(\frac{y - \mu_2}{\sigma_2}\right)^2 \, \biggr],\\$$

\begin{align*} D_2 & = \int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty [(x_1 - x_2)^2 + (y_1 - y_2)^2 ] f_1(x_1, y_1) f_2(x_2, y_2) \, dy_2dy_1 dx_2dx_1 \\ & = (\mu_{11} - \mu_{21})^2 + (\mu_{12} - \mu_{22})^2 + (\sigma_{11}^2 + \sigma_{12}^2) + (\sigma_{21}^2 + \sigma_{22}^2)\\ & = (\mu_1 - \mu_2)^{\mathrm{T}}(\mu_1 - \mu_2) + \mathrm{tr}(\Sigma_1) + \mathrm{tr}(\Sigma_2) \end{align*} where $$\mu_1 = (\mu_{11}, \mu_{12})^{\mathrm{T}}, \quad \mu_2 = (\mu_{21}, \mu_{22})^{\mathrm{T}},$$ $$\Sigma_1 = \begin{pmatrix} \sigma_{11}^2 & \rho_1\sigma_{11}\sigma_{12} \\ \rho_1\sigma_{11}\sigma_{12} & \sigma_{12}^2\\ \end{pmatrix}, \quad \Sigma_2 = \begin{pmatrix} \sigma_{21}^2 & \rho_2\sigma_{21}\sigma_{22} \\ \rho_2\sigma_{21}\sigma_{22} & \sigma_{22}^2\\ \end{pmatrix},$$ $$|\rho_1| \leq 1, \quad |\rho_2| \leq 1.$$

I am guessing $$$$D_p = (\mu_1 - \mu_2)^{\mathrm{T}}(\mu_1 - \mu_2) + \mathrm{tr}(\Sigma_1) + \mathrm{tr}(\Sigma_2). \tag1\label1$$$$

How can I prove equation \eqref{1}?

• Welcome to Math.SE! I have tried to improve the readability of your question by improving the $\rm \LaTeX$ code. It is possible that I unintentionally changed the meaning of your question. Please proofread the question to ensure this has not happened. – GNUSupporter 8964民主女神 地下教會 Dec 18 '20 at 10:14
• Thank you for your proper correction to improve readability. – Eulers Dec 18 '20 at 10:20

Suppose that $$X$$ and $$Y$$ are independent. Put $$c = \mu_1 - \mu_2$$, $$A = X-EX$$, $$B = Y-EY$$. Hence $$D_p = E(X-Y)^{T}(X-Y) =E(X-EX - (Y- EY) + c)^{T}(X-EX - (Y-EY) +c) = E(A-B+c)^{T}(A-B+c).$$
We have $$A \sim N_p(0, \Sigma_1)$$ and $$B \sim N_p(0, \Sigma_2)$$. Hence $$EA=EB=0$$, $$EA^{T}=EB^{T}=0$$, $$DA = \Sigma_1$$, $$EA^{T}A = tr(\Sigma_1), EB^{T}B = tr(\Sigma_2)$$. Thus
$$D_p = E(A-B)^{T}(A-B) + E(A-B)^{T}c + Ec^{T}(A-B)+Ec^{T}c$$ and
$$D_p - c^{T}c = E(A-B)^{T}(A-B) = EA^{T}A - EB^{T}A - EA^{T}B + E B^{T}B = EA^{T}A +E B^{T}B$$ because $$EB^{T}A = EB^{T} EA = 0$$ and $$EA^{T}B = EA^{T}EB = 0$$ by independence.
$$D_p = c^{T}c + tr(\Sigma_1) + tr(\Sigma_2) = (\mu_1 - \mu_2)^{T}(\mu_1 - \mu_2) + tr(\Sigma_1) + tr(\Sigma_2)$$ q.e.d.