simEd (version 2.0.1)

igamma: Visualization of Random Variate Generation for the Gamma Distribution

Description

Generates random variates from the Gamma distribution by inversion. Optionally graphs the population cumulative distribution function and associated random variates, the population probability density function and a histogram of the random variates, and the empirical cumulative distribution function versus the population cumulative distribution function.

Usage

igamma(
  u = runif(1),
  shape,
  rate = 1,
  scale = 1/rate,
  minPlotQuantile = 0,
  maxPlotQuantile = 0.95,
  plot = TRUE,
  showCDF = TRUE,
  showPDF = TRUE,
  showECDF = TRUE,
  show = NULL,
  maxInvPlotted = 50,
  plotDelay = 0,
  sampleColor = "red3",
  populationColor = "grey",
  showTitle = TRUE,
  respectLayout = FALSE,
  restorePar = TRUE,
  ...
)

Value

A vector of Gamma random variates

Arguments

u

vector of uniform(0,1) random numbers, or NULL to show population figures only

shape

Shape parameter

rate

Alternate parameterization for scale

scale

Scale parameter

minPlotQuantile

minimum quantile to plot

maxPlotQuantile

maximum quantile to plot

plot

logical; if TRUE (default), one or more plots will appear (see parameters below); otherwise no plots appear

showCDF

logical; if TRUE (default), cdf plot appears, otherwise cdf plot is suppressed

showPDF

logical; if TRUE (default), PDF plot appears, otherwise PDF plot is suppressed

showECDF

logical; if TRUE (default), ecdf plot appears, otherwise ecdf plot is suppressed

show

octal number (0-7) indicating plots to display; 4: CDF, 2: PDF, 1: ECDF; sum for desired combination

maxInvPlotted

number of inversions to plot across CDF before switching to plotting quantiles only

plotDelay

delay in seconds between CDF plots

sampleColor

Color used to display random sample from distribution

populationColor

Color used to display population

showTitle

logical; if TRUE (default), displays a title in the first of any displayed plots

respectLayout

logical; if TRUE (default), respects existing settings for device layout

restorePar

logical; if TRUE (default), restores user's previous par settings on function exit

...

Possible additional arguments. Currently, additional arguments not considered.

Author

Barry Lawson (),
Larry Leemis (),
Vadim Kudlay ()

Details

Generates random variates from the Gamma distribution, and optionally, illustrates

• the use of the inverse-CDF technique,

• the effect of random sampling variability in relation to the PDF and CDF.

When all of the graphics are requested,

• the first graph illustrates the use of the inverse-CDF technique by graphing the population CDF and the transformation of the random numbers to random variates,

• the second graph illustrates the effect of random sampling variability by graphing the population PDF and the histogram associated with the random variates, and

• the third graph illustrates effect of random sampling variability by graphing the population CDF and the empirical CDF associated with the random variates.

All aspects of the random variate generation algorithm are output in red by default, which can be changed by specifying sampleColor. All aspects of the population distribution are output in gray by default, which can be changed by specifying populationColor.

The gamma distribution with parameters \code{shape} = \eqn{a} and
\code{scale} = \eqn{s} has density
     \deqn{f(x) = \frac{1}{s^a\, \Gamma(a)} x^{a-1} e^{-x/s}}{
               f(x) = 1/(s^a Gamma(a)) x^(a-1) e^(-x/s)}
for \eqn{x \ge 0}, \eqn{a > 0}, and \eqn{s > 0}.
(Here \eqn{\Gamma(a)}{Gamma(a)} is the function implemented by
R's \code{\link[base:Special]{gamma}()} and defined in its help.)
The population mean and variance are \eqn{E(X) = as}
and \eqn{Var(X) = as^2}.

The algorithm for generating random variates from the gamma distribution is synchronized (one random variate for each random number) and monotone in u. This means that the variates generated here might be useful in some variance reduction techniques used in Monte Carlo and discrete-event simulation.

Values from the u vector are plotted in the cdf plot along the vertical axis as colored dots. A horizontal, dashed, colored line extends from the dot to the population cdf. At the intersection, a vertical, dashed colored line extends downward to the horizontal axis, where a second colored dot, denoting the associated gamma random variate is plotted.

This is not a particularly fast variate generation algorithm because it uses the base R qgamma function to invert the values contained in u.

All of the elements of the u vector must be between 0 and 1. Alternatively, u can be NULL in which case plot(s) of the theoretical PDF and cdf are displayed according to plotting parameter values (defaulting to display of both the PDF and cdf).

The show parameter can be used as a shortcut way to denote plots to display. The argument to show can be either:

• a binary vector of length three, where the entries from left to right correspond to showCDF, showPDF, and showECDF, respectively. For each entry, a 1 indicates the plot should be displayed, and a 0 indicates the plot should be suppressed.

• an integer in [0,7] interpreted similar to the Unix chmod command. That is, the integer's binary representation can be transformed into a length-three vector discussed above (e.g., 6 corresponds to c(1,1,0)). See examples.

Any valid value for show takes precedence over existing individual values for showCDF, showPDF, and showECDF.

If respectLayout is TRUE, the function respects existing settings for device layout. Note, however, that if the number of plots requested (either via show or via showCDF, showPMF, and showECDF) exceeds the number of plots available in the current layout (as determined by prod(par("mfrow"))), the function will display all requested plots but will also display a warning message indicating that the current layout does not permit simultaneous viewing of all requested plots. The most recent plot with this attribute can be further annotated after the call.

If respectLayout is FALSE, any existing user settings for device layout are ignored. That is, the function uses par to explicitly set mfrow sufficient to show all requested plots stacked vertically to align their horizontal axes, and then resets row, column, and margin settings to their prior state on exit.

The minPlotQuantile and maxPlotQuantile arguments are present in order to compress the plots horizontally. The random variates generated are not impacted by these two arguments. Vertical, dotted, black lines are plotted at the associated quantiles on the plots.

plotDelay can be used to slow down or halt the variate generation for classroom explanation.

In the plot associated with the PDF, the maximum plotting height is associated with 125\ that extends above this limit will have three dots appearing above it.

See Also

stats::rgamma

stats::runif, simEd::vunif

Examples

 igamma(0.5, shape = 5, scale = 3)

 set.seed(8675309)
 igamma(runif(10), 3, 2, showPDF = TRUE)

 set.seed(8675309)
 igamma(runif(10), 3, 2, showECDF = TRUE)

 set.seed(8675309)
 igamma(runif(10), 3, 2, showPDF = TRUE, showECDF = TRUE, sampleColor = "blue3")

 set.seed(8675309)
 igamma(runif(10), 3, 2, showPDF = TRUE, showCDF = FALSE)

 igamma(runif(100), 3, 2, showPDF = TRUE, minPlotQuantile = 0.02, maxPlotQuantile = 0.98)

 # plot the PDF and CDF without any variates
 igamma(NULL, 3, 2, showPDF = TRUE, showCDF = TRUE)

 # plot CDF with inversion and PDF using show
 igamma(runif(10), 3, 2, show = c(1,1,0))
 igamma(runif(10), 3, 2, show = 6)

 # plot CDF with inversion and ECDF using show, using vunif
 igamma(vunif(10), 3, 2, show = c(1,0,1))
 igamma(vunif(10), 3, 2, show = 5)

 # plot CDF with inversion, PDF, and ECDF using show
 igamma(vunif(10), 3, 2, show = c(1,1,1))
 igamma(vunif(10), 3, 2, show = 7)

 # plot three different CDF+PDF+ECDF horizontal displays,
 # with title only on the first display
 oldpar <- par(no.readonly = TRUE)
 par(mfrow = c(3,3))  # 3 rows, 3 cols, filling rows before columns
 set.seed(8675309)
 igamma(runif(20), 3, 2, show = 7, respectLayout = TRUE, restorePar = FALSE)
 igamma(runif(20), 3, 2, show = 7, respectLayout = TRUE, restorePar = FALSE, showTitle = FALSE)
 igamma(runif(20), 3, 2, show = 7, respectLayout = TRUE, restorePar = TRUE,  showTitle = FALSE)
 par(oldpar)

 # display animation of all components
 igamma(runif(10), 3, 2, show = 7, plotDelay = 0.1)

 # display animation of CDF and PDF components only
 igamma(runif(10), 3, 2, show = 5, plotDelay = 0.1)

 if (interactive()) {
   # interactive -- pause at each stage of inversion
   igamma(runif(10), 3, 2, show = 7, plotDelay = -1)
 }

  # overlay visual exploration of ks.test results
  oldpar <- par(no.readonly = TRUE)
  set.seed(54321)
  vals <- igamma(runif(10), 3, 2, showECDF = TRUE, restorePar = FALSE)
  D <- as.numeric(ks.test(vals, "pgamma", 3, 2)$statistic)
  for (x in seq(1.20, 1.60, by = 0.05)) {
    y <- pgamma(x, 3, 2)
    segments(x, y, x, y + D, col = "darkgreen", lwd = 2, xpd = NA)
  }
  par(oldpar) # restore original par values, since restorePar = FALSE above