[R] truncated normal

[Duncan Murdoch]

On 7/23/2008 4:22 PM, Duncan Murdoch wrote:
> On 7/23/2008 3:41 PM, cindy Guo wrote:
>> Hi, I want to generate random samples from truncated normal say
>> Normal(0,1)Indicator((0,1),(2,4)). It has more than one intervals. In the
>> library msm, it seems to me that the 'lower' and 'upper' arguments can only
>> be a number. I tried rtnorm(1,mean=0,sd=1, lower=c(0,2),upper=c(1,4)) and it
>> didn't work. Can you tell me  how I can do truncated normal at more than one
>> intervals?
> 
> The inverse CDF method will work.  For example, this untested code:
> 
> rtruncnorm <- function(n, mean=0, sd=1, lower=-Inf, upper=Inf) {
>    plower <- pnorm(lower, mean, sd) # get the probability values for the
>                                     # lower limits
>    pupper <- pnorm(upper, mean, sd) # ditto for the upper limits
>    p <- plower + (pupper - plower)*runif(n) # get random values between
>                                             # those
>    qnorm(p, mean, sd)  # invert the CDFs
> }
> 
> As I said, this is untested, but it should work if all of mean, sd, 
> lower, and upper are the same length, and in some cases where they aren't.
> 
> Duncan Murdoch

One case where the code above will *not* work is if your truncation is 
too far out in the tails, e.g. a standard normal, truncated to be in the 
interval [100, 101].  It is possible to do those in R, but it's tricky 
to avoid rounding error:  in the example above, both plower and pupper 
will be exactly 1 in this case.  You need to do calculations on a log 
scale, and work with upper tails instead of lower ones.  It all gets 
kind of messy.

Duncan Murdoch