Generate random numbers following a distribution within an interval

Question

I need to generate random numbers following Normal distribution within the interval $(a,b)$. (I am working in R.)

I know the function rnorm(n,mean,sd) will generate random numbers following normal distribution,but how to set the interval limits within that? Is there any particular R functions available for that?

Answer 1

It sounds like you want to simulate from a truncated distribution, and in your specific example, a truncated normal.

There are a variety of methods for doing so, some simple, some relatively efficient.

I'll illustrate some approaches on your normal example.

1. Here's one very simple method for generating one at a time (in some kind of pseudocode):

   $\tt{repeat}$ generate $x_i$ from N(mean,sd) $\tt{until}$ lower $\leq x_i\leq$ upper

   

   If most of the distribution is within the bounds, this is pretty reasonable but it can get quite slow if you nearly always generate outside the limits.

   In R you could avoid the one-at-a-time loop by computing the area within the bounds and generate enough values that you could be almost certain that after throwing out the values outside the bounds you still had as many values as needed.

2. You could use accept-reject with some suitable majorizing function over the interval (in some cases uniform will be good enough). If the limits were reasonably narrow relative to the s.d. but you weren't far into the tail, a uniform majorizing would work okay with the normal, for example.

   

3. If you have a reasonably efficient cdf and inverse cdf (such as pnorm and qnorm for the normal distribution in R) you can use the inverse-cdf method described in the first paragraph of the simulating section of the Wikipedia page on the truncated normal. [In effect this is the same as taking a truncated uniform (truncated at the required quantiles, which actually requires no rejections at all, since that's just another uniform) and apply the inverse normal cdf to that. Note that this can fail if you're far into the tail]

   

4. There are other approaches; the same Wikipedia page mentions adapting the ziggurat method, that should work for a variety of distributions.

The same Wikipedia link mentions two specific packages (both on CRAN) with functions for generating truncated normals:

The MSM package in R has a function, rtnorm, that calculates draws from a truncated normal. The truncnorm package in R also has functions to draw from a truncated normal.

---

Looking around, a lot of this is covered in answers on other questions (but not exactly duplicates since this question is more general than just the truncated normal) ... see additional discussion in

a. This answer

b. Xi'an's answer here, which has a link to his arXiv paper (along with some other worthwhile responses).

Answer 2

The quick-and-dirty approach is to use the 68-95-99.7 rule.

In a normal distribution, 99.7% of values fall within 3 standard deviations of the mean. So, if you set your mean to the middle of your desired minimum value and maximum value, and set your standard deviation to 1/3 of your mean, you get (mostly) values that fall within the desired interval. Then you can just clean up the rest.

minVal <- 0
maxVal <- 100
mn <- (maxVal - minVal)/2
# Generate numbers (mostly) from min to max
x <- rnorm(count, mean = mn, sd = mn/3)
# Do something about the out-of-bounds generated values
x <- pmax(minVal, x)
x <- pmin(maxVal, x)

I recently faced this same problem, trying to generate random student grades for test data. In the code above, I've used pmax and pmin to replace out-of-bounds values with the min or max in-bounds value. This works for my purpose, because I'm generating fairly small amounts of data, but for larger amounts it will give you noticeable bumps at the min and max values. So depending on your purposes it may be better to discard those values, replace them with NAs, or "re-roll" them until they're in-bounds.

Answer 3

None of the answers here give an efficient method of generating truncated normal variables that does not involve rejection of arbitrarily large numbers of generated values. If you want to generate values from a truncated normal distribution, with specified lower and upper bounds $a<b$, this can be done ---without rejection--- by generating uniform quantiles over the quantile range allowed by the truncation, and using inverse transformation sampling to get corresponding normal values.

Let $\Phi$ denote the CDF of the standard normal distribution. We want to generate $X_1,...,X_N$ from a truncated normal distribution (with mean parameter $\mu$ and variance parameter $\sigma^2$)$^\dagger$ with lower and upper truncation bounds $a<b$. This can be done as follows:

$$X_i = \mu + \sigma \cdot \Phi^{-1}(U_i) \quad \quad \quad U_1,...,U_N \sim \text{IID U} \Big[ \Phi \Big( \frac{a-\mu}{\sigma} \Big), \Phi \Big( \frac{b-\mu}{\sigma} \Big) \Big].$$

There is no inbuilt function for generated values from the truncated distribution, but it is trivial to program this method using the ordinary functions for generating random variables. Here is a simple R function rtruncnorm that implements this method in a few lines of code.

rtruncnorm <- function(N, mean = 0, sd = 1, a = -Inf, b = Inf) {
  if (a > b) stop('Error: Truncation range is empty');
  U <- runif(N, pnorm(a, mean, sd), pnorm(b, mean, sd));
  qnorm(U, mean, sd); }

This is a vectorised function that will generate N IID random variables from the truncated normal distribution. It would be easy to program functions for other truncated distributions via the same method. It would also not be too difficult to program associated density and quantile functions for the truncated distribution.

---

$^\dagger$ Note that the truncation alters the mean and variance of the distribution, so $\mu$ and $\sigma^2$ are not the mean and variance of the truncated distribution.

Answer 4

Three ways has worked for me:

1. using sample() with rnorm():

   sample(x=min:max, replace= TRUE, rnorm(n, mean))

2. using the msm package and the rtnorm function:

   rtnorm(n, mean, lower=min, upper=max)

3. using the rnorm() and specifying the lower and upper limits, as Hugh has posted above:

   sample <- rnorm(n, mean=mean); sample <- sample[x > min & x < max]