truncated normal

[cindy Guo]

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?
Thank you.

        [[alternative HTML version deleted]]

______________________________________________
R-help@... 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.

[Duncan Murdoch-2]

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

______________________________________________
R-help@... 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.

[Duncan Murdoch-2]

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@... 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.

[Duncan Murdoch-2]

On 23/07/2008 8:17 PM, cindy Guo wrote:
> Yes, I know. I mean if I want to generate 100 numbers from
> N(0,1)I((0,1),(5,10)). There are two intervals (0,1) and (5,10). Then the
> function will give 50 numbers in the first interval and 50 in the other.

No, it doesn't handle that case at all.  I didn't realize that's what
you wanted.  Sorry.

I've posted this back to the mailing list; maybe someone else will help
you.  (This is one reason it's not a good idea to take discussions
offline:  someone else might have recognized my misunderstanding earlier.)

Duncan Murdoch

> Cindy
>
>
> On 7/23/08, Duncan Murdoch <murdoch@...> wrote:
>> cindy Guo wrote:
>>
>>> Dear Dr. Murdoch,
>>>  First thank you very much for your reply. But it seems to me that it will
>>> generate a point with equal probability in the two intervals. I mean if I
>>> generage 100 points, then 50 will be in one interval and the other 50 in the
>>> other interval. I think more points should be in the interval closer to the
>>> mean. What do you think?
>>>
>> I think you mean something different from me when you say a "truncated"
>> distribution.  My function would generate 100 points in whatever intervals
>> you specified.
>>
>> Duncan Murdoch
>>
>>>  Thanks a lot,
>>>  Cindy
>>>
>>>  On 7/23/08, *Duncan Murdoch* <murdoch@... <mailto:
>>> murdoch@...>> wrote:
>>>
>>>    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@... 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.

[Berwin A Turlach]

G'day all,

On Wed, 23 Jul 2008 20:48:59 -0400
Duncan Murdoch <murdoch@...> wrote:

> On 23/07/2008 8:17 PM, cindy Guo wrote:
> > Yes, I know. I mean if I want to generate 100 numbers from
> > N(0,1)I((0,1),(5,10)). There are two intervals (0,1) and (5,10).
> > Then the function will give 50 numbers in the first interval and 50
> > in the other.

I can think of two strategies:

1) Use rtnorm from the msm *package* to sample from the normal
distribution truncated to the interval from your smallest lower
bound and your largest upper bound; and then use rejection sampling.

I.e. for your original example N(0,1)I((0,1),(2,4)) sample from a
N(0,1)I(0,4) distribution and reject all observations that are between
1 and 2, sample sufficiently many points until you have kept the
required number.  But this could be rather inefficient for your second
example N(0,1)I((0,1),(5,10)).

2) If I am not mistaken, truncating the normal distribution to more than
one interval essentially creates a mixture distribution where the
components of the mixture are normal distributions truncated to a
single interval.  The weights of the mixture are given by the relative
probability with which an observation from a normal distribution falls
into each of the intervals.  Thus, an alternative strategy, exemplified
using your second example, would be (code untested):

samp1 <- rtnorm(100,mean=0,sd=1, lower=0,upper=1)
samp2 <- rtnorm(100,mean=0,sd=1, lower=5,upper=10)
p1 <- pnorm(1)-pnorm(0)
p2 <- pnorm(10)-pnorm(5)  ## Take heed about Duncan's remark on
                          ## evaluating such quantities
p <- p1/(p1+p2)
ind <- runif(100) < p
finalsample <- samp1
finalsample[!ind] <- samp2[!ind]

HTH.

Cheers,

        Berwin

=========================== Full address =============================
Berwin A Turlach                            Tel.: [phone] (secr)
Dept of Statistics and Applied Probability        [phone] (self)
Faculty of Science                          FAX : [phone]      
National University of Singapore    
6 Science Drive 2, Blk S16, Level 7          e-mail: statba@...
Singapore 117546                    http://www.stat.nus.edu.sg/~statba

______________________________________________
R-help@... 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.