distributions3 (version 0.3.0)

distribution_calculate_moments: Workhorse function for numerically calculating central moments of probability distributions

Description

Method used to evaluate (approximate) the central moments (mean, variance, skewness, and kurtosis) for probability distributions for which only the cumulative distribution function (CDF) and - potentially - the quantile function is provided.

Usage

distribution_calculate_moments(
  x,
  what,
  gridsize = 500L,
  batchsize = 10000L,
  applyfun = NULL,
  cores = NULL,
  method = NULL,
  ...
)

Value

A (potentially named) numeric vector of length length(x) with the requested central moment.

Arguments

x

Object of class 'distribution'

what

integer. Controls what the C code returns: 1L (mean) 2L (variance), 3L (skewness), or 4L (kurtosis).

gridsize

positive integer. Size of the grid used to numerically calculate central moments.

batchsize

integer. Maximum batch size. Used to split the input into batches. Lower values reduce required memory but may increase computation time.

applyfun

An optional lapply-style function with arguments function(X, FUN, ...). It is used to compute the moment for each distribution in x. 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.

method

NULL (default) or character. Specifies how the grid for the calculation is set up. If NULL it is set to method = "cdf" if the distribution (y) is discrete and gridsize is large enough to span the required range for calculation. Else (including continuous distributions) method = "quantile" is used.

...

Additional arguments forwarded to distribution_calculate_moments().

Details

For discrete distributions spanning a range less than 0:gridsize the PDF is calculated at \(x = \{0, 1, 2, 3, \dots\}\) by differentiating the CDF provided, which is then used to calculate the central moments.

For continuous distributions as well as discrete distributions spanning a wide range of values (larger than 0:gridsize) a numeric grid with gridsize intervals is created. Given the distribution provides a quantile function, this grid is specified on a (mostly) uniform grid on the quantile scale. If no quantile function is provided, the 0.01 and 99.99 percentile are calculated via the CDF to serve as upper/lower bound for a uniform numeric grid. For each interval the density is approximated using numeric forward differences

$$f(x_j) = (F(x_{i+1}) - F(x_i)) / (x_{i+1} - x_i)$$

at each \(x_j = (x_{i+1} + x_i) \cdot 0.5\) with interval width of \(x_{i+1} - x_i\). The densities \(f(x_j)\) and interval mids \(x_j\) are used to calculate the weighted moments.

Examples

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

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

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

## Creating 3 distributions
dn <- MyNormal(mu = 1:3, sigma = 1:3 / 3)
dn <- setNames(dn, LETTERS[1:3])

## Calculating central moments
cbind(mean     = distribution_calculate_moments(dn, 1L),
      variance = distribution_calculate_moments(dn, 2L),
      skewness = distribution_calculate_moments(dn, 3L),
      kurtosis = distribution_calculate_moments(dn, 4L))

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

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

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

## Creating 3 distributions
dp <- MyPoisson(lambda = 1:3)
dp <- setNames(dp, LETTERS[4:6])

## Calculating central moments
cbind(mean     = distribution_calculate_moments(dp, 1L),
      variance = distribution_calculate_moments(dp, 2L),
      skewness = distribution_calculate_moments(dp, 3L),
      kurtosis = distribution_calculate_moments(dp, 4L))