Find the marginal distributions (PDFs) of a multivariate normal distribution

Question

Let $\mathbf{x}\in\Bbb{R}^n$ be a multi-variate normal vector with mean $\bar{\mathbf{x}}\in\Bbb{R}^n$ and covariance matrix $\Sigma\in\Bbb{S}_{++}^n$, where $\Bbb{S}_{++}^n$ denotes the space of all $n\times n$ symetric posotive definite matrices with real entries. Also, let $f$ denote the probability density function of $\mathbf{x}$, i.e., $$ f(\mathbf{x}) = \frac{1}{(2\pi)^{\frac{n}{2}}\vert\Sigma\vert^{\frac{1}{2}}} \exp \left( -\frac{1}{2} (\mathbf{x}-\bar{\mathbf{x}})^\top \Sigma^{-1} (\mathbf{x}-\bar{\mathbf{x}}) \right). $$ In general we suppose that $\Sigma$ is a non-diagonal matrix. Thus, the elements of $\mathbf{x}=(x_1,\ldots,x_j,\ldots,x_n)^\top$, i.e., the variables $x_j$, $j=1,\ldots,n$ are dependent normal variables.

I am interested in the computation of the marginal densities of $x_j$, $j=1,\ldots,n$. In the corresponding Wikipedia's article, it states that:

To obtain the marginal distribution over a subset of multivariate normal random variables, one only needs to drop the irrelevant variables (the variables that one wants to marginalize out) from the mean vector and the covariance matrix.

It also gives the following example:

Example: Let $\mathbf{x}=(x_1, x_2, x_3)^\top$ be multivariate normal random variables with mean vector $\bar{\mathbf{x}}=(\bar{x}_1,\bar{x}_2, \bar{x}_3)^\top$ and covariance matrix $\Sigma$ (standard parametrization for multivariate normal distributions). Then the joint distribution of $\mathbf{x}^\prime=(x_1, x_3)^\top$ is multivariate normal with mean vector $\bar{\mathbf{x}}^\prime=(\bar{x}_1, \bar{x}_3)^\top$ and covariance matrix $$ \Sigma^\prime = \begin{bmatrix} \Sigma_{11} & \Sigma_{13} \\ \Sigma_{31} & \Sigma_{33} \end{bmatrix}. $$

Does this mean that, if we want to obtain the marginal probability density function of each of $x_j$, for all $j=1,\ldots,n$, we just use the diagonal entries of $\Sigma$ as the variances of the respective random variables? What bugs me is that this would be the same as if $\Sigma$ was diagonal (and thus $x_j$'s independent). Does the solution lie in the parentheses of the above example (i.e., standard parametrization for multivariate normal distributions)?

Answer 1

Allow me to start with a simple example so I can make my point more evident. Assume you have three random variables $x, y$ and $z$, where $x$ and $y$ are dependent while $x$ and $z$ are independent. Then create two random vectors $\mathbf{x_1}=(x,y) ^\top$ and $\mathbf{x_2}=(x,z) ^\top$. If you calculate the covariance matrix of $\mathbf{x_1}$ you will get a non-diagonal covariance matrix where the off diagonal elements will indicate the relationship of $x$ and $y$. However the diagonal elements will provide information about the dispersion of $x$ and $y$ solely. Now in the case of $\mathbf{x_2}$, the matrix will be diagonal and the diagonal elements will also provide the same information as before, that is how $x$ and $y$ are spread as scalars. The fact that $x$ and $y$ are somehow connected doesn't affect the behavior of $x$ if you study it ignoring $y$.

Now in a more formal manner, each element $\Sigma_{ij}$ of the covariance matrix $\Sigma$, is calculated by $\Sigma_{ij} = \mathbb{E}[(x_i-\mu_i)(x_j-\mu_j)]$. Hence the elements of the main diagonal are calculated as $\Sigma_{ii} = \mathbb{E}[(x_i-\mu_i)(x_i-\mu_i)]$ wich is also the definition of variance.

Please note that I didn't assume that the random variables were normal.