The pdf of multivariate normal distribution with high correlation values

Question

Does the pdf of an mvn variable even exist when there is high correlation?

I want to use an algorithm (actually it is the cross-entropy method for estimating a rare-event probability) that needs the pdf value of the mvn distribution, that is mvnpdf(x,mu,Sigma).

However, my Sigma is close to singular, meaning there is a high correlation between the variables in the vector, so it is of course difficult/impossible to find the inverse of Sigma. Is there any way to overcome this problem?

Is it not true that for instance a vector [a,b,c,d] with covariance matrix [1,1,0,0;1,1,0,0;0,0,1,1;0,0,1,1] (singular!) will behave in the exact same way as the vector [a,c] with covariance matrix [1,0;0,1] (now non-singular!), while ignoring b and c? Does this mean that I can approximate the pdf value of [a,b,c,d] by the pdf of [a,c]?

Sorry if this is answered before, I've really tried to search for it.

Thank you.

EDIT: The thing is that I want to simulate samples from a multivariate normal distribution, $x \sim N(\mu,\sigma)$, in order to find the probability that ${r(x) < 1}$, where $r(x)$ returns a positive real number based on the vector $x$. So I simulate $x_i \sim N(\mu,\Sigma), i = 1 ,..., M$, and I estimate the probability by the Monte Carlo estimate $p = (1/M)\sum_i I(r(x_i)<1)$. No problems so far. However, since ${r(x_i)<1}$ is a very rare event, I want to try importance sampling instead, simulating from the distr $N(\mu_2,\Sigma_2)$ so I need the pdf value to estimate $p = (1/M)\sum{I(r(x_i)<1) pdf_1(x_i)/pdf_2(x_i)}$.

Answer 1

In your example of $[a,b,c,d]$, $a$ and $b$ are the same random variable for all intents and purposes since they take on the same value with probability $1$. Similarly, $c$ and $d$ are the same random variable. The question is what you want to do with the multivariate distribution. If you want to calculate $P\{a \in A, b \in B, c\in C, d \in D\}$ where $A$, $B$, $C$, $D$ are (measurable) sets of real numbers, you simplify it to $P\{a \in (A\cap B), c \in (C\cap D)\}$ and calculate it from the bivariate density of $a$ and $c$ which in this instance is that of two independent standard normal random variables. In other words, you don't blindly apply the formulas for the multivariate normal density of four normal random variables; you think a bit and re-state the problem to be solved in a way that it can be solved.

Answer 2

As long as the covariance matrix is symmetric and satisfies all other properties to be acovariance matrix the multivariate normal distribution exists. Now the covariance matrix can even be singular. Then of course the density cannot be written with the exponent $-(x-μ)'\Sigma^{-1}(x-μ)/2$ as it "lives" (has all its probability mass) on a lower dimensional space of rank r

Answer 3

As @Michael Chernick suggests, the Gaussian exists, but its pdf is undefined. And as @Dilip Sarwate suggests, you can find a work around.

In your case, if the correlation is high but not 1 you can still generate a Gaussian vector with a definite pdf, so let's assume that your covariance matrix is strictly singular. Because it is a symmetric real matrix, you can still diagonalize it, but at least one of the eigenvalues will be 0.

If you are using R and mat is the singular covariance matrix, you can do the following. Here I assume that you have only 1 null eigen value.

# n is the dimension of your original Gaussian vector.
n <- nrow(mat)
eigen. <- eigen(mat)
# lambas are eigen values, the last one is 0.
lambdas <- eigen.[["values"]]
# P is the matrix of eigen vectors for linear transformation.
P <- eigen.[["vectors"]]

We can now sample IID (0,1) Gaussian variables and apply the linear transform P to regenerate an n-Gaussian vector vectors. In the mean time, we keep the value of the pdf for the purposes of sampling importance.

# Simulate an (n-1)-Gaussian vector Z. and get value of pdf.
Z. <- rnorm(n-1)
pdf_Z. <- prod(dnorm(Z.))
# Regenerate an n-Gaussian vector by appending 0, scaling...
Z <- c(Z., 0) * sqrt(lambdas)
# ... and applying matrix of eigen vectors and its inverse.
X <- P %*% diag(Z, ncol=n) %*% solve(P)

Now X has the property you want, and you can use pdf_Z. in your weight.

Note: I tried to make the code above didactic to understand the point. I be it is not efficient and even if you use R you would have to recode that bit.