2
$\begingroup$

I am working on a project where I need to generate random numbers for a given task time which is normally distributed with mean = 40, and standard deviation = 150.

Because of the high SD, I will get some negative values and low values when I generated numbers directly which is unrealistic.

Is there any way where I can generate random numbers normally distributed with limitations (i.e. bounds)?

Improve this question
$\endgroup$
3
  • 6
    $\begingroup$ Whilst it is certainly quite simple to generate from a truncated normal distribution, if it is unrealistic to have negative values, you should rethink whether this distribution is appropriate. $\endgroup$ – Ben Jul 27 '20 at 11:23
  • 3
    $\begingroup$ To be even clearer: a Normal distribution is not a plausible model for positive data with such a huge SD:mean ratio. Even a truncated Normal is unlikely to be a good one for task times. This avenue of investigation is fruitless and will likely mislead you. Psychometrists tend to observe lognormal distributions of such times, suggesting that would be a better family of distributions to use. $\endgroup$ – whuber♦ Jul 27 '20 at 14:21
  • 1
    $\begingroup$ "Is there any way where I can generate random numbers normally distributed with limitations (i.e. bounds)?" This is unclear because it is not possible, but you might be willing to give up some conditions (e.g. not normally distributed), but what conditions are left over? $\endgroup$ – Sextus Empiricus Jul 27 '20 at 21:47
6
$\begingroup$

This sounds like you want to sample from truncated normal distribution. If you only want to truncate the tails of the distribution (regions with low probability), than the approach suggested by Dave is probably enough. In other cases it might however quickly get inefficient. Better approach was suggested by Christian P. Robert in

Robert, C.P. (1995). Simulation of truncated normal variables. Statistics and Computing 5(2): 121-125.

The algorithm is a bit more complicated, so I suggest you check the paper. Alternatively, if you prefer the code here you can find a C++ implementation from R package extraDistr (disclaimer: it's written by me).

Improve this answer
$\endgroup$
1
2
$\begingroup$

I could imagine something where you use some if/else logic to screen for unrealistic values. There would be some kind of recursion where you keep drawing random numbers u til you’ve gotten 1000 (or whatever) realistic values. Some pseudocode:

i=0
while i < 1000:
    x = make your draw here 
    # (np.random.normal or rnorm, for instance)
    if x is realistic:
        sample[i] = x
        increase i by 1

This will keep drawing random numbers for observation i until it gets a realistic value.

Note that you are not simulating normal data if you do this, as any real number is technically possible for any normal distribution, and you eliminate some values.

Improve this answer
$\endgroup$
8
  • 2
    $\begingroup$ This could get very inefficient, e.g. you want to sample from $\mathcal{N}(\mu, \sigma)$, but truncate it over $[\mu-4\sigma, \mu-3\sigma]$. $\endgroup$ – Tim♦ Jul 27 '20 at 10:26
  • $\begingroup$ Another problem with this approach is that if you take a normal distribution with mean $40$ and standard deviation $150$ and truncate off all the negative values then you will be left with a distribution with a mean close to $135$ and standard deviation close to $98$, which is not what you were aiming for $\endgroup$ – Henry Jul 27 '20 at 20:31
  • $\begingroup$ @Tim What would get so inefficient? $\endgroup$ – Dave Jul 27 '20 at 20:54
  • $\begingroup$ @Dave in above example, probability of observing value from the mentioned interval is 0.001, and to draw 1k samples, you need to make approximately 750k iterations. $\endgroup$ – Tim♦ Jul 27 '20 at 21:02
  • $\begingroup$ @Tim Perhaps there is a mistake in my algorithm, but the goal was to have the occasional unrealistic value trigger a re-draw. Since there would only be a few unrealistic values, there should be 1010 draws, not 750k. Where did I go awry? $\endgroup$ – Dave Jul 27 '20 at 21:05
2
$\begingroup$

Whilst it is certainly quite simple to generate from a truncated normal distribution, if it is unrealistic to have negative values, you should rethink whether this distribution is appropriate. In the present case the lower-bound for the truncation occurs very close to the mean, and so the resulting distribution does not look much like a normal distribution. In any case, in the answer below I show how you can generate values from a truncated normal distribution.

A simple and efficient method to generate random variables from truncated distributions (without discarding generated values) is to generate uniform random variables over the appropriate quantile range, and then use inverse transformation sampling to get the truncated random variables you want. To generate truncated normal random variables with minimum value $x_\min$ and maximum value $x_max$, we first compute the bounds for the quantile:

$$u_\min \equiv \Phi^{-1} \bigg( \frac{x_\min-\mu}{\sigma} \bigg) \quad \quad \quad u_\max \equiv \Phi^{-1} \bigg( \frac{x_\max-\mu}{\sigma} \bigg).$$

We generate the random quantiles $U_1,...,U_n \sim \text{IID U}(u_\min, u_max)$ and we then compute the variables:

$$X_i = \Phi \Big( (\mu + U_i) \sigma \Big).$$

The resulting values $X_1,...,X_n \sim \text{IID TruncN}(\mu, \sigma, x_\min, x_max)$ are lower-truncated normal random variables. Here is a function in R to generate these values.

#Function to generate IID values from truncated normal distribution
rtruncnorm <- function (n, mean = 0, sd = 1, xmin = -Inf, xmax = Inf) {
  
  #Check inputs
  if (!is.numeric(xmin))           { stop('Error: xmin must be numeric') }
  if (!is.vector(xmin))            { stop('Error: xmin must be a single number') }
  if (length(xmin) != 1)           { stop('Error: xmin must be a single number') }
  if (!is.numeric(xmax))           { stop('Error: xmax must be numeric') }
  if (!is.vector(xmax))            { stop('Error: xmax must be a single number') }
  if (length(xmax) != 1)           { stop('Error: xmax must be a single number') }
  if (xmin > xmax)                 { stop('Error: xmin cannot be larger than xmax') }
  
  #Generate random quantiles
  UMIN <- pnorm(xmin, mean = mean, sd = sd);
  UMAX <- pnorm(xmax, mean = mean, sd = sd);
  RAND <- runif(n = n, min = UMIN, max = UMAX);
  
  #Compute output variables
  OUT  <- qnorm(RAND, mean = mean, sd = sd, log = FALSE);
  OUT; }

In your problem you have paramaters $\mu = 40$, $\sigma = 150$, $x_\min = 0$ and $x_\max = \infty$, so here is an example of some generated values:

set.seed(1);
VALUES <- rtruncnorm(100, mean = 40, sd = 150, xmin = 0);
VALUES;
[1]  60.947620  85.841377 137.204732 278.994441  46.359263 271.360050 314.793282 163.444809 153.594112  14.403486
[11]  47.336980  40.640677 172.019700  88.705945 202.535673 116.956365 182.604309 427.443068  87.731151 205.683304
[21] 303.456097  48.742383 160.555002  29.026587  61.341039  89.188375   3.148357  88.294612 251.922414  78.320798
[31] 112.930065 144.830053 115.879226  42.838959 228.478780 165.909662 212.907688  25.004323 184.797715  95.271195
[41] 225.300966 159.110693 208.001017 131.706871 125.391533 210.766031   5.475137 111.690569 187.948715 173.943869
[51] 111.789944 246.811168 101.857892  56.198740  16.454187  23.064491  72.686603 122.442020 163.830618  94.190402
[61] 282.852262  67.424081 107.090779  76.453855 160.302723  59.227644 112.026195 201.098166  19.574651 255.456402
[71]  78.020904 234.685219  79.811798  76.777400 111.466500 266.837916 248.669289  90.119810 205.631169 337.054094
[81] 101.007310 180.790027  92.533455  74.805628 197.410947  46.589354 180.298258  28.145057  56.356945  33.072365
[91]  55.016521  13.744835 157.628333 256.063229 206.300052 214.271779 106.138756  94.981436 220.480028 146.394828

Note that when you lower-truncate so close to the mean, this means that the actual mean and standard deviation of the truncated distribution are substantially different to the pre-truncation parameter values. If you want your post-truncation mean and standard deviation to be equal to your specified parameters, you would need to change the pre-truncation values of $\mu$ and $\sigma$.

Improve this answer
$\endgroup$
16
  • 1
    $\begingroup$ I am curious about these "check inputs" tests that have appeared in your recent answers, because they don't seem to be either revealing or necessary. The best tests will catch errors that R might not detect or might not report on clearly. For instance, R will clearly tell you what the problem is if you supply a non-numerical value for n. Consider focusing your efforts on tests that help the user understand or run the code. $\endgroup$ – whuber♦ Jul 27 '20 at 13:42
  • $\begingroup$ Fair enough --- I add the checks on functions so that the user gets a clear error message when they give an incorrect input. I accept that there will sometimes be unecessary duplication when the internal functions already have suitable error messages. In my experience, the additional work this requires is often repaid by the fact that it makes debugging easier, and separating it into a section allows it to be easily ignored when parsing the code, so I don't really consider it too much of a barrier to understanding. In any case, I'm happy to hear your view on the matter. $\endgroup$ – Ben Jul 27 '20 at 13:45
  • 2
    $\begingroup$ I would like to suggest that the code will be just fine without any of the remaining checks, either: R will issue suitable warning and error messages, even when xmax is less than xmin (which, arguably, is the most informative check for the reader). There's no harm in doing all these checks, usually (they would impair the performance of any procedure that would be called many times), but it struck me that having a great many of them could distract from the statistical ideas you are conveying with the code. $\endgroup$ – whuber♦ Jul 27 '20 at 14:00
  • 1
    $\begingroup$ @whuber Indeed, I can not find any example exceeding 1. I believe (intuitively) that it is because if we cut a smaller part of the left tail, then we are approaching an exponential distribution for which the CV is 1. I have done some computations in my answer here stats.stackexchange.com/a/479394 Wherever you truncate the distribution, you do not get a CV exceeding 1. $\endgroup$ – Sextus Empiricus Jul 28 '20 at 15:14
  • 1
    $\begingroup$ I see now that in my comment from 8 hours ago I wrote " the ratio of the mean and standard deviation of a truncated normal distribution (at the point X=0) may at most be 1" this should be sd/mean<1 and not mean/sd<1 $\endgroup$ – Sextus Empiricus Jul 28 '20 at 15:27
0
$\begingroup$

While your question is not entirely clear about what you are trying to achieve (how do you wish to go from a Gaussian distribution to a distribution that is truncated at 0?)...

... I thought that it was interesting to show something about the limit of the ratio between the mean and standard deviation of a Gaussian distribution that is truncated at $x=0$. (this issue has been mentioned in some of the comments)

Below is a piece of code and a graph that shows the computation of Gaussian distributions that have been truncated at different z-values (and shifted and re-scaled appropriately in order to have the truncating occur at $x=0$ and have the population mean equal to $\bar{x} = 40$).

What we can notice is that by changing the point where we truncate the distribution we can shift from a curve that looks like a Gaussian distribution (when we cut at a low z-value) to a curve that approaches an exponential distribution (when we cut at a high z-value and only have the right tail, which approximates an exponential function).

From this display I guess, intuitively, that the ratio of the standard deviation and the mean for this truncated distribution, is not able to become larger than this ratio for an exponential distribution (for an exponential distribution this ratio is 1).

Therefore: By truncating a normal distribution such that no negative values appear, we can not get a distribution whose standard deviation is larger than it's mean. (and you are looking for sd = 150 and mean = 40, which means that truncating a normal distribution won't be able to do it)

example of different ways to truncate

library(truncnorm)

x = seq(-10^3,10^3,0.1)

### empty canvas/plot
plot(-100,-100, 
     ylim = c(0,0.025), xlim = c(0,200),
     xlab = "x", ylab = "density")

d = 20 ### number of curves
i = 0  ### counter used in for-loop

varst = rep(0,d-1)

for (trunc in qnorm(seq(1/d,1-1/d,1/d))) {
  
  ### compute truncated standard normal
  ### and it's mean and variance
  y <- dtruncnorm(x, mean = 0, sd = 1, a = trunc)
  mean = dnorm(trunc)/(1-pnorm(trunc))
  var  = (1+trunc*dnorm(trunc)/(1-pnorm(trunc)) - mean^2)
  
  ### transform such that the mean is equal to 40
  xtrans <- (x-trunc)*40/(mean-trunc)
  ytrans <- y/(40/(mean-trunc))
  
  
  ### storing variance of transformed trucated standard normal (multipliying with square of scalefactor)
  varst[i+1] = var*(40/(mean-trunc))^2
  
  ### plot
  lines(xtrans[xtrans>=0],ytrans[xtrans>=0], 
        col = hsv(0.15+i/2/d,1-i/2/d,1-(d-i)/4/d,1))
  i = i+1
}

### exponential distribution
lines(x[x>=0],dexp(x,rate=1/40)[x>=0], lty = 2)

i = 1:(d-1)
legend(200,0.025,xjust = 1, cex = 0.7,
       legend = c("exponential distribution", "normal distribution cut at 5%", "normal distribution cut at 95%"),
       lty = c(2,1,1), col = c(1, hsv(0.15+i/2/d,1-i/2/d,1-(d-i)/4/d,1)[c(1,d-1)]))

For the equations used to compute the mean and variance of the truncated normal distribution see: https://en.wikipedia.org/wiki/Truncated_normal_distribution

Improve this answer
$\endgroup$
0

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.