R package for simulating draws from a truncated normal random variable

I recently implemented an algorithm where I have to sample repeatedly from a (univariate) truncated normal random variable. This is a well-explored problem. Different accept-reject algorithms have been proposed by John Geweke, Christian Robert and Nicolas Chopin.

However, I was unable to find a state-of-the-art implementation in R - my go-to statistical software. The best implementation that I could track down was the rtruncnorm function from the truncnorm package. I was curious to see whether I could improve on the performance of rtruncnorm.

One of the fastest algorithms for drawing from the truncated normal distribution is the algorithm by Nicolas Chopin. This algorithm has already been implemented in C++ by Guillaume Dollé and Vincent Mazet and all I had to do was to adapt their code to the Rcpp environment.

The result of this effort is the new rtnorm package. It can be installed easily via the install_github command from the devtools package. Detailed installation instructions can be found on the github page of the package. The package supplies a single function - rtnorm - that draws from the truncated normal distribution.

For the most part, the rtnorm conforms with the usual expectations for R functions. In particular, the bounds of the truncation interval are allowed to be positive or negative infinity.

There are two dimensions along which the rtnorm function deviates from the expectations for R functions:

  • It draws only once from the truncated normal distribution, i.e., it is not vectorized. In most applications this is not restrictive. For example, in applications in Gibbs sampling only one draw at a time is needed.
  • It does not fill in default values for the mean and the standard deviation of the normal distribution (before truncation). The only way I could make default arguments work was through an additional wrapper function in pure R. This generated quite a lot of overhead that simply did not seem worth it.

Here are some benchmarks for 10e4 random draws from a standard normal distribution truncated to the the interval [1, 2]

         test replications elapsed relative user.self sys.self user.child sys.child
2      rtnorm            1    0.86    1.000      0.86        0         NA        NA
1   truncnorm            1    6.59    7.663      6.60        0         NA        NA
3 msm::rtnorm            1   13.87   16.128     13.87        0         NA        NA

and the interval [5, 7]

         test replications elapsed relative user.self sys.self user.child sys.child
2      rtnorm            1    0.92    1.000      0.93        0         NA        NA
1   truncnorm            1    6.43    6.989      6.44        0         NA        NA
3 msm::rtnorm            1   13.23   14.380     13.24        0         NA        NA

These benchmarks are measured on my Windows computer. On my Mac the truncnorm package performs considerably better, though still not anywhere close to as fast as the new rtnorm package.

The results of the second benchmark are actually quite puzzling, since for truncation to the tails the algorithms used in rtnorm and truncnorm are quite similar. My conjecture is that rtnorm is better at avoiding unnecessary overhead from communication between R and C/C++. However, I am not sufficiently well versed in how R communicates with dynamic C/C++ libraries to investigate this conjecture more closely.

The package is already fully functional but there is still some work left do. I have informally tested the package against the truncnorm package but I have not written a full suite of automated unit tests. I plan to add this in the future.

Putting together the package was a great learning experience and reminded me to use Rcpp more often in my work.

Categories: ,

Updated: