[R] Simulation using parts of density function
Thür Brigitte
Brigitte.Thuer at swica.ch
Wed May 2 09:08:40 CEST 2007
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.
--------------------------------------------------------------------
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 ------------------------------
*** eSafe at SWICA scanned this email for malicious content and found it to be
clean ***
More information about the R-help
mailing list