[R] generating multivariate normal random numbers in C

Prof Brian Ripley ripley at stats.ox.ac.uk
Tue Jun 26 18:46:01 CEST 2001

On Tue, 26 Jun 2001, Mathieu Ros wrote:

> >>>>> "TL" == Thomas Lumley <tlumley at u.washington.edu> disait:
>     TL> On Tue, 26 Jun 2001, Mathieu Ros wrote:
>     >> hello all, I'm looking for a C/C++ library (or algorithm) to
>     >> compute draws from a multivariate normal distribution in order
>     >> to do block sampling.  (I'm plugging my results in an R
>     >> function so this is not really off topic ;)
>     TL> In that case, why not use mvrnorm() in the MASS package? It's
>     TL> not in C, but it doesn't need to be.
> Well, as I just replied to prof. Rossini, I'm doing about 50000
> iterations (up to 100000) with at least 88 parameters to update at each turn
> (at least 2x4 + 2x40 if I do block sampling).
> I don't know much about R/C interface but I believe it will slow down
> my program if I call an R function to do the job.
> My elementwise sampling algorithm (C code) takes between 3 and 4
> hours to run (with intermediate data sets).

Well, it's an elementary computation once you have a random normal
generator and a numerical linear algebra suite.  For 88 dimensions (if I
understand you aright) the issue is going to be the matrix decomposition.
Is this the same multivariate normal distribution at each draw?  If so
just write a binary file using writeBin: it will not slow you down.  If it
is different each time the linear algebra computations will be slow enough
to swamp the communications overhead.

Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch

More information about the R-help mailing list