R/diagnostics.R

In strakaps/CTRE: Thresholding Bursty Time Series

#' Mittag-Leffler QQ Plot
#'
#' Generates a QQ plot for assessing the fit of a Mittag-Leffler
#' distribution.
#'
#' @param x
#'     A vector of data to be compared against the Mittag-Leffler
#'     distribution.
#' @param tail
#'    Tail parameter of the Mittag-Leffler population. Default is
#'    \code{1}, i.e. the exponential distribution.
#' @param scale
#'    Scale parameter of the Mittag-Leffler population, if known.
#' @param ...
#'     Additional plotting arguments, e.g. \code{log = 'xy'}.
#' @examples
#'   library(magrittr)
#'   flares %>% ctre() %>% thin(k=200) %>% interarrival() %>% mlqqplot(tail = 1, log = 'xy')
#'   flares %>% ctre() %>% thin(k=200) %>% interarrival() %>% mlqqplot(tail = 0.8, log = 'xy')
#'
#'   seaquakes %>% ctre() %>% thin(k=150) %>% interarrival() %>% mlqqplot(tail = 0.9, log = 'xy')
#' @export
mlqqplot <- function(x,
                     tail = 1,
                     scale = NULL,
                     ...) {
  n <- length(x)
  if (is.null(scale))
    scale <- 1
  pop <-
    MittagLeffleR::qml(p = stats::ppoints(n),
                       tail = tail,
                       scale = scale)
  graphics::plot(pop,
                 sort(x),
                 xlab = "Population Quantiles",
                 ylab = "Sample Quantiles",
                 ...)
}

#' Static QQ Plot estimator
#'
#' Generates a static QQ plot for Pareto tail estimates
#'
#' @param data
#'     A vector of data.
#' @param top_k
#'     Only use the top_k largest values (in the tail) for the plot
#' @param plot_me
#'     Should a plot be produced? If not, only the estimate is returned.
#' @param ...
#'     Additional plotting arguments
#' @return
#'     An estimate of the Pareto tail exponent (invisible if plotted).
#' @export

qqestplot_static <- function(data, top_k = NULL, plot_me = TRUE, ...) {
  n <- length(data)
  if (is.null(top_k))
    top_k <- n
  else if (top_k > n)
    stop("Can't choose top k order statistics, only have", n)
  x <- -log(sort(stats::ppoints(top_k)))
  y <- log(sort(data, decreasing = TRUE)[1:top_k])
  l <- stats::lm(y~x)
  if (plot_me) {
    graphics::plot(x, y, ...)
    graphics::abline(l, col = 2)
    return (1 / l$coefficients[2])
  } else
    invisible(1 / l$coefficients[2])
}



#' Autocorrelation function
#'
#' @param x time series or ctre object.
#' @param ... Additional arguments passed to \code{stats::\link[stats]{acf}}
#' @seealso \code{\link{acf.ctre}}
#' @export
acf <- function(x, ...)
  UseMethod("acf", x)


#' @export
acf.default <- function(x, ...) stats::acf(x)

#' Autocorrelation function
#'
#' Calculates and plots the autocorrelation function for the bivariate
#' time series of interarrival times and magnitudes.
#'
#' @param x An object of class \code{\link[CTRE]{ctre}}
#' @param OCTRE
#'     If FALSE (default), each magnitude is matched with its preceding
#'     interarrival time. If TRUE, each magnitude is matched with its
#'     succeeding interarrival time.
#' @param ...
#'     Additional arguments passed to \code{\link[stats]{acf}}
#' @examples
#'   library(magrittr)
#'   flares %>% ctre() %>% thin(k=150) %>% acf()
#' @export
acf.ctre <- function(x, OCTRE = FALSE, ...){
  T_ell <- interarrival(x)
  X_ell <- magnitudes(x)
  n <- length(x)
  assertthat::are_equal(length(T_ell), n)
  if (OCTRE)
    T_ell <- T_ell[-1]
  else
    X_ell <- X_ell[-n]
  acf(cbind(T_ell, X_ell), ...)
}

#' Plot empirical copula
#'
#' Plots the ranks of the magnitudes against the ranks of the preceding
#' (or succeeding) interarrival times.
#'
#' @param ctre A \code{\link{ctre}} object
#' @param OCTRE
#'     Shall each magnitude be matched with the preceding interarrival
#'     time (FALSE) or the succeeding interarrival time (TRUE)?
#' @param ... Additional plotting arguments
#' @examples
#'   library(magrittr)
#'   flares %>% ctre() %>% thin(k = 300) %>% empcopula(pch = '*')
#' @export
empcopula <- function(ctre, OCTRE = FALSE, ...){
  T_ell <- interarrival(ctre)
  X_ell <- magnitudes(ctre)
  n <- length(ctre)
  assertthat::are_equal(length(X_ell), n)
  if (OCTRE)
    X_ell <- X_ell[-1]
  else
    X_ell <- X_ell[-n]
  TT <- rank(T_ell)/n
  JJ <- rank(X_ell)/n
  ggplot2::ggplot(data.frame(TT, JJ),
                  ggplot2::aes(x=TT, y=JJ)) +
    ggplot2::geom_point() +
    ggplot2::labs(ggplot2::ggtitle("Emp. Copula: Exc & Exc Time")) +
    ggplot2::xlab("rank(Exc. Time)") +
    ggplot2::ylab("rank(Exc.)")
  #graphics::plot(rank(T_ell)/n, rank(X_ell)/n, main = "Emp. Copula (Exc & Exc Time)", ...)
}