Find the marginal distributions (PDFs) of a multivariate normal distribution 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)?
 A: The easiest way to convince yourself that the statement from Wikipedia you highlighted above is true is by use of the characteristic function. As we know it is defined as follows:
\begin{equation}
\kappa_\vec{X}(\vec{k}) := E\left[ e^{\imath \vec{k} \cdot \vec{X}} \right]
\end{equation}
Now in case $\vec{X}$ is multivariate normal with means $\vec{\mu}$ and covariances ${\bf \Sigma}$,i.e. $\vec{X} \sim N(\vec{\mu},{\bf \Sigma})$ then from https://en.wikipedia.org/wiki/Multivariate_normal_distribution we have:
\begin{equation}
\kappa_\vec{X}(\vec{k}) = e^{\imath \vec{\mu}^T \cdot \vec{k} - \frac{1}{2} \vec{k}^T \cdot {\bf \Sigma} \cdot \vec{k}}
\end{equation}
Now to marginalize in $\vec{k}$-space is easy. All we need to do is to set the respective components, i.e. those corresponding to variables we want to marginalize over in real space, of the $\vec{k}$-vector to zero. therefore we take ${\vec{k}_1} := \left( 0,\cdots 0, k_j, 0, \cdots, 0 \right)$ where the non-zero element sits at the jth position. Then clearly we have:
\begin{equation}
\kappa_\vec{X}(\vec{k}_1) = e^{\imath \mu_j k_j - \frac{1}{2} k_j {\bf \Sigma}_{j,j} k_j}
\end{equation}
Now in order to get the marginal distribution all we need to do is to invert the Fourier transform. We have:
\begin{eqnarray}
\rho^{(marg)}_j(x_j) &:=& \frac{1}{2\pi} \int\limits_{-\infty}^{\infty} e^{-\imath k_j x_j} \kappa_\vec{X}(\vec{k}_1) d k_j \\
&=& \frac{1}{2\pi} \int\limits_{\mathbb R} e^{-\imath k_j (x_j-\mu_j) - \frac{1}{2} {\bf \Sigma}_{j,j} k_j^2} dk_j \\
&=&e^{-\frac{1}{2} \frac{(x_j-\mu_j)^2}{{\bf \Sigma}_{j,j}} } \cdot 
\frac{1}{2\pi} \cdot \int\limits_{\mathbb R} e^{-\frac{1}{2} {\bf \Sigma}_{j,j} k_j^2 }dk_j \\
&=&\frac{1}{\sqrt{2 \pi {\bf \Sigma}_{j,j}}} e^{-\frac{1}{2} \frac{(x_j-\mu_j)^2}{{\bf \Sigma}_{j,j}} }
\end{eqnarray}
as expected.
A: 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. 
