In [1] is noted, that a covariance matrix is "positive- semi definite and symmetric". I wonder if it is possible to a multivariate normal distribution with a covariance matrix that is only positive semi-definite but not positive definite?
If yes, how can the density be evaluated, since it involves the inverse and the 1/determinant of the covariance matrix.
As the commenters have already mentioned, there isn't a probability density function in the case where the covariance matrix is singular. Rather, you have a distribution that lives on a lower dimensional subspace of $R^n$.
For example, suppose $X_{1} \sim N(0,1)$, and $X_{2}=-X_{1}$. The covariance of $X_{1}$ and $X_{2}$ is -1, and the variances of $X_{1}$ and $X_{2}$ are both +1. This covariance matrix is positive semidefinite, but singular.
Since $x_{2}=-x_{1}$, the "probability density" must be 0 everywhere off of this line. However, you still need the probability distribution to integrate out to 1. No function from $R^2$ to $R$ can do this, so there isn't actually a probability density. Rather, you have delta function like distribution that lives on the line $ x_{2}=-x_{1}$.
If you haven't studied enough analysis to work with such distributions, then be very careful about this. Even if you have studied enough analysis to understand this, beware that doing anything numerically with such a covariance matrix can be very difficult.
Happens all the time in ANOVA and regression problems. Here's a really simple one: Suppose $$ Y_i = \alpha_0 x_i + \alpha_1 + \varepsilon_i $$ where $\varepsilon_i \sim N(0,\sigma^2)$ and the $\varepsilon$s are independent.
Now let $\hat{\alpha}_0$ and $\hat{\alpha}_1$ be the least-squares estimates. Let $\hat{Y}_i = \hat{\alpha}_0 x_i + \hat\alpha_1$ be the fitted values. Let $\hat\varepsilon_i= Y_i-\hat Y_i$ be the residuals.
Now notice the difference between the errors $\varepsilon_i$ and the residuals $\hat\varepsilon_i$. The former are independent, so the probability that their sum is $0$ is $0$. The latter necessarily sum to $0$ and also satisfy the identity $\sum_i x_i\hat\varepsilon_i=0$. So they're not independent. Thus there are $n-2$ degrees of freedom for error. (They're also not identically distributed, by the way. That's one thing that complicates the studentization of residuals.)
The vector of errors is distributed as $N(0,\sigma^2 I_n)$ where $I_n$ is the $n\times n$ identity matrix and the $0$ is an $n\times 1$ column vector. The vector of residuals has the same expected-value vector, but its variance is a singular matrix of rank $n-2$. It's $\sigma^2$ times an orthogonal projection matrix onto the orthogonal complement of the space spanned by $(1,\ldots,1)^T$ and $(x_1,\ldots,x_n)^T$.
Multivariate normal distributions with singular variances arise in ways like this incessantly in statistics.
More simply, suppose $A$ is any non-negative-definite symmetric real matrix. The finite-dimensional spectral theorem says $A$ has a non-negative-definite square root $A^{1/2}$. Let $Z = (Z_1,\ldots,Z_n)^T$ be a vector of i.i.d. standard normals. Then the variance (or "covariance matrix", if you like) of $A^{1/2}Z$ is $A$.
So every non-negative-definite symmetric real matrix is realizable as the variance of some vector-valued random variable.
BTW, shall we be clear about the definition? A random vector has a multivariate normal distribution if its dot-product with every constant vector has a univariate normal distribution. (I mention this in part because the question about densities made me wonder if some people think densities are essential to the concept. They're not part of any of the standard definitions.)
