A Problem with mvnrnd and mvnpdf function

when I was using mvnrnd and mvnpdff function, the code is like this:

1. Mu=[3100,1600,2.28];
2. Cx=[61723.4468937876,43206.4128256513,7.89178356713429;
3. 43206.4128256513,30244.4889779559,5.52424849699400;
4. 7.89178356713429,5.52424849699400,0.00102240480961];
5. f1 = mvnrnd(Mu,Cx,1);
6. p1 = mvnpdf(f1,Mu,Cx);

When I run this code, I got the result like this: f1 = [3332,1762,2.322],and p1 = 55.035. So what I want to know is the result of p1 means? and the probability of f1.

Answer 1

[Roger Stafford]

You have created only a single 1-by-3 random sample, f1, with 'mvnrnd' and p1 is simply the 3D probability density, not probability, at that particular point for your given Mu and Cx values. Normal stochastic processes are distributed continuously over the entire 3D space and consequently the probability of each individual point is zero. The fact that p1 is greater than one is not a surprise. Probability density is a measure of probability per unit 3D volume and this can easily exceed one. It is the triple integral of the probability density over any given volume, which is to say the probability of a random triple falling in that volume, that must not exceed one.

[John M]

Thank you so much! I thought the function was wrong when I got that answer! you said "probability of each individual point is zero". but when we use randn to generate like 10000 numbers. it's clearly that the probability to get a 0.5 is much greater than 5. how we use math or number to express this kind of relationship? "Probability density is a measure of probability per unit 3D volume". So can we say if Probability density is higher and the Probability is higher too?

[Roger Stafford]

In the ideal mathematical world where there are infinitely many points on the real line, the probability of each individual point with any normal (gaussian) distribution must be said to be zero. However, the probability of belonging to some interval of non-zero length will be greater than zero. This sounds paradoxical but that is the way things are when dealing with infinities.

In the practical world of computers, there are only a finite number of possible values for a normal distribution, so each individual value does admittedly have a non-zero probability, but for double precision values that value is extremely small, so small that for all practical purposes it can be considered to be zero. That is because there are so many possible double precision numbers. For example, using the 'randn' function, try repeating it until you get an exact zero, which is its most likely value. You will have to wait a very long time to see this happen. I started it on my computer just now and had time to eat breakfast and it still hadn't found zero after many millions of trials.

Thus it is best to speak of the probabilities of intervals or of sets of values which have a non-zero measure. For example the probability that randn gives an answer between 0 and 1 is 0.34134474606854. On the other hand there is real meaning to the probability density at a single value. The probability density of randn at x is the limit as d approaches zero of the ratio between the probability of lying between x-d and x+d and then that probability divided by 2*d. It is in fact 1/sqrt(2*pi)*exp(-x^2/2). For normal distributions with a small variance that density can easily be greater than zero, as you have found. The determinant of your matrix Cx was extremely small, 3.8561e-12, almost singular, which accounts for the high density you saw.