thinning: Thinning Algorithm Visualization

Description

This function animates the "thinning" approach the generation of the random event times for a non-homogeneous Poisson process with a specified intensity function, given a majorizing function that dominates the intensity function.

Usage

thinning(
  maxTime = 24,
  intensityFcn = function(x) (5 - sin(x/0.955) - (4 * cos(x/3.82)))/0.5,
  majorizingFcn = NULL,
  majorizingFcnType = NULL,
  seed = NA,
  maxTrials = Inf,
  plot = TRUE,
  showTitle = TRUE,
  plotDelay = plot * -1
)

Value

returns a vector of the generated random event times

Arguments

maxTime: maximum time of the non-homogeneous Poisson process. (The minimum time is assumed to be zero.)
intensityFcn: intensity function corresponding to rate of arrivals across time.
majorizingFcn: majorizing function. Default value is NULL, corresponding to a constant majorizing function that is 1.01 times the maximum value of the intensity function. May alternatively be provided as a user-specified function, or as a data frame requiring additional notation as either piecewise-constant or piecewise-linear. See examples.
majorizingFcnType: used to indicate whether a majorizing function that is provided via data frame is to be interpreted as either piecewise-constant ("pwc") or piecewise-linear ("pwl"). If the majorizing function is either the default or a user-specified function (closure), the value of this parameter is ignored.
seed: initial seed for the uniform variates used during generation.
maxTrials: maximum number of accept-reject trials; infinite by default.
plot: if TRUE, visual display will be produced. If FALSE, generated event times will be returned without visual display.
showTitle: if TRUE, display title in the main plot.
plotDelay: wait time, in seconds, between plots; -1 (default) for interactive mode, where the user is queried for input to progress.

Details

There are three modes for visualizing Lewis and Shedler's thinning algorithm for generating random event times for a non-homogeneous Poisson process with a particular intensity function:

interactive advance (plotDelay = -1), where pressing the 'ENTER' key advances to the next step (an accepted random variate) in the algorithm, typing 'j #' jumps ahead # steps, typing 'q' quits immediately, and typing 'e' proceeds to the end;
automatic advance (plotDelay > 0); or
final visualization only (plotDelay = 0).

As an alternative to visualizing, event times can be generated

References

Lewis, P.A.W. and Shedler, G.S. (1979). Simulation of non-homogeneous Poisson processes by thinning. Naval Research Logistics, 26, 403–413.

Examples

Run this code


nhpp <- thinning(maxTime = 12, seed = 8675309, plotDelay = 0)
nhpp <- thinning(maxTime = 24, seed = 8675309, plotDelay = 0)

# \donttest{
nhpp <- thinning(maxTime = 48, seed = 8675309, plotDelay = 0)

# thinning with custom intensity function and default majorizing function
intensity <- function(x) { 
    day <- 24 * floor(x/24)
    return(80 * (dnorm(x, day + 6,    2.5) + 
                 dnorm(x, day + 12.5, 1.5) + 
                 dnorm(x, day + 19,   2.0)))
}
nhpp <- thinning(maxTime = 24, plotDelay = 0, intensityFcn = intensity)

# thinning with custom intensity and constant majorizing functions
major <- function(x) { 25 }
nhpp <- thinning(maxTime = 24, plotDelay = 0, intensityFcn = intensity,
                 majorizingFcn = major)

# piecewise-constant data.frame for bounding default intensity function
fpwc <- data.frame(
    x = c(0, 2, 20, 30, 44, 48),
    y = c(5, 5, 20, 12, 20,  5)
)
nhpp <- thinning(maxTime = 24, plotDelay = 0, majorizingFcn = fpwc, majorizingFcnType = "pwc")

# piecewise-linear data.frame for bounding default intensity function
fpwl <- data.frame(
    x = c(0, 12, 24, 36, 48),
    y = c(5, 25, 10, 25, 5)
)
nhpp <- thinning(maxTime = 24, plotDelay = 0, majorizingFcn = fpwl, majorizingFcnType = "pwl")

# piecewise-linear closure/function bounding default intensity function
fclo <- function(x) {
    if (x <= 12) (5/3)*x + 5
    else if (x <= 24) 40 - (5/4)*x
    else if (x <= 36) (5/4)*x - 20
    else 85 - (5/3) * x
}
nhpp <- thinning(maxTime = 48, plotDelay = 0, majorizingFcn = fclo)

# thinning with fancy custom intensity function and default majorizing
intensity <- function(x) { 
    day <- 24 * floor(x/24)
    return(80 * (dnorm(x, day + 6,    2.5) + 
                 dnorm(x, day + 12.5, 1.5) + 
                 dnorm(x, day + 19,   2.0)))
}
nhpp <- thinning(maxTime = 24, plotDelay = 0, intensityFcn = intensity)

# piecewise-linear data.frame for bounding custom intensity function
fpwl <- data.frame(
    x = c(0,  6,  9, 12, 16, 19, 24, 30, 33, 36, 40, 43, 48),
    y = c(5, 17, 12, 28, 14, 18,  7, 17, 12, 28, 14, 18,  7)
)
nhpp <- thinning(maxTime = 48, plotDelay = 0, intensityFcn = intensity,
          majorizingFcn = fpwl, majorizingFcnType = "pwl")

# thinning with simple custom intensity function and custom majorizing
intensity <- function(t) {
    if      (t < 12) t
    else if (t < 24) 24 - t
    else if (t < 36) t - 24
    else             48 - t
}
majorizing <- data.frame(
    x = c(0, 12, 24, 36, 48),
    y = c(1, 13,  1, 13,  1))
times <- thinning(plotDelay = 0, intensityFcn = intensity,
    majorizingFcn = majorizing , majorizingFcnType = "pwl", maxTime = 48)
# }

Run the code above in your browser using DataLab