Recently I wanted to generate random multivariate normal numbers. I was a bit surprised--even, I'll admit, mildly irritated--to discover that there was not, among the many billion pages searched by Google, a simple, free, easily-used computer program for this purpose. This is despite (1) the computational simplicity of the task, and (2) the commonness of the need.
There are several algorithms in source code for this on the web. However there's a difference between having a source code algorithm that can *theoretically* be simply complied and run on an ordinary PC, and having an actual, working program. Producing the latter from the former too often requires specialized knowledge of computers, programming languages, and, in some cases numerical analysis.
As this is something of potential interest to scientists and researchers in many fields--biology, earth sciences, physical sciences, engineering, statistics, and psychology, to name a few--and since I had to write a program for this anyway--it seemed a good idea to make the program available to others. So...here it is.
To download the program, click this link: mvn.zip
In any case, running this program requires no more than clicking on the program icon. A Command Prompt window will open automatically in which the program will run. (Note: This is the easiest, but not necessarily the best way; better is to open the Command Prompt window first, navigate to the folder with the MVN program, then type the command, MVN, and press the enter key. This will ensure the window stays open in the event of any errors.)
To produce random multivariate normal numbers, MVN first generates random univariate normal numbers using the TOMS Algorithm 712 by JL Leva. The full reference is: Leva JL. Algorithm 712. A normal random number generator. ACM Transactions on Mathematical Software (TOMS), v.18 n.4, pp. 454-455, Dec. 1992. (Note this reference applies only to the generation of random univariate normal numbers, which is a step in the generation of random multivariate normal numbers.)
The algorithm above uses the ratio of uniforms method of AJ Kinderman and JF Monahan augmented with quadratic bounding curves (citation needed).
Uniform random numbers, used by this algorithm, are supplied by the default random number function of the Fortran 90 compiler (Absoft Pro Fortran 90, v. 7).
The random multivariate normal numbers are produced by pre-multiplying a vector of random univariate normal numbers by the Cholesky decomposition of the correlation matrix according to the formula:
Y = L X
where
Here the Cholesky decomposition is stored in the lower triangle and main diagonal of a square matrix; elements in the upper triangle of the matrix are 0.
Standard deviations are then multiplied and/or means added per the user specifications.
Because the aim is to provide something convenient and helpful, I am especially interested in learning of any program bugs promptly. In short, if you discover a bug, I'll try to drop everything else and correct it immediately. I just ask that you first consult the Troubleshooting section of the Readme.txt document that accompanies the program.
Last updated: 21 August 2006 (added counter)