Simulating a Poisson Process in R by calling C Code over .Call

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Simulating a Poisson Process in R by calling C Code over .Call

Fabian Zäpernick
Hi

I want to write a C function for the R Code below and call it with .Call:

SimPoisson <- function(lambda,tgrid,T2M)
#Simulation eines Poissonprozesses
{
        NT <- 0
        Ni <- rep(0,length(tgrid))
        tau <- 0
        sign <- 0
        if(lambda != 0)
        {
                i=1
                j=1
                while (1)
                {
    EVar <- rexp(1,lambda)
                        sign <- sign + EVar
                        if (sign > T2M)
                        {
                                break
                        }
                        tau[i] <- sign
                        i = i+1
                        for (j in j:length(tgrid))
                        {
                                if (tgrid[j] < sign)
                                {
                                        Ni[j] <- NT
                                }else
                                {
                                        break
                                }
                        }
                        NT <- NT + 1
                }
                for (j in j:length(tgrid))
                {
                        Ni[j] <- NT
                }
  }
        return(list(NT=NT,Ni=Ni,tau=tau))
}

I read the manual "writing R extensions" over and over again, but i have
no idea, how to solve the problem with tau because i dont no the length
of tau at the begining of the function

Fabian

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Simulating a Poisson Process in R by calling C Code over .Call

Thomas Lumley
On Sun, 13 Jun 2010, Fabian Zäpernick wrote:

> Hi
>
> I want to write a C function for the R Code below and call it with .Call:
>
> SimPoisson <- function(lambda,tgrid,T2M)
> #Simulation eines Poissonprozesses
<snip>
> return(list(NT=NT,Ni=Ni,tau=tau))
> }
>
> I read the manual "writing R extensions" over and over again, but i have
> no idea, how to solve the problem with tau because i dont no the length
> of tau at the begining of the function
>

The standard approach is to start off with some reasonable guess at the length of tau and then if the vector fills up, allocate one that is twice as large and copy the current values of tau into the new vector.

In this case you can do better, since you have a very good initial guess for how long tau will be.  If you make the vector of length qpois(0.9999, lambda), then there is a 99.99% chance that it is long enough.  Using lambda+4*sqrt(lambda) is almost as good.

Yet another approach is to take advantage of the memoryless property of the Poisson process: set tau to some reasonable length, then if it fill up, return and call the function again to handle the remainder of the duration.

     -thomas

Thomas Lumley Assoc. Professor, Biostatistics
[hidden email] University of Washington, Seattle

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.