Skip to content

truncated normal

1 message · David Brahm

#
Scott Summerill <scott.ctr.summerill at tc.faa.gov> wrote:
I asked as similar question on S-news a while ago, and this is what I put
together from the replies.  Note "sd" is the standard deviation of the
underlying normal, not of the final truncated distribution.  "lo" and "hi" are
the cutoff points on either side.


rnbnd <- function(N, mn=0, sd=1, lo=-Inf, hi=Inf, tol=3) {
  if (hi-mn < tol*sd || mn-lo < tol*sd)               # Alan Zaslavsky's method
    return(mn + sd * qnorm(runif(N, pnorm((lo-mn)/sd), pnorm((hi-mn)/sd))))
  x <- rep(NA, N)                                     # Rejection method
  bnd <- function(z, lo, hi) ifelse(z<lo | z>hi, NA, z)
  while (any(q <- is.na(x))) x[q] <- bnd(rnorm(sum(q), mn, sd), lo, hi)
  x
}


For efficiency, this function chooses between 2 methods.  If the truncation is
"severe", I generate the distribution from qnorm on a uniform distribution
(thanks originally to Alan Zaslavsky, and this was Thomas Lumley's approach as
well).  If the truncation is not severe (more than 3 sigma out), I use the
rejection method, just throwing out points in the truncated tails.  Good luck!