17
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Why do you want to do this? If it's bounded then it can't really be normal. What are you trying to achieve? $\endgroup$ – gung - Reinstate Monica Aug 26 '14 at 1:34
  • $\begingroup$ x <- rnorm(n, mean, sd); x <- x[x > lower.limit & x < upper.limit] $\endgroup$ – Hugh Aug 26 '14 at 10:31
  • 3
    $\begingroup$ @Hugh that's great ... as long as you don't care how many random values you get. $\endgroup$ – Glen_b -Reinstate Monica Aug 26 '14 at 10:45
31
$\begingroup$

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

    enter image description here

    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.

    enter image description here

  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]

    enter image description here

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

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Why bother doing this? It is so simple to generate normal random numbers and drop those that need truncation that it isn't necessary to be complicated about it unless the truncation desired is a close to 100% of the area of the density. $\endgroup$ – Carl Nov 1 '17 at 0:10
  • 2
    $\begingroup$ Perhaps I'm misinterpreting the original question. I came across this question while trying to figure out how to achieve a not-directly-stats-related programming task in R, and I've only now noticed this page is a statistics stackexchange, not a programming stackexchange. :) In my case, I wanted to generate a specific amount of random integers, with values ranging from 0 to 100, and I wanted the generated values to fall on a nice bell curve across that range. Since writing this I've realized that sample(x=min:max, prob=dnorm(...)) is maybe an easier way to do that. $\endgroup$ – Aaron Wells Nov 2 '17 at 2:57
  • $\begingroup$ @Glen_b Aaron Wells mentions sample(x=min:max, prob=dnorm(...)) which seems a bit shorter than your answer. $\endgroup$ – Carl Nov 2 '17 at 3:04
  • $\begingroup$ But note that the sample() trick is only useful if you're trying to pick random integers, or some other set of discrete, predefined values. $\endgroup$ – Aaron Wells Nov 14 '17 at 0:21
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.