R/pmfcm.R

In eFCM: Exponential Factor Copula Model

#' CDF of the exponential Factor Copula Model (vector input)
#'
#' Computes the eFCM-based \eqn{P(W \leq w)} for a single \eqn{d}-dimensional vector \eqn{w}.
#'
#' @param w Numeric vector of length \eqn{d}.
#' @param lambda,delta Positive scalars: common-factor rate \eqn{\lambda} and range \eqn{\delta}.
#' @param dist Optional \eqn{d\times d} distance matrix. If \code{NULL}, provide \code{coord}.
#' @param coord Optional two-column matrix/data.frame of coordinates (lon, lat) to build \code{dist}.
#' @param smooth Matérn smoothness \eqn{\nu} (default \code{0.5}).
#' @param abseps,releps Absolute/relative tolerances for the MVN CDF.
#' @param maxpts Maximum number of function evaluations for the MVN CDF.
#' @param miles Logical; passed to \code{fields::rdist.earth()} if \code{coord} is used.
#'
#' @return A single numeric CDF value in \eqn{[0,1]}.
#' @examples
#' data(LonLat)
#' d <- 2
#' w <- rep(0.3, d)
#' pmfcm(w, lambda = 2, delta = 100, coord = LonLat[1:2, ])
#' @importFrom fields rdist.earth
#' @export
pmfcm <- function(w, lambda, delta, dist = NULL, coord = NULL,
                  smooth = 0.5, abseps = 1e-5, releps = 1e-5, maxpts = 25000,
                  miles = FALSE) {
  stopifnot(is.numeric(w), is.vector(w), all(is.finite(w)))
  stopifnot(is.numeric(lambda), length(lambda) == 1, lambda > 0,
            is.numeric(delta), length(delta) == 1, delta > 0,
            is.numeric(smooth), length(smooth) == 1, smooth > 0)

  if (is.null(dist)) {
    if (is.null(coord)) stop("Provide either `dist` or `coord`.", call. = FALSE)
    if (!requireNamespace("fields", quietly = TRUE))
      stop("Package 'fields' is required when `coord` is supplied.", call. = FALSE)
    dist <- fields::rdist.earth(as.matrix(coord), miles = miles)
  } else {
    dist <- as.matrix(dist)
  }

  if (nrow(dist) != ncol(dist))
    stop("`dist` must be a square matrix.", call. = FALSE)
  if (length(w) != nrow(dist))
    stop("length(w) must equal nrow(dist).", call. = FALSE)

  pmfcm_rcpp(
    w        = as.numeric(w),
    lambda     = lambda,
    delta    = delta,
    dist     = dist,
    smooth   = smooth,
    abseps   = abseps,
    releps   = releps,
    maxpts   = as.integer(maxpts)
  )
}