Incorrect calculation of logpdf for multivariate normal distribution about distributions.jl HOT 10 CLOSED

I made a quick correction at 19d958e

from distributions.jl.

Thanks.

from distributions.jl.

Yes, we badly need tests. I would say that everything that's not being brought in from R should be treated as alpha or even pre-alpha software.

from distributions.jl.

@lindahua I don't think the old version was wrong. Please note that the mean is standardized with respect to the cholesky factor of the covariance matrix and not the covariance itself. Also, in the new version the result of your own example is

julia> logpdf(g,x)
-21.291568715206875

I have just tried to reproduce the MATLAB result and first i couldn't, but then I tried

julia> log(1/(sqrt(2pi)^3*sqrt(det(c)))*exp(-0.5*x'*(c\x)))
1-element Float64 Array:
 -15.7554

The problem is that the pdf is tiny, 1.43721e-7, and hence the calculation should be done in logs. MATLAB and Mathematica don't do that and get results with large errors.

from distributions.jl.

Okay, please ignore last part of my message. I should have thought a little more about it. The new and my old versions are still wrong but my idea that MATLAB and Mathematica had it wrong was not that smart. I'll think it through after lunch.

from distributions.jl.

The problem was the transpose of the cholesky factor. I have pushed a fix 19d958e.

from distributions.jl.

@andreasnoackjensen Oh, it is my fault I did not test the correction myself. The problem was that I did not notice it was d.covchol.LR instead of d.covchol. The correct version of the quadratic part should be
(x - d.mean) * (d.covchol \ (x - d.mean)). Instead of (x - d.mean) * (d.covchol.LR \ (x - d.mean))

from distributions.jl.

The pdf is not very large (at the level of 1.4e-7, but this is not enough to produce such big difference. Actually, a correct implementation will produce very accurate result as long as the pdf is not over/underflow.

from distributions.jl.

You are right about the 1.4e-7. I should have waited half an hour and looked at it again before posting. However, hopefully we now have a correct version. You are also completely right about the tests. It is so easy to get these things wrong.

from distributions.jl.

@andreasnoackjensen I made a revision at 819db1f

It is now correct (I added a test to ensure this) and much faster (see my gist ).

from distributions.jl.