thinning {simEd}

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
)

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

Value

returns a vector of the generated random event times

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

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


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)