Re: [R] truncated normal

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Wed, 23 Jul 2008 16:34:06 -0400

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



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 23 Jul 2008 - 20:55:36 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 23 Jul 2008 - 21:32:54 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive