[R] truncated normal

[Duncan Murdoch]

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 at stats.uwo.ca> 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 at stats.uwo.ca <mailto:
>>> murdoch at stats.uwo.ca>> 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
>>>
>>>
>>>
>