R/quantiledensity.R

In spatstat.univar: One-Dimensional Probability Distribution Support for the 'spatstat' Family

#'
#'   quantiledensity.R
#'
#'  quantile method for class 'density'
#'
#'  Also a CDF from a 'density'
#' 
#'  $Revision: 1.5 $ $Date: 2025/04/01 02:43:51 $

quantile.density <- local({

  quantile.density <- function(x, probs = seq(0, 1, 0.25), names = TRUE, ...,
                               warn=TRUE) {
    stopifnot(inherits(x, "density"))
    #' check whether density estimate was restricted to an interval
    if(warn && is.call(cl <- x$call) && any(c("from", "to") %in% names(cl)))
      warning(paste("Density was normalised within the computed range",
                    "of x values", prange(c(cl$from, cl$to))),
              call.=FALSE)
    #' validate probs
    eps <- 100 * .Machine$double.eps
    if(any((p.ok <- !is.na(probs)) & (probs < -eps | probs > 1 + eps))) 
      stop("'probs' outside [0,1]")
    if (na.p <- any(!p.ok)) {
      o.pr <- probs
      probs <- probs[p.ok]
      probs <- pmax(0, pmin(1, probs))
    }
    np <- length(probs)
    qs <- rep(NA_real_, np)
    if (np > 0) {
      #' extract density values 
      xx <- x$x
      yy <- x$y
      nn <- length(xx)
      #' integrate, normalise
      Fx <- cumsum(yy * c(0, diff(xx)))
      Fx <- Fx/Fx[nn]
      #' quantile
      for(j in 1:np) {
        ii <- min(which(Fx >= probs[j]))
        if(!is.na(ii) && ii >= 1 && ii <= nn) 
          qs[j] <- xx[ii]
      }
      if (names && np > 0L) {
        names(qs) <- format_perc(probs)
      }
    }
    if (na.p) {
      o.pr[p.ok] <- qs
      if(names) {
        names(o.pr) <- rep("", length(o.pr))
        names(o.pr)[p.ok] <- names(qs)
      }
      return(o.pr)
    } else return(qs)
  }

  format_perc <- function (x, digits = max(2L, getOption("digits")),
                           probability = TRUE, use.fC = length(x) < 100, ...) {
    if (length(x)) {
      if (probability) x <- 100 * x
      paste0(if (use.fC) 
             formatC(x, format = "fg", width = 1, digits = digits)
      else format(x, trim = TRUE, digits = digits, ...), "%")
    }
    else character(0)
  }

  quantile.density
})


CDF <- function(f, ...) {
  UseMethod("CDF")
}

CDF.density <- function(f, ..., warn=TRUE) {
  stopifnot(inherits(f, "density"))
  #' check whether density estimate was restricted to an interval
  if(warn && is.call(cl <- f$call) && any(c("from", "to") %in% names(cl)))
    warning(paste("Density was normalised within the computed range",
                  "of x values", prange(c(cl$from, cl$to))),
            call.=FALSE)
  #' integrate
  xx <- f$x
  yy <- f$y
  nn <- length(xx)
  Fx <- cumsum(yy * c(0, diff(xx)))
  #' normalise
  Fx <- Fx/Fx[nn]
  #' 
  FF <- approxfun(xx, Fx, method="linear", rule=2)
  class(FF) <- c("interpolatedCDF", class(FF))
  return(FF)
}

print.interpolatedCDF <- function(x, ...)  {
  splat("Smooth estimate of cumulative distribution function")
  xrange <- range(get("x", envir=environment(x)))
  splat("Domain of x values:", prange(signif(xrange, 6)))
  NextMethod("print")
}

plot.interpolatedCDF <- function(x, ...)  {
  xname <- short.deparse(substitute(x))
  xrange <- range(get("x", envir=environment(x)))
  do.call(plot.function,
          resolve.defaults(x=quote(x),
                           list(...),
                           list(xlim=xrange,
                                xlab=expression(x),
                                ylab=expression(F(x)),
                                main=xname)))
}