A number of fantastical solutions have been proposed to the problem of
generating data from a truncated normal distribution. Brian Ripley
and Nick Ellis have put forward one of the CORRECT solutions, i.e. to
draw from an untruncated normal and reject values outside the range
of truncation.
Brian and Nick's solution is the simplest, and will be useful and even
fast if the range of truncation is not too narrow. It begins to break
down as you truncate to a range with small probability content under
the untruncated distribution, since then you will have a correspondingly
small acceptance rate.
The other general approach to this problem is to generate from the
appropriately truncated uniform and use the inverse probability function
to transform back to the truncated range of the desired distribution.
For example:
trunc.rnorm <- function(n,mean,sd,lower,upper)
{
zlower <- (lower-mean)/sd ## transform bounds to z-scores
zupper <- (upper-mean)/sd
u <- runif(n,pnorm(zlower),pnorm(zupper)) ## truncated uniform
mean + sd * qnorm(u) ## transform to normal
}
The computations here are slower per draw than rnorm(), but nothing is
discarded. If you go far into the tail, the numerical approximations in
pnorm() and qnorm() begin to break down and then you need to use other
approximations.
This method obviously generalizes to any distribution for which you have
a "p---" and "q---" function.
|