# 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.