[R] Viewing function source

[(Ted Harding)]

On 26-Aug-03 Thomas Lumley wrote:
> You can generate Poisson random numbers from a Poisson process like
> this:
> 
> rfishy<-function(lambda){
>     t <- 0
>     i <- -1
>     while(t<=lambda){
>         t<-t-log(runif(1))
>         i<-i+1
>     }
>     return(i)
> }

You can of course easily vectorise this:

  rfishy<-function(n,lambda){
    t<-0*rep(1,n)
    i<-(-1)*rep(1,n)
    while(any(t<=lambda)){
      u<-(t<=lambda)
      t<-t-log(runif(n))*u
      i<-i+1*u
    }
    return(i)
  }

> The name of this generator is descriptive, not a pub. It is very slow
> for large lambda, and incorrect for extremely large lambda (and
> possibly for extremely small lambda).

Not sure why it should be incorrect. It's certainly not theoretically
incorrect. Rounding errors? Problem with runif()?

> If you only wanted a fairly small number of random variates with, say,
> 1e-6<lambda<100, then it's not too bad.

In the above vectorised form, if it's fast for n=1 it stays pretty fast
for large n: try it with z<-rfishy(10000,5); even rfishy(100000,5) only
takes a few seconds.

Best wishes,
Ted.




--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 26-Aug-03                                       Time: 23:38:09
------------------------------ XFMail ------------------------------