A fix to the mvrnorm function

Recently when our Bayesian Data Analysis class was doing Gibbs sampling in R, some students noticed that they get missing values when sampling from multivariate normal distribution using the function mvrnorm in R (in the package MASS). The function seems to fail in some cases. This is due to some weirdness in the basic R function eigen that computes eigenvalues, so mvrnorm itself is not to blame. However mvrnorm might be the function you use more often so I wrote a simple patched version of mvrnorm which should fix this problem. It will also alert the user if the patch fails.

Details and more are here.

3 thoughts on “A fix to the mvrnorm function

  1. I can't reproduce the problem with R 2.0.0 on

    a Intel Xeon system. However, I do recall that

    one of my PhD students had a similar problem with

    eigen, that could be solved by using the option

    EISPACK = TRUE to force use of the older eigenvalue

    routines.

  2. Dear Radford,

    thank you for the hint. EISPACK=TRUE does indeed fix the problem (at least for this matrix) on my computer. I changed my patch to try symmetric=TRUE, EISPACK=TRUE if eigen fails.

    I also announced this in the R-help mailing list.

  3. I learned that this problem is taken care of in the latest versions of mvrnorm by setting EISPACK=TRUE. Just upgrade R to the latest version of R and you'll get the latest version of mvrnorm along with it. Our mvrnorm was an old one.

    Anyway for anyone computing eigenvectors, knowing that eigen() may fail is an useful thing to know. I've left the instructions on my web page.

    The eigen() function doesn't seem to fail always, and it doesn't fail always on all computers, but that's the most annoying case of software: ideally software should either fail always or work always :-)

Comments are closed.