[R] Simulation using parts of density function

[(Ted Harding)]

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.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 01-May-07                                       Time: 19:18:01
------------------------------ XFMail ------------------------------