Nothing
#' Calculate fire intervals from a `composite`
#'
#' @param comp A `composite` instance, usually output from [composite()].
#' Should contain only one series.
#' @param densfun String giving desired distribution to fit. Either "weibull"
#' or "lognormal". Default is "weibull".
#'
#' @return An `intervals` object. `intervals` have components:
#' * "intervals" an integer vector giving the actual fire intervals.
#' * "fitdistr" a `fitdistr` object from [MASS::fitdistr()] representing the
#' density function fit.
#' * "densfun" a string giving the name of the density function used.
#' * "kstest" an `htest` object from [stats::ks.test()] giving the result of a
#' one-sample Kolmogorov-Smirnov test.
#' * "shapirotest" an `htest` object from [stats::shapiro.test()] giving the
#' result of a Shapiro-Wilk normality test.
#' * "comp_name" a string giving the name of the interval's input composite.
#' * "event_range" an integer vector giving the year range (min, max) of
#' events used to create this intervals.
#'
#' @seealso
#' * [composite()] to create a `composite` object.
#' * [mean.intervals()] gets mean fire interval.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [plot_intervals_dist()] plots `intervals`.
#' * [min.intervals()] gives the minimum fire interval.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @examples
#' data(pgm)
#' interv <- intervals(composite(pgm))
#' print(interv)
#'
#' mean(interv) # Mean interval
#'
#' # Now fit log-normal distribution instead of Weibull.
#' intervals(composite(pgm), densfun = "lognormal")
#' \dontrun{
#' # Boxplot of fire interval distribution.
#' boxplot(intervals(composite(pgm))$intervals)
#' }
#'
#' @export
intervals <- function(comp, densfun = "weibull") {
stopifnot(is_fhx(comp))
stopifnot(densfun %in% c("weibull", "lognormal"))
if (length(series_names(comp)) > 1) {
stop("Found multiple series in `comp`. There can only be one.")
}
dens2cum <- list(weibull = stats::pweibull, lognormal = stats::plnorm)
x <- list()
x$intervals <- diff(get_event_years(comp)[[1]])
if (length(x$intervals) < 2) {
stop("Too few fire events to compute intervals")
}
class(x) <- c("intervals")
x$fitdistr <- MASS::fitdistr(x$intervals, densfun)
x$densfun <- densfun
p_densfun <- dens2cum[[densfun]]
kstest_args <- c(
list(x = x$intervals, y = p_densfun),
x$fitdistr$estimate
)
x$kstest <- do.call(stats::ks.test, kstest_args)
x$shapirotest <- stats::shapiro.test(x$intervals)
x$comp_name <- series_names(comp)
x$event_range <- range(unlist(get_event_years(comp)))
x
}
#' Check if object is fire `intervals`
#'
#' @param x An R object.
#'
#' @return Boolean indicating whether x is an `intervals` object.
#'
#' @seealso [intervals()] creates an `intervals` object.
#' @export
is_intervals <- function(x) inherits(x, "intervals")
#' Alias to [is_intervals()]
#'
#' @inherit is_intervals
#'
#' @export
is.intervals <- function(x) is_intervals(x)
#' Fire `intervals` arithmetic mean
#'
#' @param x An `intervals` object.
#' @param ... Additional arguments passed to [mean()].
#'
#' @return Numeric or NA.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [min.intervals()] gives the minimum fire interval.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @export
mean.intervals <- function(x, ...) {
mean(x$intervals, ...)
}
#' Fire `intervals` median
#'
#' @param x An `intervals` object.
#' @param ... Additional arguments passed to [stats::median()].
#'
#' @return Numeric or NA.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [mean.intervals()] gets mean fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [min.intervals()] gives the minimum fire interval.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @importFrom stats median
#' @export
median.intervals <- function(x, ...) {
median(x$intervals, ...)
}
#' Minimum interval in fire `intervals`
#'
#' @param x An `intervals` object.
#' @param ... Additional arguments passed to [min()].
#'
#' @return Numeric or NA.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [mean.intervals()] gets median fire interval.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @export
min.intervals <- function(x, ...) {
min(x$intervals)
}
#' Maximum interval in fire `intervals`
#'
#' @param x An `intervals` object.
#' @param ... Additional arguments passed to [max()].
#'
#' @return Numeric or NA.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [mean.intervals()] gets median fire interval.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [min.intervals()] gives the minimum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @export
max.intervals <- function(x, ...) {
max(x$intervals)
}
#' Fit distribution quantiles to fire `intervals`
#'
#' @param x An `intervals` object.
#' @param q Vector giving the desired quantiles.
#' @param ... Additional arguments passed to the [quantile()] method for the
#' fit distribution.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [mean.intervals()] gets median fire interval.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [min.intervals()] gives the minimum fire interval.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @examples
#' data(pgm)
#' intervs <- intervals(composite(pgm))
#' quantile(intervs)
#'
#' # Or you can pass in your own quantiles:
#' quantile(intervs, q = c(0.25, 0.5, 0.75))
#' @importFrom stats quantile
#' @export
quantile.intervals <- function(x, q = c(0.125, 0.5, 0.875), ...) {
dens2cum <- list(weibull = stats::qweibull, lognormal = stats::qlnorm)
q_densfun <- dens2cum[[x$densfun]]
quant_args <- c(
list(q, ...),
x$fitdistr$estimate
)
quants <- do.call(q_densfun, quant_args)
quants
}
#' Plot a fire `intervals` object
#'
#' @param ... Arguments passed to [plot_intervals_dist()].
#'
#' @seealso [plot_intervals_dist()] plot `intervals` distributions.
#'
#' @examples
#' data(pgm)
#' interv <- intervals(composite(pgm))
#'
#' plot(interv, binwidth = 5)
#' @export
plot.intervals <- function(...) {
print(plot_intervals_dist(...))
}
#' Basic fire `intervals` distribution plot
#'
#' @param x An `intervals` object, from [intervals()].
#' @inheritParams ggplot2::geom_histogram
#'
#' @return A ggplot object.
#'
#' @seealso
#' * [intervals()] to create a fire `intervals` object.
#' * [mean.intervals()] gets mean fire interval.
#' * [median.intervals()] gets median fire interval.
#' * [quantile.intervals()] get fit distribution quantiles.
#' * [min.intervals()] gives the minimum fire interval.
#' * [max.intervals()] gives the maximum fire interval.
#' * [print.intervals()] prints common fire-interval summary statistics.
#'
#' @export
plot_intervals_dist <- function(x, binwidth = NULL) {
p <- ggplot2::ggplot(
data.frame("intervals" = x[["intervals"]]),
ggplot2::aes_string(x = "intervals")
)
if (is.null(binwidth)) {
p <- (p + ggplot2::geom_histogram())
} else {
p <- (p + ggplot2::geom_histogram(binwidth = binwidth))
}
p <- (p + ggplot2::geom_rug() + ggplot2::theme_bw())
p
}
#' Print a fire `intervals` object
#'
#' @param x An `intervals` object.
#' @param ... Additional arguments that are tossed.
#'
#' @seealso [intervals()] to create a fire `intervals` object.
#'
#' @examples
#' data(pgm)
#' interv <- intervals(composite(pgm))
#' print(interv)
#'
#' # Note, you can also catch the printed table:
#' summary_stats <- print(interv)
#'
#' @export
print.intervals <- function(x, ...) {
wfit <- MASS::fitdistr(x$intervals, "weibull")
weib_invshape <- 1 / wfit$estimate["shape"]
weibull_median <- wfit$estimate["scale"] * log(2)^weib_invshape # nolint
quants <- quantile(x, q = c(0.847, 0.5, 0.125))
cat(strwrap("Interval Analysis", prefix = "\t"), sep = "\n")
cat(strwrap("=================", prefix = "\t"), sep = "\n")
cat("\n")
cat(paste0("Composite name: ", x$comp_name, "\n"))
cat(paste0(
"Events range: ", x$event_range[1],
" to ", x$event_range[2], "\n"
))
cat("\n")
cat(paste0("\tTotal intervals: ", length(x$intervals), "\n"))
cat(paste0("\tMean interval: ", round(mean(x), 1), "\n"))
cat(paste0("\tMedian interval: ", round(median(x), 1), "\n"))
cat(paste0("\tWeibull median: ", round(weibull_median, 1), "\n"))
cat(paste0("\tStandard deviation: ", round(stats::sd(x$intervals), 1), "\n"))
cat(paste0("\tMinimum interval: ", min(x), "\n"))
cat(paste0("\tMaximum interval: ", max(x), "\n"))
cat("\n\n")
cat(strwrap("Theoretical distribution"), sep = "\n")
cat(strwrap("------------------------"), sep = "\n")
cat("\n")
cat(paste0("Fit distribution: ", x$densfun, "\n\n"))
print(x$fitdistr)
cat("\n\n")
cat(strwrap("Percentiles"), sep = "\n")
cat(strwrap("-----------"), sep = "\n")
cat("\n")
cat(strwrap(paste0(
"lower (12.5%): ", round(quants[1], 1), " | ", x$densfun, " ",
"median (50.0%): ", round(quants[2], 1), " | ",
"upper (87.5%): ", round(quants[3], 1), "\n"
)))
cat("\n\n")
cat(strwrap(x$kstest$method), sep = "\n")
cat(strwrap(strrep("-", nchar(x$kstest$method))), sep = "\n")
cat("\n")
cat(paste0(
"D^- = ", round(x$kstest$statistic, 5),
", p = ", round(x$kstest$p.value, 5), "\n"
))
null_msg <- paste0(
"Null hypothesis: ",
"The intervals were sampled from the theoretical distribution."
)
cat(strwrap(null_msg, exdent = 4), sep = "\n")
alt_msg <- paste0(
"Alt. hypothesis: ",
"The intervals were not sampled from the theoretical distribution."
)
cat(strwrap(alt_msg, exdent = 4), sep = "\n")
cat("\n\n")
invisible(x)
}
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.