Look into switching from `MASS::mvrnorm()` to `mvtnorm::rmvnorm()`

[john-harrold]

Refer to the discussion here:

https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/

[mattfidler]

As a note, rxode2 also uses the same methodology as mvtnorm I believe, ported and parallelized to use sitmo engine.

It is also parallelized based on the id of the simulation since sitmo can use serialized numbers.

@djnavarro didn't even cover that some implementations of random number generators may be different on different platforms since they use different compilers, nor the oddness of random numbers with low seeds in certain random number generators, or even the birthday problem when applied to seeds.

There are many pitfalls you may randomly fall into 😄

[djnavarro]

Agreed, and given that the original context was a PMX problem, it's good to know that rxode2 uses this approach. Something that John Lamping (I think on here as @johnlamping) later pointed out to me (and will eventually appear in the blog post once I recover from covid brain fog) is that the mvtnorm default is a matrix square root method based on eigendecomposition, which is why it ends up robust to those sign flips. But yes, there's so many other ways we can run into difficulties (I didn't dare talk about compilers in that post but yeah... compilers)

[mattfidler]

In particular clang ang gcc use different ways to generate random numbers of different distributions. While neither approach is wrong per se, it means that code compiled with GCC's normal generator and clang's normal generator would be different even if the underlying uniform pseudo random-number generator is the same. It is more than simply differences in how floating point numbers change things.

The solution (in reprospect) is to use random number generators from an external library like boost to keep the numbers platform/compiler independent.

The finding is important because R random numbers cannot be parallelized on the C/C++ level which is how rxode2 and nlmixr2 achieve parallelization. It is also how other packages have parallel random number generation and can lead to much larger platform inconsistencies.

That also means that the results from rxode2 would likely be different from the mvtnorm package since the random number generators and the way to get normal variables in general are different than the R ones. Still, the method once these numbers are generated are the same.

Unfortunately, like you, most of my knowledge comes from problems observed and reported. In my case it is in developing rxode2 and nlmixr2.

In this particular case, I used mvtnorm simply because I could follow the code well and test against specific parts of the code to make sure I was doing the right thing.

[john-harrold]

I think I found the function in rxode2 that does this, and I couldn't really tell exactly what it was doing by glancing at it. So I was going to bring this up in our next development meeting :).

[mattfidler]

I was wrong, I use the methodology here

https://github.com/lbelzile/TruncatedNormal

Instead, though the cholesky decomposition is the way that it is handled 😄