Nothing
# =================================== ett =====================================
#' Extremal Types Theorem (ETT)
#'
#' A movie to illustrate the extremal types theorem, that is, convergence
#' of the distribution of the maximum of a random sample of size \eqn{n}
#' from certain distributions to a member of the Generalized Extreme Value
#' (GEV) family, as \eqn{n} tends to infinity.
#' Samples of size \eqn{n} are simulated repeatedly from the chosen
#' distribution. The distributions (simulated empirical and true) of the
#' sample maxima are compared to the relevant GEV limit.
#'
#' @param n An integer scalar. The size of the samples drawn from the
#' distribution chosen using \code{distn}. \code{n} must be no smaller
#' than 2.
#' @param distn A character scalar specifying the distribution from which
#' observations are sampled. Distributions \code{"beta"},
#' \code{"cauchy"}, \code{"chisq"}, \code{"chi-squared"},
#' \code{"exponential"}, \code{"f"}, \code{"gamma"}, \code{"gp"},
#' \code{"lognormal"}, \code{"log-normal"}, \code{"ngev"}, \code{"normal"},
#' \code{"t"}, \code{"uniform"} and \code{"weibull"} are recognised, case
#' being ignored.
#'
#' If \code{distn} is not supplied then \code{distn = "exponential"}
#' is used.
#'
#' The \code{"gp"} case uses the \code{\link[revdbayes]{gp}}
#' distributional functions in the
#' \code{\link[revdbayes]{revdbayes}} package.
#'
#' The \code{"ngev"} case is a negated GEV(1 / \eqn{\xi}, 1, \eqn{\xi})
#' distribution, for \eqn{\xi} > 0, and uses the \code{\link[revdbayes]{gev}}
#' distributional functions in the
#' \code{\link[revdbayes]{revdbayes}} package.
#' If \eqn{\xi} = 1 then this coincides with Example 1.7.5 in Leadbetter,
#' Lindgren and Rootzen (1983).
#'
#' The other cases use the distributional functions in the
#' \code{\link[stats]{stats-package}}.
#' If \code{distn = "gamma"} then the \code{(shape, rate)}
#' parameterisation is used. If \code{scale} is supplied via \code{params}
#' then \code{rate} is inferred from this.
#' If \code{distn = "beta"} then \code{ncp} is forced to be zero.
#' @param params A named list of additional arguments to be passed to the
#' density function associated with distribution \code{distn}.
#' The \code{(shape, rate)} parameterisation is used for the gamma
#' distribution (see \code{\link[stats]{GammaDist}}) even if the value of
#' the \code{scale} parameter is set using \code{params}.
#'
#' If a parameter value is not supplied then the default values in the
#' relevant distributional function set using \code{distn} are used,
#' except for
#' \code{"beta"} (\code{shape1 = 2, shape2 = 2}),
#' \code{"chisq"} (\code{df = 4}),
#' \code{"f"} (\code{df1 = 4, df2 = 8}),
#' \code{"ngev"} (\code{shape = 0.2}).
#' \code{"gamma"} (\code{shape = 2},
#' \code{"gp"} (\code{shape = 0.1}),
#' \code{"t"} (\code{df = 4}) and
#' \code{"weibull"} (\code{shape = 2}).
#' @param panel_plot A logical parameter that determines whether the plot
#' is placed inside the panel (\code{TRUE}) or in the standard graphics
#' window (\code{FALSE}). If the plot is to be placed inside the panel
#' then the tkrplot library is required.
#' @param hscale,vscale Numeric scalars. Scaling parameters for the size
#' of the plot when \code{panel_plot = TRUE}. The default values are 1.4 on
#' Unix platforms and 2 on Windows platforms.
#' @param n_add An integer scalar. The number of simulated datasets to add
#' to each new frame of the movie.
#' @param delta_n A numeric scalar. The amount by which n is increased
#' (or decreased) after one click of the + (or -) button in the parameter
#' window.
#' @param arrow A logical scalar. Should an arrow be included to show the
#' simulated sample maximum from the top plot being placed into the
#' bottom plot?
#' @param leg_cex The argument \code{cex} to \code{\link[graphics]{legend}}.
#' Allows the size of the legend to be controlled manually.
#' @param ... Additional arguments to the rpanel functions
#' \code{\link[rpanel]{rp.button}} and
#' \code{\link[rpanel]{rp.doublebutton}}, not including \code{panel},
#' \code{variable}, \code{title}, \code{step}, \code{action}, \code{initval},
#' \code{range}.
#' @details Loosely speaking, a consequence of the
#' \href{https://en.wikipedia.org/wiki/Extreme_value_theory#Univariate_theory}{Extremal Types Theorem}
#' is that, in many situations, the maximum of a \emph{large number}
#' \eqn{n} of independent random variables has \emph{approximately} a
#' GEV(\eqn{\mu, \sigma, \xi)}) distribution, where \eqn{\mu} is a location
#' parameter, \eqn{\sigma} is a scale parameter and \eqn{\xi} is a shape
#' parameter. See Coles (2001) for an introductory account and
#' Leadbetter et al (1983) for greater detail and more examples.
#' The Extremal Types Theorem is an asymptotic result that considers the
#' possible limiting distribution of linearly normalised maxima
#' as \eqn{n} tends to infinity.
#' This movie considers examples where this limiting result holds and
#' illustrates graphically the closeness of the limiting approximation
#' provided by the relevant GEV limit to the true finite-\eqn{n}
#' distribution.
#'
#' Samples of size \code{n} are repeatedly simulated from the distribution
#' chosen using \code{distn}. These samples are summarized using a histogram
#' that appears at the top of the movie screen. For each sample the maximum
#' of these \code{n} values is calculated, stored and added to another plot,
#' situated below the first plot.
#' A \code{\link[graphics]{rug}} is added to a histogram provided that it
#' contains no more than 1000 points.
#' This plot is either a histogram or an empirical c.d.f., chosen using a
#' radio button.
#'
#' The probability density function (p.d.f.) of the original
#' variables is superimposed on the top histogram.
#' There is a checkbox to add to the bottom plot the exact p.d.f./c.d.f. of
#' the sample maxima and an approximate (large \code{n}) GEV p.d.f./c.d.f.
#' implied by the ETT.
#' The GEV shape parameter \eqn{\xi} that applies in the limiting
#' case is used. The GEV location \eqn{\mu} and scale
#' \eqn{\sigma} are set based on constants used to normalise the maxima
#' to achieve the GEV limit.
#' Specifically, \eqn{\mu} is set at the 100(1-1/\eqn{n})\% quantile of the
#' distribution \code{distn} and \eqn{\sigma} at
#' (1 / \eqn{n}) / \eqn{f(\mu)}, where \eqn{f} is the
#' density function of the distribution \code{distn}.
#'
#' Once it starts, four aspects of this movie are controlled by the user.
#' \itemize{
#' \item There are buttons to increase (+) or decrease (-) the sample
#' size, that is, the number of values over which a maximum is
#' calculated.
#' \item Each time the button labelled "simulate another \code{n_add}
#' samples of size n" is clicked \code{n_add} new samples are simulated
#' and their sample maxima are added to the bottom histogram.
#' \item There is a button to switch the bottom plot from displaying
#' a histogram of the simulated maxima, the exact p.d.f. and the
#' limiting GEV p.d.f. to the empirical c.d.f. of the simulated data,
#' the exact c.d.f. and the limiting GEV c.d.f.
#' \item There is a box that can be used to display only the bottom
#' plot. This option is selected automatically if the sample size
#' \eqn{n} exceeds 100000.
#' \item There is a box that can be used to display only the bottom
#' plot. This option is selected automatically if the sample size
#' \eqn{n} exceeds 100000.
#' }
#' For further detail about the examples specified by \code{distn}
#' see Chapter 1 of Leadbetter et al. (1983) and Chapter 3 of
#' Coles (2001). In many of these examples
#' (\code{"exponential"}, \code{"normal"}, \code{"gamma"},
#' \code{"lognormal"}, \code{"chi-squared"}, \code{"weibull"}, \code{"ngev"})
#' the limiting GEV distribution has a shape
#' parameter that is equal to 0. In the \code{"uniform"} case the limiting
#' shape parameter is -1 and in the \code{"beta"} case it is
#' -1 / \code{shape2}, where \code{shape2} is the
#' second parameter of the \code{\link[stats]{Beta}} distribution.
#' In the other cases the limiting shape parameter is positive,
#' with respective values \code{shape}
#' (\code{"gp"}, see \code{\link[revdbayes]{gp}}),
#' 1 / \code{df} (\code{"t"}, see \code{\link[stats]{TDist}}),
#' 1 (\code{"cauchy"}, see \code{\link[stats]{Cauchy}}),
#' 2 / \code{df2} (\code{"f"}, see \code{\link[stats]{FDist}}).
#' @return Nothing is returned, only the animation is produced.
#' @references Coles, S. G. (2001) \emph{An Introduction to Statistical
#' Modeling of Extreme Values}, Springer-Verlag, London.
#' \doi{10.1007/978-1-4471-3675-0_3}
#' @references Leadbetter, M., Lindgren, G. and Rootzen, H. (1983)
#' \emph{Extremes and Related Properties of Random Sequences and Processes.}
#' Springer-Verlag, New York.
#' \doi{10.1007/978-1-4612-5449-2}
#' @seealso \code{\link{movies}}: a user-friendly menu panel.
#' @seealso \code{\link{smovie}}: general information about smovie.
#' @examples
#' # Exponential data: xi = 0
#' ett()
#'
#' # Uniform data: xi =-1
#' ett(distn = "uniform")
#'
#' # Student t data: xi = 1 / df
#' ett(distn = "t", params = list(df = 5))
#' @export
ett <- function(n = 20, distn, params = list(), panel_plot = TRUE, hscale = NA,
vscale = hscale, n_add = 1, delta_n = 1, arrow = TRUE,
leg_cex = 1.25, ...) {
# Check that revdbayes is installed
if (!requireNamespace("revdbayes", quietly = TRUE)) {
stop("the revdbayes package is needed. Please install it.",
call. = FALSE)
}
temp <- set_scales(hscale, vscale)
hscale <- temp$hscale
vscale <- temp$vscale
# To add another distribution
# 1. misc.R: add code to set_fun_args(), set_top_range(), set_leg_pos()
# 2. add lines to rfun, dfun, qfun, pfun
# 3. ett_movie_plot(): add to the_distn and gev_pars
if (!is.wholenumber(n) | n < 2) {
stop("n must be an integer that is no smaller than 2")
}
if (!is.wholenumber(n_add) | n_add < 1) {
stop("n_add must be an integer that is no smaller than 1")
}
if (!is.wholenumber(delta_n) | delta_n < 1) {
stop("delta_n must be an integer that is no smaller than 1")
}
if (!is.list(params)) {
stop("params must be a named list")
}
if (missing(distn)) {
distn <- "exponential"
}
if (distn == "ngev") {
if (!is.null(params$loc) | !is.null(params$scale)) {
warning("In the negated GEV case you cannot set loc or scale")
}
if (!is.null(params$shape)){
if (params$shape <= 0) {
stop("the shape parameter must be positive in the negated GEV case")
}
}
}
xlab <- "x"
#
distn <- tolower(distn)
if (distn == "log-normal") {
distn <- "lognormal"
}
if (distn == "chisq") {
distn <- "chi-squared"
}
# Set the density, distribution, quantile and simulation functions
# "rngev" is included because it is an example that is in the domain of
# attraction of the Gumbel case but the upper endpoint is finite.
#
rfun <-
switch(distn,
"exponential" = stats::rexp,
"uniform" = stats::runif,
"gp" = revdbayes::rgp,
"normal" = stats::rnorm,
"beta" = stats::rbeta,
"t" = stats::rt,
"gamma" = stats::rgamma,
"lognormal" = stats::rlnorm,
"cauchy" = stats::rcauchy,
"chi-squared" = stats::rchisq,
"f" = stats::rf,
"weibull" = stats::rweibull,
"ngev" = rngev,
NULL)
if (is.null(rfun)) {
stop("Unsupported distribution")
}
dfun <-
switch(distn,
"exponential" = stats::dexp,
"uniform" = stats::dunif,
"gp" = revdbayes::dgp,
"normal" = stats::dnorm,
"beta" = stats::dbeta,
"t" = stats::dt,
"gamma" = stats::dgamma,
"lognormal" = stats::dlnorm,
"cauchy" = stats::dcauchy,
"chi-squared" = stats::dchisq,
"f" = stats::df,
"weibull" = stats::dweibull,
"ngev" = dngev)
qfun <-
switch(distn,
"exponential" = stats::qexp,
"uniform" = stats::qunif,
"gp" = revdbayes::qgp,
"normal" = stats::qnorm,
"beta" = stats::qbeta,
"t" = stats::qt,
"gamma" = stats::qgamma,
"lognormal" = stats::qlnorm,
"cauchy" = stats::qcauchy,
"chi-squared" = stats::qchisq,
"f" = stats::qf,
"weibull" = stats::qweibull,
"ngev" = qngev)
pfun <-
switch(distn,
"exponential" = stats::pexp,
"uniform" = stats::punif,
"gp" = revdbayes::pgp,
"normal" = stats::pnorm,
"beta" = stats::pbeta,
"t" = stats::pt,
"gamma" = stats::pgamma,
"lognormal" = stats::plnorm,
"cauchy" = stats::pcauchy,
"chi-squared" = stats::pchisq,
"f" = stats::pf,
"weibull" = stats::pweibull,
"ngev" = pngev)
# Set the arguments to the distributional functions
fun_args <- set_fun_args(distn, dfun, fun_args, params)
# Set sensible scales for the plots
if (distn == "t") {
if (fun_args$df < 2) {
top_p_vec <- c(0.05, 0.95)
bottom_p_vec <- c(0.01, 0.7)
} else if (fun_args$df < 3) {
top_p_vec <- c(0.01, 0.99)
bottom_p_vec <- c(0.01, 0.9)
} else {
top_p_vec <- c(0.001, 0.999)
bottom_p_vec <- c(0.001, 0.999)
}
} else if (distn == "cauchy"){
top_p_vec <- c(0.05, 0.95)
bottom_p_vec <- c(0.01, 0.7)
} else if (distn == "gp") {
if (fun_args$shape > 0.3) {
top_p_vec <- c(0.001, 0.95)
bottom_p_vec <- c(0.001, 0.7)
} else {
top_p_vec <- c(0.001, 0.999)
bottom_p_vec <- c(0.001, 0.999)
}
} else {
top_p_vec <- c(0.001, 0.999)
bottom_p_vec <- c(0.001, 0.999)
}
# Set the range for the top plot
top_range <- set_top_range(distn, p_vec = top_p_vec, fun_args, qfun)
# Set the legend position
leg_pos <- set_leg_pos(distn, fun_args)
top_leg_pos <- leg_pos$top_leg_pos
bottom_leg_pos <- leg_pos$bottom_leg_pos
# Set a unique panel name to enable saving of objects to the correct panel
now_time <- strsplit(substr(date(), 12, 19), ":")[[1]]
now_time <- paste(now_time[1], now_time[2], now_time[3], sep = "")
my_panelname <- paste("ett_", now_time, sep = "")
old_n <- 0
# Create buttons for movie
show_dens <- FALSE
show_dens_only <- FALSE
pdf_or_cdf <- "pdf"
old_y <- NULL
save_last_y <- NULL
sample_maxima <- NULL
ett_panel <- rpanel::rp.control("extremal types theorem",
panelname = my_panelname, n = n,
n_add = n_add, dfun = dfun, qfun = qfun,
rfun = rfun, pfun = pfun,
fun_args = fun_args, distn = distn,
top_range = top_range, top_p_vec = top_p_vec,
bottom_p_vec = bottom_p_vec,
show_dens = show_dens,
show_dens_only = show_dens_only,
pdf_or_cdf = pdf_or_cdf,
top_leg_pos = top_leg_pos,
bottom_leg_pos = bottom_leg_pos,
xlab = xlab, arrow = arrow,
old_n = old_n, old_pdf_or_cdf = pdf_or_cdf,
old_show_dens = show_dens,
old_show_dens_only = show_dens_only,
old_y = old_y, save_last_y = save_last_y,
sample_maxima = sample_maxima,
leg_cex = leg_cex)
#
redraw_plot <- NULL
panel_redraw <- function(panel) {
rpanel::rp.tkrreplot(panel = panel, name = redraw_plot)
# rp.tkrreplot() doesn't update the panel automatically, so do it manually
# Get ...
panel$sample_maxima <- rpanel::rp.var.get(my_panelname, "sample_maxima")
panel$old_n <- rpanel::rp.var.get(my_panelname, "old_n")
panel$old_pdf_or_cdf <- rpanel::rp.var.get(my_panelname, "old_pdf_or_cdf")
panel$old_show_dens <- rpanel::rp.var.get(my_panelname, "old_show_dens")
panel$old_show_dens_only <- rpanel::rp.var.get(my_panelname, "old_show_dens_only")
panel$old_y <- rpanel::rp.var.get(my_panelname, "old_y")
panel$save_last_y <- rpanel::rp.var.get(my_panelname, "save_last_y")
# Put ...
rpanel::rp.control.put(my_panelname, panel)
return(panel)
}
if (panel_plot & !requireNamespace("tkrplot", quietly = TRUE)) {
warning("tkrplot is not available so panel_plot has been set to FALSE.")
panel_plot <- FALSE
}
if (panel_plot) {
rpanel::rp.tkrplot(panel = ett_panel, name = redraw_plot,
plotfun = ett_movie_plot, pos = "right",
hscale = hscale, vscale = vscale, background = "white")
action <- panel_redraw
} else {
action <- ett_movie_plot
}
# Check whether or not n = 1 will work
n_check <- 1
b1 <- do.call(qfun, c(list(p = 1 - 1 / n_check), fun_args))
a1 <- (1 / n_check) / do.call(dfun, c(list(x = b1), fun_args))
if (any(is.infinite(c(a1, b1)))) {
n_lower <- 2
} else {
n_lower <- 1
}
rpanel::rp.doublebutton(panel = ett_panel, variable = n, step = delta_n,
title = "sample size, n",
action = action, initval = n,
range = c(n_lower, NA), showvalue = TRUE, ...)
if (n_add == 1) {
my_title <- paste("simulate another sample")
} else {
my_title <- paste("simulate another", n_add, "samples")
}
dlist <- list(...)
# If the user hasn't set either repeatdelay or repeatinterval then set them
# to the default values in rp.doublebutton (100 milliseconds)
if (is.null(dlist$repeatdelay) & is.null(dlist$repeatinterval)) {
rpanel::rp.button(panel = ett_panel, action = action, title = my_title,
repeatdelay = 100, repeatinterval = 100, ...)
} else {
rpanel::rp.button(panel = ett_panel, action = action, title = my_title,
...)
}
rpanel::rp.radiogroup(panel= ett_panel, pdf_or_cdf, c("pdf", "cdf"),
title = "pdf or cdf in bottom plot",
action = action)
rpanel::rp.checkbox(panel = ett_panel, show_dens,
labels = "show exact and GEV pdf/cdf",
action = action)
rpanel::rp.checkbox(panel = ett_panel, show_dens_only,
labels = "show only exact and GEV pdf/cdf",
action = action)
if (!panel_plot) {
rpanel::rp.do(panel = ett_panel, action = action)
}
return(invisible())
}
# Function to be called by ett().
ett_movie_plot <- function(panel) {
oldpar <- graphics::par(mfrow = c(2, 1), oma = c(0, 0, 0, 0),
mar = c(4, 4, 2, 2) + 0.1)
on.exit(graphics::par(oldpar))
# To please R CMD check
n <- distn <- fun_args <- pdf_or_cdf <- show_dens <- n_add <- rfun <-
qfun <- pfun <- top_range <- dfun <- xlab <- top_leg_pos <- arrow <-
bottom_p_vec <- bottom_leg_pos <- leg_cex <- NULL
panel <- within(panel, {
# Don't simulate very large samples (only show pdfs or cdfs)
if (n > 100000) {
show_dens_only <- TRUE
}
# Don't add the rug in the top plot if n is large
if (n > 1000) {
show_rug <- FALSE
} else {
show_rug <- TRUE
}
if (show_dens_only) {
graphics::par(mfrow = c(1, 1), oma = c(0, 0, 0, 0),
mar = c(4, 4, 2, 2) + 0.1)
} else {
graphics::par(mfrow = c(2, 1), oma = c(0, 0, 0, 0),
mar = c(4, 4, 2, 2) + 0.1)
}
# Set the range of values for the x-axis of the bottom plot
for_qfun <- c(list(p = bottom_p_vec ^ (1 / n)), fun_args)
bottom_range <- do.call(qfun, for_qfun)
# Do the simulation (if required)
if (!show_dens_only) {
sim_list <- c(list(n = n), fun_args)
if (old_pdf_or_cdf == pdf_or_cdf & old_show_dens == show_dens &
old_show_dens_only == show_dens_only) {
temp <- as.matrix(replicate(n_add, do.call(rfun, sim_list)))
max_y <- apply(temp, 2, max)
# Extract the last dataset and the last maximum (for drawing the arrow)
y <- temp[, n_add]
old_y <- y
rm(temp)
last_y <- max_y[n_add]
save_last_y <- last_y
} else {
max_y <- NULL
y <- old_y
last_y <- save_last_y
}
if (n != old_n) {
sample_maxima <- max_y
} else {
sample_maxima <- c(sample_maxima, max_y)
}
}
#
n_x_axis <- 501
# Top plot --------
#
# Set range for x-axis
x <- seq(top_range[1], top_range[2], len = n_x_axis)
# Calculate the density over this range
dens_list <- c(list(x = x), fun_args)
ydens <- do.call(dfun, dens_list)
# Remove any infinite values
finite_vals <- is.finite(ydens)
ydens <- ydens[finite_vals]
x <- x[finite_vals]
# Set the top of the y-axis
ytop <- max(ydens) * 1.2
# Extract the distribution name and parameters
the_distn <-
switch(distn,
"exponential" = paste(distn, "(", fun_args$rate, ")"),
"uniform" = paste(distn, "(", fun_args$min, ",", fun_args$max, ")"),
"gp" = paste("GP", "(", fun_args$loc, ",", fun_args$scale, ",",
fun_args$shape, ")"),
"normal" = paste(distn, "(", fun_args$mean, ",", fun_args$sd, ")"),
"beta" = paste(distn, "(", fun_args$shape1, ",", fun_args$shape2, ")"),
"t" = paste("Student t", "(", fun_args$df, ")"),
"gamma" = paste(distn, "(", fun_args$shape, ",", fun_args$rate, ")"),
"lognormal" = paste(distn, "(", fun_args$meanlog, ",", fun_args$sdlog,
")"),
"cauchy" = paste("Cauchy", "(", fun_args$location, ",", fun_args$scale,
")"),
"chi-squared" = paste(distn, "(", fun_args$df, ",", fun_args$ncp, ")"),
"f" = paste("F", "(", fun_args$df1, ",", fun_args$df2, ",",
fun_args$ncp, ")"),
"weibull" = paste("Weibull", "(", fun_args$shape, ",", fun_args$scale,
")"),
"ngev" = paste("negated GEV", "(", fun_args$loc, ",", fun_args$scale,
",", fun_args$shape, ")")
)
if (!show_dens_only) {
my_xlim <- pretty(c(y, top_range))
my_xlim <- my_xlim[c(1, length(my_xlim))]
# Histogram with rug
graphics::hist(y, col = 8, probability = TRUE, axes = FALSE,
xlab = xlab, ylab = "density", main = "",
xlim = my_xlim, ylim = c(0, ytop))
graphics::lines(x, ydens, xpd = TRUE, lwd = 2, lty = 2)
graphics::axis(2)
graphics::axis(1, line = 0.5)
graphics::title(main = paste(the_distn, ", n = ", n))
graphics::legend(top_leg_pos, legend = expression(f(x)),
col = 1, lwd = 2, lty = 2, box.lty = 0, cex = leg_cex)
u_t <- graphics::par("usr")
if (arrow) {
graphics::segments(last_y, u_t[3], last_y, -10, col = "red", xpd = TRUE,
lwd = 2, lty = 2)
}
if (show_rug) {
graphics::rug(y, line = 0.5, ticksize = 0.05)
}
graphics::rug(last_y, line = 0.5, ticksize = 0.05, col = "red", lwd = 2)
u_t <- my_xlim
}
#
# Bottom plot --------
#
my_xlab <- paste("sample maximum of", n, "values")
bn <- do.call(qfun, c(list(p = 1 - 1 / n), fun_args))
an <- (1 / n) / do.call(dfun, c(list(x = bn), fun_args))
# Set the limiting GEV parameters
gev_pars <-
switch(distn,
"exponential" = list(loc = bn, scale = an, shape = 0),
"uniform" = list(loc = bn, scale = an, shape = -1),
"gp" = list(loc = bn, scale = an, shape = fun_args$shape),
"normal" = list(loc = bn, scale = an, shape = 0),
"beta" = list(loc = bn, scale = an, shape = -1 / fun_args$shape2),
"t" = list(loc = bn, scale = an, shape = 1 / fun_args$df),
"gamma" = list(loc = bn, scale = an, shape = 0),
"lognormal" = list(loc = bn, scale = an, shape = 0),
"cauchy" = list(loc = bn, scale = an, shape = 1),
"chi-squared" = list(loc = bn, scale = an, shape = 0),
"f" = list(loc = bn, scale = an, shape = 2 / fun_args$df2),
"weibull" = list(loc = bn, scale = an, shape = 0),
"ngev" = list(loc = bn, scale = an, shape = 0)
)
for_qgev <- c(list(p = bottom_p_vec), gev_pars)
gev_bottom_range <- do.call(revdbayes::qgev, for_qgev)
if (!show_dens_only) {
bottom_range <- range(c(bottom_range, gev_bottom_range), sample_maxima)
} else {
bottom_range <- range(c(bottom_range, gev_bottom_range))
}
# Set range for x-axis
x <- seq(bottom_range[1], bottom_range[2], len = n_x_axis)
# Calcuate the density over this range
if (pdf_or_cdf == "pdf") {
dens_list <- c(list(x = x), gev_pars)
ygev <- do.call(revdbayes::dgev, dens_list)
} else {
dens_list <- c(list(q = x), gev_pars)
ygev <- do.call(revdbayes::pgev, dens_list)
}
p_list <- c(list(q = x), fun_args, list(log.p = TRUE))
d_list <- c(list(x = x), fun_args)
if (pdf_or_cdf == "pdf") {
temp <- exp((n - 1) * do.call(pfun, p_list))
ytrue <- n * temp * do.call(dfun, d_list)
my_ylab <- "pdf"
# Set the top of the y-axis
if (n > 1){
ytop <- max(ygev, ytrue) * 1.5
}
} else{
ytrue <- exp(n * do.call(pfun, p_list))
my_ylab <- "cdf"
# Set the top of the y-axis
ytop <- 1
}
# Histogram with rug
y <- sample_maxima
if (length(sample_maxima) > 1000) {
show_bottom_rug <- FALSE
} else {
show_bottom_rug <- TRUE
}
if (!show_dens_only) {
my_xlim <- pretty(c(y, bottom_range))
} else {
my_xlim <- pretty(bottom_range)
}
my_xlim <- my_xlim[c(1, length(my_xlim))]
my_col <- 8
if (!show_dens_only) {
if (pdf_or_cdf == "pdf") {
graphics::hist(y, col = my_col, probability = TRUE, las = 1,
axes = FALSE, xlab = my_xlab, ylab = my_ylab, main = "",
xpd = TRUE, xlim = my_xlim, ylim = c(0, ytop))
} else {
ecdfy <- stats::ecdf(y)
graphics::plot(ecdfy, col = my_col, las = 1, main = "",
axes = FALSE, xlab = my_xlab, ylab = my_ylab,
xpd = TRUE, xlim = my_xlim, ylim = c(0, ytop),
col.01line = 0)
}
if (show_dens) {
graphics::lines(x, ygev, xpd = TRUE, lwd = 2, lty = 2)
graphics::lines(x, ytrue, xpd = TRUE, lwd = 2, lty = 2, col = "red")
}
} else {
graphics::matplot(x, cbind(ygev, ytrue), col = c("black", "red"),
lwd = 2, lty = 2, ylab = my_ylab, las = 1,
xlab = my_xlab, xlim = my_xlim, ylim = c(0, ytop),
axes = FALSE, type = "l")
}
graphics::axis(2)
graphics::axis(1, line = 0.5)
if (show_dens_only) {
graphics::title(paste(the_distn, ", n = ", n))
}
if (!show_dens_only) {
if (show_bottom_rug) {
graphics::rug(y, line = 0.5, ticksize = 0.05)
}
graphics::rug(last_y, line = 0.5, ticksize = 0.05, col = "red", lwd = 2)
}
u_b <- my_xlim
my_leg_2 <- paste("GEV (", signif(gev_pars$loc, 2), ",",
signif(gev_pars$scale, 2), ",",
signif(gev_pars$shape, 2), ")" )
if (pdf_or_cdf == "pdf") {
my_leg_true <- expression(n * F ^ {n-1} * f)
if (show_dens || show_dens_only) {
graphics::legend(bottom_leg_pos, legend = c(my_leg_2, my_leg_true),
col = 1:2, lwd = 2, lty = 2, box.lty = 0,
cex = leg_cex)
}
} else {
my_leg_true <- expression(F ^ n)
if (!show_dens_only) {
if (show_dens) {
graphics::legend("topleft",
legend = c(my_leg_2, my_leg_true, "empirical cdf"),
col = c(1:2, 8), lwd = 2, lty = c(2, 2, -1),
pch = c(-1, -1, 16), box.lty = 0, cex = leg_cex)
} else {
graphics::legend("topleft",
legend = c(my_leg_2, my_leg_true, "empirical cdf"),
col = c(0, 0, 8), lwd = 2, lty = c(2, 2, -1),
pch = c(-1, -1, 16), box.lty = 0,
text.col = c(0, 0, 1), cex = leg_cex)
}
} else {
graphics::legend("topleft",
legend = c(my_leg_2, my_leg_true),
col = 1:2, lwd = 2, lty = 2, box.lty = 0,
cex = leg_cex)
}
}
if (!show_dens_only) {
top_ratio <- (last_y - u_t[1]) / (u_t[2] - u_t[1])
top_loc <- u_b[1] + (u_b[2] - u_b[1]) * top_ratio
if (arrow) {
graphics::segments(top_loc, ytop * 2, top_loc, ytop, col = "red",
xpd = TRUE, lwd = 2, lty = 2)
graphics::arrows(top_loc, ytop, last_y, 0, col = "red", lwd = 2, lty = 2,
xpd = TRUE, code = 2)
}
}
old_n <- n
old_pdf_or_cdf <- pdf_or_cdf
old_show_dens <- show_dens
old_show_dens_only <- show_dens_only
})
return(panel)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.