[R] Simulation using parts of density function

[Prof Brian Ripley]

Please do not send everything twice: you are using R-help in both the To: 
and Cc: fields.

I disagree with Ted: it _is_ much easier to create a generator for this 
purpose.

Consider

rtgamma <- function(n, ..., tr = log(5000000))
{
     p <- pgamma(tr, ...)
     qgamma(p*runif(n), ...)
}

as inversion (especially at C level) is plenty fast enough.


On Wed, 2 May 2007, Thür Brigitte wrote:

>
> Thanks for your code! It is not exactly what I really want - but it is my fault, because my description was wrong...
>
> It is not "sim" but rhater exp(rgamma(...)) that should not exceed 5000000. So I tried to modify your code but it doesn't really work. "sim.test" returns just 1 value and not 999.
>
> My modified code:
>
> sim.test <- NULL
> for(i in 1:999)
> {sim<-NULL
>  remain <- rpois(1,2000)
>  x <- remain
>  while(remain>0){
>    sim0<-replicate(10*remain,
>       exp(rgamma(1, scale = 0.5, shape = 12))
>       )
>    sim<-c(sim,sim0[sim0<=5000000])
>    remain<-(x - length(sim))
>  }
>  sim<-sim[1:x]}
> sim.test <- rbind(sim.test, c(value=sum(sim)))
>
> Thanks for any help,
> Brigitte
>
>
> -----Ursprüngliche Nachricht-----
> Von: ted.harding at nessie.mcc.ac.uk [mailto:ted.harding at nessie.mcc.ac.uk]
> Gesendet: Dienstag, 1. Mai 2007 20:18
> An: Thür Brigitte
> Cc: r-help at stat.math.ethz.ch
> Betreff: RE: [R] Simulation using parts of density function
>
> On 01-May-07 17:03:46, Thür Brigitte wrote:
>>
>> Hi
>>
>> My simulation with the followin R code works perfectly:
>> sim <- replicate(999, sum(exp(rgamma(rpois(1,2000),
>> scale = 0.5, shape = 12))))
>>
>> But now I do not want to have values in object "sim" exceeding
>> 5'000'000, that means that I am just using the beginning of
>> densitiy function gamma x < 15.4. Is there a possibility to
>> modify my code in an easy way?
>>
>> Thanks for any help!
>>
>> Regards, Brigitte
>
> A somewhat extreme problem!
>
> The easiest way to modify the code is as below -- certiainly easier
> than writing a special function to draw random samples from the
> truncated gamma distribution.
>
> A bit of experimentation shows that, from your code above, about
> 10% of the results are <= 5000000. So:
>
>  sim<-NULL
>  remain <- 999
>  while(remain>0){
>    sim0<-replicate(10*remain,
>       sum(exp(rgamma(rpois(1,2000), scale = 0.5, shape = 12)))
>       )
>    sim<-c(sim,sim0[sim0<=5000000])
>    remain<-(999 - length(sim))
>  }
>  sim<-sim[1:999]
>
> Results of a run:
>
>  sum(sim>5000000)
>  [1] 0
>
>  max(sim)
>  [1] 4999696
>
>  length(sim)
>  [1] 999
>
> It may be on the slow side (though not hugely -- on a quite slow
> machine the above run was completed in 2min 5sec, while the
> 999-replicate in your original took 15sec. So about 8 times as long.
> Most of this, of course, is taken up with the first round.
>
> Hoping this helps,
> Ted.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595