simPois {IHSEP} | R Documentation |
Simulate a Poisson process
Description
Simulate an (inhomogeneous) Poisson process with a given intensity/rate
function over the interval
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, |
int.M |
A positive scalar. Maximum value of the intensity function over |
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 , and then thinning the process
with retention probability function
.
When generating the homoneous Poisson process, it first generates
about exponential variables, then,
if the exponential variables are not enough (their sum has not reached
the censoring time
yet), generates exponential variables in
blocks of size
B
until the total of all the generated
exponentials exceeds . Then
cumsum
s of the exponentials
that are smaller than or equal to 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 , in ascending order.
Author(s)
Feng Chen <feng.chen@unsw.edu.au>
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))]
}