pdf.distribution: Methods for numerically approximating distribution functions

Description

S3 methods for distribution objects used if certain (analytic) distribution functions are not available. Allows to use all distributions3 methods even if there are no dedicated methods for a certain distribution. See section 'Details' for more information.

Usage

# S3 method for distribution
pdf(
  d,
  x,
  drop = TRUE,
  elementwise = NULL,
  log = FALSE,
  applyfun = NULL,
  cores = NULL,
  ...
)
# S3 method for distribution
log_pdf(d, x, ...)
# S3 method for distribution
quantile(
  x,
  probs,
  drop = TRUE,
  elementwise = NULL,
  lower = -1/sqrt(.Machine$double.eps),
  upper = +1/sqrt(.Machine$double.eps),
  tol = .Machine$double.eps^0.5,
  maxit = 1000,
  ...
)
# S3 method for distribution
cdf(d, x, drop = TRUE, elementwise = NULL, lower.tail = TRUE, ...)

Arguments

d: An object of class distribution.
x: Either a numeric vector of probabilities to be evaluated (if pdf() is called), or an object of class distributions when calling the quantile() function.
drop: logical. Should the result be simplified to a vector if possible?
elementwise: logical. Should each distribution (in d/x) be evaluated at all elements in x (when pdf() is called) or probs (if quantile() is called)? By default (if elementwise = NULL) it is set to elementwise = TRUE if the lengths match, else elementwise is set FALSE.
log: logical. If TRUE, probabilities are given as log(p).
applyfun: An optional lapply-style function with arguments function(X, FUN, ...). It is used to compute the CRPS for each element of y. The default is to use the basic lapply function unless the cores argument is specified (see below).
cores: NULL or positive integer. If set to an integer the applyfun is set to mclapply with the desired number of cores, except on Windows where parLapply with makeCluster(cores) is used.
...: Currently ignored.
probs: Numeric vector of probabilities with values in [0,1].
lower, upper: numeric. Lower and upper end points for the interval to be searched, forwarded to uniroot.
tol: numeric. Desired accuracy for uniroot.
maxit: integer. Maximum number of iterations used when iteratively evaluating quantiles based on a pdf (discrete distributions only). If maxit is reached before the quantile has been found, an error will be thrown.
lower.tail: logical. If TRUE (default), probabilities are \(P[X \le x]\), else \(P[X \ge x]\).

Details

For distribution classes that do not (or not yet) provide all S3 methods for the usual generic functions (pdf, cdf, quantile, random), a fallback method for the general distribution class is provided. This fallback method numerically approximates the probability density function (PDF), cumulative distribution function (CDF), quantile function, or random number generation, provided that at least some of the methods are available (see below).

S3 methods for is_discrete and support are mandatory. Besides these two (fairly simple) methods, the following methods are required.

Continuous distributions:

provide a cdf.* method,
or a pdf.* and quantile.* method.

Discrete distributions:

provide a pdf.* method,
or a cdf.* method.

Any additional S3 method provided on top of the minimal requirements will always be leveraged if available (i.e., the more dedicated/analytical methods provided the better).

Examples

Run this code

## ------------- custom Normal distribution (MyNormal1) ---------------
## -------------------------- using cdf -------------------------------

## Constructor function for new 'MyNormal1' distribution
MyNormal1 <- function(mu, sigma) {
    d <- data.frame(mu = mu, sigma = sigma)
    class(d) <- c("MyNormal1", "distribution")
    return(d)
}

## Additionally required S3 methods (borrowed from class Normal)
registerS3method("cdf", "MyNormal1",
                 getS3method("cdf", class = "Normal"))
registerS3method("is_discrete", "MyNormal1",
                 getS3method("is_discrete", class = "Normal"))
registerS3method("support", "MyNormal1",
                 getS3method("support", class = "Normal"))

## Constructing objects; three normal distributions with
## mean c(1, 2, 3) and standard deviation c(1, 2.5, 5).
## mn3: Based on MyNormal1 where only cdf, id_discrete, and support
##      are defined.
## n3:  Analytic solution (Normal()) with analytic functions
##      for all distribution functions (pdf, cdf, quantile) as well as for
##      the first four central moments (mean, variance, skewness, kurtosis)
mn3 <- MyNormal1(mu = 1:3, sigma = c(1, 2.5, 5))
n3  <- Normal(mu = 1:3, sigma = c(1, 2.5, 5))

## Class 'MyNormal1' knows the analytic cdf:
cdf(mn3, x = 2)
identical(cdf(mn3, x = 2), cdf(n3, x = 2))

## Calculating probability at x = 2
pdf(mn3, x = 2) ## Numeric approximation
pdf(n3,  x = 2) ## Analytic solution
pdf(mn3, x = 2) - pdf(n3, x = 2) ## Pairwise differences/precision

## Calculating quantiles
probs <- c(0.0, 0.01, 0.25, 0.5, 0.75, 0.99, 1.0)
quantile(mn3, probs = probs) ## Numeric approximation
quantile(n3,  probs = probs) ## Analytic solution

probs2 <- seq(0.01, 0.99, by = 0.01)
qmn3   <- quantile(mn3, probs = probs2) ## Numeric approximation
qn3    <- quantile(n3,  probs = probs2) ## Analytic solution
range(qmn3 - qn3) ## Range of pairwise differences/precision

## Central moments
cbind(mean     = mean(mn3),     variance = variance(mn3),
      skewness = skewness(mn3), kurtosis = kurtosis(mn3))

## Visual comparison: density
x    <- seq(-3, 5, by = 0.1)
pmn3 <- pdf(mn3[1], x = x)
pn3  <- pdf(n3[1],  x = x)
dpdf <- data.frame(x          = x,
                   analytical = pdf(n3[1],  x = x),
                   numerical  = pdf(mn3[1], x = x))

matplot(dpdf[, 1], dpdf[, -1],  type = "l", col = 1:2,
        xlab = NA, ylab = "density", lty = 1, lwd = 2:1,
        main = "Density function (mean = 1, sigma = 1)")
legend("bottom", legend = names(dpdf)[-1], bty = "n",
       col = 1:2, lty = 1, lwd = 2:1)

## Visual comparison: quantiles
probs     <- c(0.001, seq(0.01, 0.99, by = 0.01), 0.999)
dquantile <- data.frame(probs      = probs,
                        analytical = quantile(n3[1],  probs = probs),
                        numerical  = quantile(mn3[1], probs = probs))

matplot(dquantile[, -1], dquantile[, 1], type = "l", col = 1:2,
        xlab = NA, ylab = "probability", lty = 1, lwd = 2:1,
        main = "Quantile function (mean = 1, sigma = 1)")
legend("bottom", legend = names(dquantile)[-1], bty = "n",
       col = 1:2, lwd = 2:1)

## Visual comparison: quantile-quantile plot
plot(dquantile[, "analytical"], dquantile[, "analytical"],
     xlab = "numerically approximated quantiles",
     ylab = "theoretical quantiles",
     main = "QQ-plot Normal vs. MyNormal1")
abline(0, 1, col = 2, lty = 2)

## Drawing random numbers (500 on third distribution)
set.seed(6020); rmn <- random(mn3[3], 500L)
set.seed(6020); rn  <- random(n3[3],  500L)

hmn <- hist(rmn, breaks = 15L, plot = FALSE)
hn  <- hist(rn,  breaks = 15L, plot = FALSE)

plot(hmn$mids, hmn$density, type = "l",
     xlab = NA, ylab = "density",
     main = "Density of random numbers: Normal vs. MyNormal1")
lines(hn$mids, hn$density, col = 2)

## ------------ custom Normal distribution (MyNormal2) ----------------
## ------------------- using pdf and quantile -------------------------

## Constructor function for new 'MyNormal2' distribution
MyNormal2 <- function(mu, sigma) {
    d <- data.frame(mu = mu, sigma = sigma)
    class(d) <- c("MyNormal2", "distribution")
    return(d)
}

## Additionally required S3 methods (borrowed from class Normal)
registerS3method("pdf", "MyNormal2",
                 getS3method("pdf", class = "Normal"))
registerS3method("quantile", "MyNormal2",
                 getS3method("quantile", class = "Normal"))
registerS3method("is_discrete", "MyNormal2",
                 getS3method("is_discrete", class = "Normal"))
registerS3method("support", "MyNormal2",
                 getS3method("support", class = "Normal"))

## Constructing objects; creating three (named) 'MyNormal2' distributions
## with mean c(1, 2, 3) and standard deviation c(1, 2.5, 5) for which only
## pdf, quantile, is_discrete, and support are defined.
mn3 <- MyNormal2(mu = 1:3, sigma = c(1, 2.5, 5))
mn3 <- setNames(mn3, LETTERS[1:3])

random(mn3, n = 3L)
cdf(mn3, x = 2)
pdf(mn3, x = 2)
quantile(mn3, 0.5)
cbind(mean     = mean(mn3),     variance = variance(mn3),
      skewness = skewness(mn3), kurtosis = kurtosis(mn3))

# \donttest{

## ------------ custom Poisson distribution (MyPoisson1) --------------
## -------------------------- using pdf -------------------------------

## Custom constructor function for the 'MyPoisson1' distribution
MyPoisson1 <- function(lambda) {
    d <- data.frame(lambda = lambda)
    class(d) <- c("MyPoisson1", "distribution")
    return(d)
}

## Additionally required S3 methods (borrowed from class Poisson)
registerS3method("pdf", "MyPoisson1",
                 getS3method("pdf", class = "Poisson"))
registerS3method("is_discrete", "MyPoisson1",
                 getS3method("is_discrete", class = "Poisson"))
registerS3method("support", "MyPoisson1",
                 getS3method("support", class = "Poisson"))

## Constructing objects; three normal distributions with
## parameter lambda = c(1, 2.5, 5).
## mp3: Based on MyPoisson1 where only pdf, id_discrete, and support
##      are defined.
## p3:  Analytic solution (Poisson()) with analytic
##      functions for all distribution functions (pdf, cdf, quantile)
##      as well as for the first four central moments (mean,
##      variance, skewness, kurtosis)
mp3 <- MyPoisson1(lambda = c(1, 2.5, 5))
p3  <- Poisson(lambda = c(1, 2.5, 5))

## Class 'MyPoisson1' knows the analytic cdf:
pdf(mp3, x = 2)
identical(pdf(mp3, x = 2), pdf(p3, x = 2))

## Calculating distribution at x = 2
cdf(mp3, x = 2) ## Numeric approximation
cdf(mp3, x = 2) ## Analytic solution
cdf(mp3, x = 2) - cdf(p3, x = 2) ## Pairwise differences/precision

## Calculating quantiles
probs <- c(0.0, 0.01, 0.25, 0.5, 0.75, 0.99, 1.0)
quantile(mp3, probs = probs) ## Numeric approximation
quantile(p3,  probs = probs) ## Analytic solution

probs2 <- seq(0.01, 0.99, by = 0.01)
qmp3   <- quantile(mp3, probs = probs2) ## Numeric approximation
qp3    <- quantile(p3,  probs = probs2) ## Analytic solution
range(qmp3 - qp3) ## Range of pairwise differences/precision

## Central moments
cbind(mean     = mean(mp3),     variance = variance(mp3),
      skewness = skewness(mp3), kurtosis = kurtosis(mp3))

## Visual comparison: distribution function
x    <- seq(-5, 20, by = 1)
dcdf <- data.frame(x          = x,
                   analytical = cdf(p3[2L], x = x),
                   numerical  = cdf(mp3[2L], x = x))

matplot(dcdf[, 1], dcdf[, -1], type = "s", col = 1:2,
        xlab = NA, ylab = "probability",
        lty = 1, lwd = 2:1, main = "Distribution function (lambda = 2.5)")
legend("bottomright", legend = names(dcdf)[-1], lty = 1,
       col = 1:2, lwd = 2:1, bty = "n")


## Visual comparison: quantile function
probs     <- seq(0.01, 0.99, by = 0.01)
dquantile <- data.frame(probs = probs,
                        analytical = quantile(p3[2L],  probs = probs),
                        numerical  = quantile(mp3[2L], probs = probs))

matplot(dquantile[, -1], dquantile[, 1], type = "s", col = 1:2,
        xlab = NA, ylab = "probability",
        lty = 1, lwd = 2:1, main = "Quantile function (lambda = 2.5)")
legend("bottomright", legend = names(dquantile)[-1], lty = 1,
       col = 1:2, lwd = 2:1, bty = "n")

## Quantile-Quantile plot
probs <- seq(0.01, 0.99, by = 0.01)
plot(quantile(p3[2], probs), quantile(mp3[2], probs),
     main = "QQ-plot Poisson vs. MyPoisson1",
     xlab = "approximated quantiles", ylab = "theoretical quantiles")
abline(0, 1, col = 2, lty = 2)

## Drawing random numbers (200 on third distribution)
set.seed(6020); rmp <- random(mp3[3], 100L)
set.seed(6020); rp  <- random(p3[3],  100L)

hmp <- hist(rmp, breaks = 0:15, plot = FALSE)
hp  <- hist(rp,  breaks = 0:15, plot = FALSE)

plot(hmp$breaks[-1], hmp$density, type = "s",
     xlab = NA, ylab = "density",
     main = "Density of random numbers: Poisson vs. MyPoisson1")
lines(hp$breaks[-1], hp$density, col = 2, type = "s")


## ------------ custom Poisson distribution (MyPoisson2) --------------
## -------------------------- using cdf -------------------------------

## Custom constructor function for the 'MyPoisson2' distribution
MyPoisson2 <- function(lambda) {
    d <- data.frame(lambda = lambda)
    class(d) <- c("MyPoisson2", "distribution")
    return(d)
}

## Additionally required S3 methods (borrowed from class Poisson)
registerS3method("cdf", "MyPoisson2",
                 getS3method("cdf", class = "Poisson"))
registerS3method("is_discrete", "MyPoisson2",
                 getS3method("is_discrete", class = "Poisson"))
registerS3method("support", "MyPoisson2",
                 getS3method("support", class = "Poisson"))


## Constructing objects; creating three (named) 'MyPoisson2' distributions
## with lambda c(1, 2.5, 5) for which only cdf, is_discrete,
## and support are defined.
mp3 <- MyPoisson2(lambda = c(1, 2.5, 5))
mp3 <- setNames(mp3, LETTERS[4:6])

random(mp3, n = 3L)
cdf(p3, x = 2)
pdf(p3, x = 2)
quantile(p3, 0.5)
cbind(mean     = mean(p3),     variance = variance(p3),
      skewness = skewness(p3), kurtosis = kurtosis(p3))

# }

Run the code above in your browser using DataLab