simPois {IHSEP}

Simulate a Poisson process

Description

Simulate an (inhomogeneous) Poisson process with a given intensity/rate function over the interval [0,T]

Usage

simPois(int=NULL,cens=1,int.M=optimize(int,c(0,cens),maximum=TRUE)$obj*1.1,
        B=max(as.integer(sqrt(int.M * cens)),10), exp.int=FALSE,par=c(1,2))

Arguments

int

A (vectorized) positive function. The intensity/rate function.

cens

A positive scalar. The censoring time, or time of termination of observations, T.

int.M

A positive scalar. Maximum value of the intensity function over [0,T], or a larger value.

B

An integer. The block size to be used in generating exponential random variables in blocks.

exp.int

Logical. Set to TRUE, if the intensity function is exponential, i.e. a*exp(-b*x). If set to TRUE, the parameters a and b should also be supplied via the argument par.

par

A numerical vector of two elements, giving the parameters a and b of the exponential intensity function. The values are not ignored if exp.int is set to FALSE.

Details

The function works by first generating a Poisson process with constant rate int.M oever [0,T], and then thinning the process with retention probability function

p(x)=\code{int}(x)/\code{int.M}

.

When generating the homoneous Poisson process, it first generates about \Lambda+1.96*\sqrt{\Lambda} exponential variables, then, if the exponential variables are not enough (their sum has not reached the censoring time T yet), generates exponential variables in blocks of size B until the total of all the generated exponentials exceeds T. Then cumsums of the exponentials that are smaller than or equal to T are retained as the event times of the homogeneous Poisson process. This method apparantly does not produce tied event times.

Value

A numerical vector giving the event/jump times of the Poisson process in [0,T], in ascending order.

Author(s)

Feng Chen <[e-mail]>

See Also

simPois0

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (int, cens = 1, int.M = optimize(int, c(0, cens), maximum = TRUE)$obj * 
    1.1, B = max(as.integer(sqrt(int.M * cens)), 10)) 
{
    tms <- rexp(as.integer(int.M * cens + 2 * sqrt(int.M * cens)), 
        rate = int.M)
    while (sum(tms) < cens) tms <- c(tms, rexp(B, rate = int.M))
    cumtms <- cumsum(tms)
    tms <- cumtms[cumtms <= cens]
    N <- length(tms)
    if (N == 0) 
        return(numeric(0))
    tms[as.logical(mapply(rbinom, n = 1, size = 1, prob = int(tms)/int.M))]
  }