R Under development (unstable) (2023-11-07 r85491 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ### actuar: Actuarial Functions and Heavy Tailed Distributions
> ###
> ### Tests of functions for continuous and discrete probability
> ### distributions.
> ###
> ### Despite the name of the file, the tests are for [dpqrm,lev]<dist>
> ### functions (for continuous distributions):
> ###
> ###   d<dist>: density or probability mass
> ###   p<dist>: cumulative distribution
> ###   q<dist>: quantile
> ###   r<dist>: random number generation
> ###   m<dist>: moment
> ###   lev<dist>: limited moment
> ###
> ### Distributions are classified and sorted as in appendix A and
> ### appendix B of the 'distributions' package vignette.
> ###
> ### Inspired by (and some parts taken from) `tests/d-p-q-r-tests.R` in
> ### R sources.
> ###
> ### AUTHOR: Vincent Goulet <vincent.goulet@act.ulaval.ca> with
> ###         indirect help from the R Core Team
> 
> ## Load the package
> library(actuar)

Attaching package: 'actuar'

The following objects are masked from 'package:stats':

    sd, var

The following object is masked from 'package:grDevices':

    cm

> library(expint)                         # for gammainc
> 
> ## Define a "local" version of the otherwise non-exported function
> ## 'betaint'.
> betaint <- actuar:::betaint
> 
> ## No warnings, unless explicitly asserted via tools::assertWarning.
> options(warn = 2)
> assertWarning <- tools::assertWarning
> 
> ## Special values and utilities. Taken from `tests/d-p-q-r-tests.R`.
> Meps <- .Machine$double.eps
> xMax <- .Machine$double.xmax
> xMin <- .Machine$double.xmin
> All.eq <- function(x, y)
+ {
+     all.equal.numeric(x, y, tolerance = 64 * .Machine$double.eps,
+                       scale = max(0, mean(abs(x), na.rm = TRUE)))
+ }
> 
> if(!interactive())
+     set.seed(123)
> 
> ###
> ### CONTINUOUS DISTRIBUTIONS
> ###
> 
> ##
> ## FELLER-PARETO AND PARETO II, III, IV DISTRIBUTIONS
> ##
> 
> ## When reasonable, we also test consistency with the special cases
> ## min = 0:
> ##
> ##   Feller-Pareto -> Transformated beta
> ##   Pareto IV     -> Burr
> ##   Pareto III    -> Loglogistic
> ##   Pareto II     -> Pareto
> 
> ## Density: first check that functions return 0 when scale = Inf, and
> ## when x = scale = Inf.
> stopifnot(exprs = {
+     dfpareto(c(42, Inf), min = 1, shape1 = 2, shape2 = 3, shape3 = 4, scale = Inf) == c(0, 0)
+     dpareto4(c(42, Inf), min = 1, shape1 = 2, shape2 = 3, scale = Inf) == c(0, 0)
+     dpareto3(c(42, Inf), min = 1, shape  = 3, scale = Inf) == c(0, 0)
+     dpareto2(c(42, Inf), min = 1, shape  = 2, scale = Inf) == c(0, 0)
+ })
> 
> ## Next test density functions for an array of standard values.
> nshpar <- 3                     # (maximum) number of shape parameters
> min <- round(rnorm(30, 2), 2)
> shpar <- replicate(30, rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> for (i in seq_along(min))
+ {
+     m <- min[i]
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, t)
+     for (s in scpar)
+     {
+         x <- rfpareto(100, min = m, shape1 = a, shape2 = g, shape3 = t, scale = s)
+         y <- (x - m)/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             all.equal(d1 <- dfpareto(x, min = m,
+                                      shape1 = a, shape2 = g, shape3 = t,
+                                      scale = s),
+                       d2 <- dfpareto(y, min = 0,
+                                      shape1 = a, shape2 = g, shape3 = t,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       dtrbeta(y,
+                               shape1 = a, shape2 = g, shape3 = t,
+                               scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * y^(g*t - 1)/(s * Be * (1 + y^g)^(a + t)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * u^t * (1 - u)^a/((x - m) * Be),
+                       tolerance = 1e-10)
+         })
+         x <- rpareto4(100, min = m, shape1 = a, shape2 = g, scale = s)
+         y <- (x - m)/s
+         u <- 1/(1 + y^g)
+         stopifnot(exprs = {
+             all.equal(d1 <- dpareto4(x, min = m,
+                                      shape1 = a, shape2 = g,
+                                      scale = s),
+                       d2 <- dpareto4(y, min = 0,
+                                      shape1 = a, shape2 = g,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       dburr(y,
+                             shape1 = a, shape2 = g,
+                             scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * g * y^(g - 1)/(s * (1 + y^g)^(a + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * g * u^a * (1 - u)/(x - m),
+                       tolerance = 1e-10)
+         })
+         x <- rpareto3(100, min = m, shape = g, scale = s)
+         y <- (x - m)/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             all.equal(d1 <- dpareto3(x, min = m,
+                                      shape = g,
+                                      scale = s),
+                       d2 <- dpareto3(y, min = 0,
+                                      shape = g,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       dllogis(y,
+                               shape = g,
+                               scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * y^(g - 1)/(s * (1 + y^g)^2),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * u * (1 - u)/(x - m),
+                       tolerance = 1e-10)
+         })
+         x <- rpareto2(100, min = m, shape = a, scale = s)
+         y <- (x - m)/s
+         u <- 1/(1 + y)
+         stopifnot(exprs = {
+             all.equal(d1 <- dpareto2(x, min = m,
+                                      shape = a,
+                                      scale = s),
+                       d2 <- dpareto2(y, min = 0,
+                                      shape = a,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       dpareto(y,
+                               shape = a,
+                               scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a/(s * (1 + y)^(a + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * u^a * (1 - u)/(x - m),
+                       tolerance = 1e-10)
+         })
+     }
+ }
> 
> ## Tests on the cumulative distribution function.
> ##
> ## Note: when shape1 = shape3 = 1, the underlying beta distribution is
> ## a uniform. Therefore, pfpareto(x, min, 1, shape2, 1, scale) should
> ## return the value of u = v/(1 + v), v = ((x - min)/scale)^shape2.
> ##
> ## x = 2/Meps = 2^53 (with min = 0, shape2 = scale = 1) is the value
> ## where the cdf would jump to 1 if we weren't using the trick to
> ## compute the cdf with pbeta(1 - u, ..., lower = FALSE).
> scLrg <- 1e300 * c(0.5, 1, 2)
> m <- rnorm(1)
> stopifnot(exprs = {
+     pfpareto(Inf,         min = 10,   1, 2, 3, scale = xMax) == 1
+     pfpareto(2^53,        min = 0,    1, 1, 1, scale = 1) != 1
+     pfpareto(2^53 + xMax, min = xMax, 1, 1, 1, scale = 1) != 1
+     all.equal(pfpareto(xMin + m, min = m, 1, 1, 1, scale = 1), xMin)
+     all.equal(y <- pfpareto(1e300 + m, min = m,
+                             shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                             shape3 = 1,
+                             scale = scLrg, log = TRUE),
+               ptrbeta(1e300,
+                       shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                       shape3 = 1,
+                       scale = scLrg, log = TRUE))
+     all.equal(y,
+               c(pbeta(c(2/3, 1/2), 1, 3, lower.tail = TRUE,  log = TRUE),
+                 pbeta(2/3, 3, 1, lower.tail = FALSE, log = TRUE),
+                 pbeta(c(4/5, 1/2), 1, 3, lower.tail = TRUE,  log = TRUE),
+                 pbeta(4/5, 3, 1, lower.tail = FALSE, log = TRUE)))
+ })
> stopifnot(exprs = {
+     ppareto4(Inf,         min = 10,   1, 3, scale = xMax) == 1
+     ppareto4(2^53,        min = 0,    1, 1, scale = 1) != 1
+     ppareto4(2^53 + xMax, min = xMax, 1, 1, scale = 1) != 1
+     all.equal(ppareto4(xMin + m, min = m, 1, 1, scale = 1), xMin)
+     all.equal(y <- ppareto4(1e300 + m, min = m,
+                             shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                             scale = scLrg, log = TRUE),
+               pburr(1e300,
+                     shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                     scale = scLrg, log = TRUE))
+     all.equal(y,
+               c(log(1 - c(1/3, 1/2, 2/3)^3),
+                 log(1 - c(1/5, 1/2, 4/5)^3)))
+ })
> stopifnot(exprs = {
+     ppareto3(Inf,         min = 10,   3, scale = xMax) == 1
+     ppareto3(2^53,        min = 0,    1, scale = 1) != 1
+     ppareto3(2^53 + xMax, min = xMax, 1, scale = 1) != 1
+     all.equal(ppareto3(xMin + m, min = m, 1, scale = 1), xMin)
+     all.equal(y <- ppareto3(1e300 + m, min = m,
+                             shape = rep(c(1, 2), each = length(scLrg)),
+                             scale = scLrg, log = TRUE),
+               pllogis (1e300,
+                        shape = rep(c(1, 2), each = length(scLrg)),
+                        scale = scLrg, log = TRUE))
+     all.equal(y,
+               c(log(c(2/3, 1/2, 1/3)),
+                 log(c(4/5, 1/2, 1/5))))
+ })
> stopifnot(exprs = {
+     ppareto2(Inf,         min = 10,   3, scale = xMax) == 1
+     ppareto2(2^53,        min = 0,    1, scale = 1) != 1
+     ppareto2(2^53 + xMax, min = xMax, 1, scale = 1) != 1
+     all.equal(ppareto2(xMin + m, min = m, 1, scale = 1), xMin)
+     all.equal(y <- ppareto2(1e300 + m, min = m,
+                             shape = 3,
+                             scale = scLrg, log = TRUE),
+               ppareto (1e300,
+                        shape = 3,
+                        scale = scLrg, log = TRUE))
+     all.equal(y,
+               c(log(1 - c(1/3, 1/2, 2/3)^3)))
+ })
> 
> ## Also check that distribution functions return 0 when scale = Inf.
> stopifnot(exprs = {
+     pfpareto(x, min = m, shape1 = a, shape2 = g, shape3 = t, scale = Inf) == 0
+     ppareto4(x, min = m, shape1 = a, shape2 = g, scale = Inf) == 0
+     ppareto3(x, min = m, shape  = g, scale = Inf) == 0
+     ppareto2(x, min = m, shape  = a, scale = Inf) == 0
+ })
> 
> ## Tests for first three (positive) moments
> ##
> ## Simulation of new parameters ensuring that the first three moments
> ## exist.
> set.seed(123)                   # reset the seed
> nshpar <- 3                     # (maximum) number of shape parameters
> min <- round(rnorm(30, 2), 2)
> shpar <- replicate(30, c(3, 3, 0) + rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> for (i in seq_along(min))
+ {
+     m <- min[i]
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, t)
+     Ga <- gamma(a)
+     for (s in scpar)
+     {
+         stopifnot(exprs = {
+             All.eq(mfpareto(1, min = m,
+                             shape1 = a, shape2 = g, shape3 = t,
+                             scale = s),
+                    m * (Be + (s/m) * beta(t + 1/g, a - 1/g))/Be)
+             All.eq(mfpareto(2, min = m,
+                             shape1 = a, shape2 = g, shape3 = t,
+                             scale = s),
+                    m^2 * (Be + 2 * (s/m) * beta(t + 1/g, a - 1/g)
+                        + (s/m)^2 * beta(t + 2/g, a - 2/g))/Be)
+             All.eq(mfpareto(3, min = m,
+                             shape1 = a, shape2 = g, shape3 = t,
+                             scale = s),
+                    m^3 * (Be + 3 * (s/m) * beta(t + 1/g, a - 1/g)
+                        + 3 * (s/m)^2 * beta(t + 2/g, a - 2/g)
+                        + (s/m)^3 * beta(t + 3/g, a - 3/g))/Be)
+         })
+         stopifnot(exprs = {
+             All.eq(mpareto4(1, min = m,
+                             shape1 = a, shape2 = g,
+                             scale = s),
+                    m * (Ga + (s/m) * gamma(1 + 1/g) * gamma(a - 1/g))/Ga)
+             All.eq(mpareto4(2, min = m,
+                             shape1 = a, shape2 = g,
+                             scale = s),
+                    m^2 * (Ga + 2 * (s/m) * gamma(1 + 1/g) * gamma(a - 1/g)
+                        + (s/m)^2 * gamma(1 + 2/g) * gamma(a - 2/g))/Ga)
+             All.eq(mpareto4(3, min = m,
+                             shape1 = a, shape2 = g,
+                             scale = s),
+                    m^3 * (Ga + 3 * (s/m) * gamma(1 + 1/g) * gamma(a - 1/g)
+                        + 3 * (s/m)^2 * gamma(1 + 2/g) * gamma(a - 2/g)
+                        + (s/m)^3 * gamma(1 + 3/g) * gamma(a - 3/g))/Ga)
+         })
+         stopifnot(exprs = {
+             All.eq(mpareto3(1, min = m,
+                             shape = g,
+                             scale = s),
+                    m * (1 + (s/m) * gamma(1 + 1/g) * gamma(1 - 1/g)))
+             All.eq(mpareto3(2, min = m,
+                             shape = g,
+                             scale = s),
+                    m^2 * (1 + 2 * (s/m) * gamma(1 + 1/g) * gamma(1 - 1/g)
+                        + (s/m)^2 * gamma(1 + 2/g) * gamma(1 - 2/g)))
+             All.eq(mpareto3(3, min = m,
+                             shape = g,
+                             scale = s),
+                    m^3 * (1 + 3 * (s/m) * gamma(1 + 1/g) * gamma(1 - 1/g)
+                        + 3 * (s/m)^2 * gamma(1 + 2/g) * gamma(1 - 2/g)
+                        + (s/m)^3 * gamma(1 + 3/g) * gamma(1 - 3/g)))
+         })
+         stopifnot(exprs = {
+             All.eq(mpareto2(1, min = m,
+                             shape = a,
+                             scale = s),
+                    m * (Ga + (s/m) * gamma(1 + 1) * gamma(a - 1))/Ga)
+             All.eq(mpareto2(2, min = m,
+                             shape = a,
+                             scale = s),
+                    m^2 * (Ga + 2 * (s/m) * gamma(1 + 1) * gamma(a - 1)
+                        + (s/m)^2 * gamma(1 + 2) * gamma(a - 2))/Ga)
+             All.eq(mpareto2(3, min = m,
+                             shape = a,
+                             scale = s),
+                    m^3 * (Ga + 3 * (s/m) * gamma(1 + 1) * gamma(a - 1)
+                        + 3 * (s/m)^2 * gamma(1 + 2) * gamma(a - 2)
+                        + (s/m)^3 * gamma(1 + 3) * gamma(a - 3))/Ga)
+         })
+     }
+ }
> 
> ## Tests for first three limited moments
> ##
> ## Limits are taken from quantiles of each distribution.
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)
> for (i in seq_along(min))
+ {
+     m <- min[i]
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Ga <- gamma(a)
+     Gt <- gamma(t)
+     for (s in scpar)
+     {
+         limit <- qfpareto(q, min = m,
+                           shape1 = a, shape2 = g, shape3 = t,
+                           scale = s)
+         y <- (limit - m)/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             All.eq(levfpareto(limit, order = 1, min = m,
+                               shape1 = a, shape2 = g, shape3 = t,
+                               scale = s),
+                    m * (betaint(u, t, a) + (s/m) * betaint(u, t + 1/g, a - 1/g))/(Ga * Gt) +
+                    limit * pbeta(u, t, a, lower = FALSE))
+             All.eq(levfpareto(limit, order = 2, min = m,
+                               shape1 = a, shape2 = g, shape3 = t,
+                               scale = s),
+                    m^2 * (betaint(u, t, a) + 2 * (s/m) * betaint(u, t + 1/g, a - 1/g)
+                        + (s/m)^2 * betaint(u, t + 2/g, a - 2/g))/(Ga * Gt) +
+                    limit^2 * pbeta(u, t, a, lower = FALSE))
+             All.eq(levfpareto(limit, order = 3, min = m,
+                               shape1 = a, shape2 = g, shape3 = t,
+                               scale = s),
+                    m^3 * (betaint(u, t, a) + 3 * (s/m) * betaint(u, t + 1/g, a - 1/g)
+                        + 3 * (s/m)^2 * betaint(u, t + 2/g, a - 2/g)
+                        + (s/m)^3 * betaint(u, t + 3/g, a - 3/g))/(Ga * Gt) +
+                    limit^3 * pbeta(u, t, a, lower = FALSE))
+         })
+         limit <- qpareto4(q, min = m,
+                           shape1 = a, shape2 = g,
+                           scale = s)
+         y <- (limit - m)/s
+         u <- 1/(1 + y^g)
+         u1m <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             All.eq(levpareto4(limit, order = 1, min = m,
+                               shape1 = a, shape2 = g,
+                               scale = s),
+                    m * (betaint(u1m, 1, a) + (s/m) * betaint(u1m, 1 + 1/g, a - 1/g))/Ga +
+                    limit * u^a)
+             All.eq(levpareto4(limit, order = 2, min = m,
+                               shape1 = a, shape2 = g,
+                               scale = s),
+                    m^2 * (betaint(u1m, 1, a) + 2 * (s/m) * betaint(u1m, 1 + 1/g, a - 1/g)
+                        + (s/m)^2 * betaint(u1m, 1 + 2/g, a - 2/g))/Ga +
+                    limit^2 * u^a)
+             All.eq(levpareto4(limit, order = 3, min = m,
+                               shape1 = a, shape2 = g,
+                               scale = s),
+                    m^3 * (betaint(u1m, 1, a) + 3 * (s/m) * betaint(u1m, 1 + 1/g, a - 1/g)
+                        + 3 * (s/m)^2 * betaint(u1m, 1 + 2/g, a - 2/g)
+                        + (s/m)^3 * betaint(u1m, 1 + 3/g, a - 3/g))/Ga +
+                    limit^3 * u^a)
+         })
+         limit <- qpareto3(q, min = m,
+                           shape = g,
+                           scale = s)
+         y <- (limit - m)/s
+         u <- 1/(1 + y^(-g))
+         u1m <- 1/(1 + y^g)
+         stopifnot(exprs = {
+             All.eq(levpareto3(limit, order = 1, min = m,
+                               shape = g,
+                               scale = s),
+                    m * (u + (s/m) * betaint(u, 1 + 1/g, 1 - 1/g)) +
+                    limit * u1m)
+             All.eq(levpareto3(limit, order = 2, min = m,
+                               shape = g,
+                               scale = s),
+                    m^2 * (u + 2 * (s/m) * betaint(u, 1 + 1/g, 1 - 1/g)
+                        + (s/m)^2 * betaint(u, 1 + 2/g, 1 - 2/g)) +
+                    limit^2 * u1m)
+             All.eq(levpareto3(limit, order = 3, min = m,
+                               shape = g,
+                               scale = s),
+                    m^3 * (u + 3 * (s/m) * betaint(u, 1 + 1/g, 1 - 1/g)
+                        + 3 * (s/m)^2 * betaint(u, 1 + 2/g, 1 - 2/g)
+                        + (s/m)^3 * betaint(u, 1 + 3/g, 1 - 3/g)) +
+                    limit^3 * u1m)
+         })
+         limit <- qpareto2(q, min = m,
+                           shape = a,
+                           scale = s)
+         y <- (limit - m)/s
+         u <- 1/(1 + y)
+         u1m <- 1/(1 + y^(-1))
+         stopifnot(exprs = {
+             All.eq(levpareto2(limit, order = 1, min = m,
+                               shape = a,
+                               scale = s),
+                    m * (betaint(u1m, 1, a) + (s/m) * betaint(u1m, 1 + 1, a - 1))/Ga +
+                    limit * u^a)
+             All.eq(levpareto2(limit, order = 2, min = m,
+                               shape = a,
+                               scale = s),
+                    m^2 * (betaint(u1m, 1, a) + 2 * (s/m) * betaint(u1m, 1 + 1, a - 1)
+                        + (s/m)^2 * betaint(u1m, 1 + 2, a - 2))/Ga +
+                    limit^2 * u^a)
+             All.eq(levpareto2(limit, order = 3, min = m,
+                               shape = a,
+                               scale = s),
+                    m^3 * (betaint(u1m, 1, a) + 3 * (s/m) * betaint(u1m, 1 + 1, a - 1)
+                        + 3 * (s/m)^2 * betaint(u1m, 1 + 2, a - 2)
+                        + (s/m)^3 * betaint(u1m, 1 + 3, a - 3))/Ga +
+                    limit^3 * u^a)
+         })
+     }
+ }
> 
> ##
> ## TRANSFORMED BETA FAMILY
> ##
> 
> ## Density: first check that functions return 0 when scale = Inf, and
> ## when x = scale = Inf.
> stopifnot(exprs = {
+     dtrbeta      (c(42, Inf), shape1 = 2, shape2 = 3, shape3 = 4, scale = Inf) == c(0, 0)
+     dburr        (c(42, Inf), shape1 = 2, shape2 = 3, scale = Inf) == c(0, 0)
+     dllogis      (c(42, Inf), shape  = 3, scale = Inf) == c(0, 0)
+     dparalogis   (c(42, Inf), shape  = 2, scale = Inf) == c(0, 0)
+     dgenpareto   (c(42, Inf), shape1 = 2, shape2 = 4, scale = Inf) == c(0, 0)
+     dpareto      (c(42, Inf), shape  = 2, scale = Inf) == c(0, 0)
+     dinvburr     (c(42, Inf), shape1 = 4, shape2 = 3, scale = Inf) == c(0, 0)
+     dinvpareto   (c(42, Inf), shape  = 4, scale = Inf) == c(0, 0)
+     dinvparalogis(c(42, Inf), shape  = 4, scale = Inf) == c(0, 0)
+ })
> 
> ## Next test density functions for an array of standard values.
> set.seed(123)                   # reset the seed
> nshpar <- 3                     # (maximum) number of shape parameters
> shpar <- replicate(30, rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, t)
+     for (s in scpar)
+     {
+         x <- rtrbeta(100, shape1 = a, shape2 = g, shape3 = t, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             all.equal(d1 <- dtrbeta(x, shape1 = a, shape2 = g, shape3 = t,
+                                     scale = s),
+                       d2 <- dtrbeta(y, shape1 = a, shape2 = g, shape3 = t,
+                                     scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * y^(g*t - 1)/(s * Be * (1 + y^g)^(a + t)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * u^t * (1 - u)^a/(x * Be),
+                       tolerance = 1e-10)
+         })
+         x <- rburr(100, shape1 = a, shape2 = g, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^g)
+         stopifnot(exprs = {
+             all.equal(d1 <- dburr(x, shape1 = a, shape2 = g,
+                                   scale = s),
+                       d2 <- dburr(y, shape1 = a, shape2 = g,
+                                   scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * g * y^(g - 1)/(s * (1 + y^g)^(a + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * g * u^a * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rllogis(100, shape = g, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             all.equal(d1 <- dllogis(x, shape = g,
+                                     scale = s),
+                       d2 <- dllogis(y, shape = g,
+                                     scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * y^(g - 1)/(s * (1 + y^g)^2),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       g * u * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rparalogis(100, shape = a, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^a)
+         stopifnot(exprs = {
+             all.equal(d1 <- dparalogis(x, shape = a,
+                                        scale = s),
+                       d2 <- dparalogis(y, shape = a,
+                                        scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a^2 * y^(a - 1)/(s * (1 + y^a)^(a + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a^2 * u^a * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rgenpareto(100, shape1 = a, shape2 = t, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-1))
+         stopifnot(exprs = {
+             all.equal(d1 <- dgenpareto(x, shape1 = a, shape2 = t,
+                                        scale = s),
+                       d2 <- dgenpareto(y, shape1 = a, shape2 = t,
+                                        scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       y^(t - 1)/(s * Be * (1 + y)^(a + t)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       u^t * (1 - u)^a/(x * Be),
+                       tolerance = 1e-10)
+         })
+         x <- rpareto(100, shape = a, scale = s)
+         y <- x/s
+         u <- 1/(1 + y)
+         stopifnot(exprs = {
+             all.equal(d1 <- dpareto(x, shape = a,
+                                     scale = s),
+                       d2 <- dpareto(y, shape = a,
+                                     scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a/(s * (1 + y)^(a + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       a * u^a * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rpareto1(100, min = s, shape = a)
+         stopifnot(exprs = {
+             all.equal(d1 <- dpareto1(x, min = s, shape = a),
+                       a * s^a/(x^(a + 1)),
+                       tolerance = 1e-10)
+         })
+         x <- rinvburr(100, shape1 = t, shape2 = g, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-g))
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvburr(x, shape1 = t, shape2 = g,
+                                      scale = s),
+                       d2 <- dinvburr(y, shape1 = t, shape2 = g,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * g * y^(g*t - 1)/(s * (1 + y^g)^(t + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * g * u^t * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rinvpareto(100, shape = t, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-1))
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvpareto(x, shape = t,
+                                  scale = s),
+                       d2 <- dinvpareto(y, shape = t,
+                                        scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * y^(t - 1)/(s * (1 + y)^(t + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * u^t * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+         x <- rinvparalogis(100, shape = t, scale = s)
+         y <- x/s
+         u <- 1/(1 + y^(-t))
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvparalogis(x, shape = t,
+                                           scale = s),
+                       d2 <- dinvparalogis(y, shape = t,
+                                           scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t^2 * y^(t^2 - 1)/(s * (1 + y^t)^(t + 1)),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t^2 * u^t * (1 - u)/x,
+                       tolerance = 1e-10)
+         })
+     }
+ }
> 
> ## Tests on the cumulative distribution function.
> ##
> ## Note: when shape1 = shape3 = 1, the underlying beta distribution is
> ## a uniform. Therefore, ptrbeta(x, 1, shape2, 1, scale) should return
> ## the value of u = v/(1 + v), v = (x/scale)^shape2.
> ##
> ## x = 2/Meps = 2^53 (with, shape2 = scale = 1) is the value where the
> ## cdf would jump to 1 if we weren't using the trick to compute the
> ## cdf with pbeta(1 - u, ..., lower = FALSE).
> scLrg <- 1e300 * c(0.5, 1, 2)
> stopifnot(exprs = {
+     ptrbeta(Inf,  1, 2, 3, scale = xMax) == 1
+     ptrbeta(2^53, 1, 1, 1, scale = 1) != 1
+     all.equal(ptrbeta(xMin, 1, 1, 1, scale = 1), xMin)
+     all.equal(ptrbeta(1e300,
+                       shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                       shape3 = 1,
+                       scale = scLrg, log = TRUE),
+               c(pbeta(c(2/3, 1/2), 1, 3, lower.tail = TRUE,  log = TRUE),
+                 pbeta(2/3, 3, 1, lower.tail = FALSE, log = TRUE),
+                 pbeta(c(4/5, 1/2), 1, 3, lower.tail = TRUE,  log = TRUE),
+                 pbeta(4/5, 3, 1, lower.tail = FALSE, log = TRUE)))
+ })
> stopifnot(exprs = {
+     pburr(Inf,  1, 3, scale = xMax) == 1
+     pburr(2^53, 1, 1, scale = 1) != 1
+     all.equal(pburr(xMin, 1, 1, scale = 1), xMin)
+     all.equal(pburr(1e300,
+                     shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                     scale = scLrg, log = TRUE),
+               c(log(1 - c(1/3, 1/2, 2/3)^3),
+                 log(1 - c(1/5, 1/2, 4/5)^3)))
+ })
> stopifnot(exprs = {
+     pllogis(Inf,  3, scale = xMax) == 1
+     pllogis(2^53, 1, scale = 1) != 1
+     all.equal(pllogis(xMin, 1, scale = 1), xMin)
+     all.equal(pllogis(1e300,
+                       shape = rep(c(1, 2), each = length(scLrg)),
+                       scale = scLrg, log = TRUE),
+               c(log(c(2/3, 1/2, 1/3)),
+                 log(c(4/5, 1/2, 1/5))))
+ })
> stopifnot(exprs = {
+     pparalogis(Inf,  3, scale = xMax) == 1
+     pparalogis(2^53, 1, scale = 1) != 1
+     all.equal(pparalogis(xMin, 1, scale = 1), xMin)
+     all.equal(pparalogis(1e300,
+                          shape = rep(c(2, 3), each = length(scLrg)),
+                          scale = scLrg, log = TRUE),
+               c(log(1 - c(1/5, 1/2, 4/5)^2),
+                 log(1 - c(1/9, 1/2, 8/9)^3)))
+ })
> stopifnot(exprs = {
+     pgenpareto(Inf,  1, 3, scale = xMax) == 1
+     pgenpareto(2^53, 1, 1, scale = 1) != 1
+     all.equal(pgenpareto(xMin, 1, 1, scale = 1), xMin)
+     all.equal(pgenpareto(1e300,
+                          shape1 = 3, shape2 = 1,
+                          scale = scLrg, log = TRUE),
+               c(pbeta(c(2/3, 1/2), 1, 3, lower.tail = TRUE,  log = TRUE),
+                 pbeta(2/3, 3, 1, lower.tail = FALSE, log = TRUE)))
+ })
> stopifnot(exprs = {
+     ppareto(Inf,  3, scale = xMax) == 1
+     ppareto(2^53, 1, scale = 1) != 1
+     all.equal(ppareto(xMin, 1, scale = 1), xMin)
+     all.equal(ppareto(1e300,
+                       shape = 3,
+                       scale = scLrg, log = TRUE),
+               c(log(1 - c(1/3, 1/2, 2/3)^3)))
+ })
> stopifnot(exprs = {
+     ppareto1(Inf,  3, min = xMax) == 1
+     ppareto1(2^53, 1, min = 1) != 1
+     all.equal(ppareto1(xMin, 1, min = 1), xMin)
+     all.equal(ppareto1(1e300,
+                       shape = 3,
+                       min = 1e300 * c(0.001, 0.1, 0.5), log = TRUE),
+               c(log(1 - c(0.001, 0.1, 0.5)^3)))
+ })
> stopifnot(exprs = {
+     pinvburr(Inf,  1, 3, scale = xMax) == 1
+     pinvburr(2^53, 1, 1, scale = 1) != 1
+     all.equal(pinvburr(xMin, 1, 1, scale = 1), xMin)
+     all.equal(pinvburr(1e300,
+                        shape1 = 3, shape2 = rep(c(1, 2), each = length(scLrg)),
+                        scale = scLrg, log = TRUE),
+               c(log(c(2/3, 1/2, 1/3)^3),
+                 log(c(4/5, 1/2, 1/5)^3)))
+ })
> stopifnot(exprs = {
+     pinvpareto(Inf,  3, scale = xMax) == 1
+     pinvpareto(2^53, 1, scale = 1) != 1
+     all.equal(pinvpareto(xMin, 1, scale = 1), xMin)
+     all.equal(pinvpareto(1e300,
+                          shape = 3,
+                          scale = scLrg, log = TRUE),
+               c(log(c(2/3, 1/2, 1/3)^3)))
+ })
> stopifnot(exprs = {
+     pinvparalogis(Inf,  3, scale = xMax) == 1
+     pinvparalogis(2^53, 1, scale = 1) != 1
+     all.equal(pinvparalogis(xMin, 1, scale = 1), xMin)
+     all.equal(pinvparalogis(1e300,
+                             shape = rep(c(2, 3), each = length(scLrg)),
+                             scale = scLrg, log = TRUE),
+               c(log(c(4/5, 1/2, 1/5)^2),
+                 log(c(8/9, 1/2, 1/9)^3)))
+ })
> 
> ## Also check that distribution functions return 0 when scale = Inf.
> stopifnot(exprs = {
+     ptrbeta      (x, shape1 = a, shape2 = g, shape3 = t, scale = Inf) == 0
+     pburr        (x, shape1 = a, shape2 = g, scale = Inf) == 0
+     pllogis      (x, shape  = g, scale = Inf) == 0
+     pparalogis   (x, shape  = a, scale = Inf) == 0
+     pgenpareto   (x, shape1 = a, shape2 = t, scale = Inf) == 0
+     ppareto      (x, shape  = a, scale = Inf) == 0
+     pinvburr     (x, shape1 = t, shape2 = g, scale = Inf) == 0
+     pinvpareto   (x, shape  = t, scale = Inf) == 0
+     pinvparalogis(x, shape  = t, scale = Inf) == 0
+ })
> 
> ## Tests for first three positive moments and first two negative
> ## moments.
> ##
> ## Simulation of new parameters ensuring that said moments exist.
> set.seed(123)                   # reset the seed
> nshpar <- 3                     # (maximum) number of shape parameters
> shpar <- replicate(30, c(3, 3, 3) + rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> k <- c(-2, -1, 1, 2, 3)         # orders
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, t)
+     Ga <- gamma(a)
+     for (s in scpar)
+     {
+         stopifnot(exprs = {
+             All.eq(mtrbeta(k, shape1 = a, shape2 = g, shape3 = t, scale = s),
+                    s^k * beta(t + k/g, a - k/g)/Be)
+             All.eq(mburr(k, shape1 = a, shape2 = g, scale = s),
+                    s^k * gamma(1 + k/g) * gamma(a - k/g)/Ga)
+             All.eq(mllogis(k, shape = g, scale = s),
+                    s^k * gamma(1 + k/g) * gamma(1 - k/g))
+             All.eq(mparalogis(k, shape = a, scale = s),
+                    s^k * gamma(1 + k/a) * gamma(a - k/a)/Ga)
+             All.eq(mgenpareto(k, shape1 = a, shape2 = t, scale = s),
+                    s^k * beta(t + k, a - k)/Be)
+             All.eq(mpareto(k[k > -1], shape = a, scale = s),
+                    s^k[k > -1] * gamma(1 + k[k > -1]) * gamma(a - k[k > -1])/Ga)
+             All.eq(mpareto1(k, shape = a, min = s),
+                    s^k * a/(a - k))
+             All.eq(minvburr(k, shape1 = a, shape2 = g, scale = s),
+                    s^k * gamma(a + k/g) * gamma(1 - k/g)/Ga)
+             All.eq(minvpareto(k[k < 1], shape = a, scale = s),
+                 s^k[k < 1] * gamma(a + k[k < 1]) * gamma(1 - k[k < 1])/Ga)
+             All.eq(minvparalogis(k, shape = a, scale = s),
+                    s^k * gamma(a + k/a) * gamma(1 - k/a)/Ga)
+         })
+     }
+ }
> 
> ## Tests for first three positive limited moments and first two
> ## negative limited moments.
> ##
> ## Limits are taken from quantiles of each distribution.
> order <- c(-2, -1, 1, 2, 3)             # orders
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)     # quantiles
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Ga <- gamma(a)
+     Gt <- gamma(t)
+     for (s in scpar)
+     {
+         limit <- qtrbeta(q, shape1 = a, shape2 = g, shape3 = t, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-g))
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levtrbeta(limit, order = k, shape1 = a, shape2 = g, shape3 = t, scale = s),
+                        s^k * betaint(u, t + k/g, a - k/g)/(Ga * Gt) +
+                        limit^k * pbeta(u, t, a, lower = FALSE))
+             })
+         limit <- qburr(q, shape1 = a, shape2 = g, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^g)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levburr(limit, order = k, shape1 = a, shape2 = g, scale = s),
+                        s^k * betaint(1 - u, 1 + k/g, a - k/g)/Ga +
+                        limit^k * u^a)
+             })
+         limit <- qllogis(q, shape = g, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-g))
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levllogis(limit, order = k, shape = g, scale = s),
+                        s^k * betaint(u, 1 + k/g, 1 - k/g) +
+                        limit^k * (1 - u))
+             })
+         limit <- qparalogis(q, shape = a, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^a)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levparalogis(limit, order = k, shape = a, scale = s),
+                        s^k * betaint(1 - u, 1 + k/a, a - k/a)/Ga +
+                        limit^k * u^a)
+             })
+         limit <- qgenpareto(q, shape1 = a, shape2 = t, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-1))
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levgenpareto(limit, order = k, shape1 = a, shape2 = t, scale = s),
+                        s^k * betaint(u, t + k, a - k)/(Ga * Gt) +
+                        limit^k * pbeta(u, t, a, lower = FALSE))
+             })
+         limit <- qpareto(q, shape = a, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y)
+         for (k in order[order > -1])
+             stopifnot(exprs = {
+                 All.eq(levpareto(limit, order = k, shape = a, scale = s),
+                        s^k * betaint(1 - u, 1 + k, a - k)/Ga +
+                        limit^k * u^a)
+         })
+         limit <- qpareto1(q, shape = a, min = s)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levpareto1(limit, order = k, shape = a, min = s),
+                        s^k * a/(a - k) - k * s^a/((a - k) * limit^(a - k)))
+             })
+         limit <- qinvburr(q, shape1 = a, shape2 = g, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-g))
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvburr(limit, order = k, shape1 = a, shape2 = g, scale = s),
+                        s^k * betaint(u, a + k/g, 1 - k/g)/Ga +
+                        limit^k * (1 - u^a))
+             })
+         limit <- qinvpareto(q, shape = a, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-1))
+         for (k in order[order < 1])
+             stopifnot(exprs = {
+                 All.eq(levinvpareto(limit, order = k, shape = a, scale = s),
+                        s^k * a *
+                        sapply(u,
+                               function(upper)
+                                   integrate(function(x) x^(a+k-1) * (1-x)^(-k),
+                                             lower = 0, upper = upper)$value) +
+                        limit^k * (1 - u^a))
+             })
+         limit <- qinvparalogis(q, shape = a, scale = s)
+         y <- limit/s
+         u <- 1/(1 + y^(-a))
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvparalogis(limit, order = k, shape = a, scale = s),
+                        s^k * betaint(u, a + k/a, 1 - k/a)/Ga +
+                        limit^k * (1 - u^a))
+         })
+     }
+ }
> 
> ##
> ## TRANSFORMED GAMMA AND INVERSE TRANSFORMED GAMMA FAMILIES
> ##
> 
> ## Density: first check that functions return 0 when scale = Inf, and
> ## when x = scale = Inf (transformed gamma), or when scale = 0 and
> ## when x = scale = 0 (inverse distributions).
> stopifnot(exprs = {
+     dtrgamma   (c(42, Inf), shape1 = 2, shape2 = 3, scale = Inf) == c(0, 0)
+     dinvtrgamma(c(42, 0), shape1 = 2, shape2 = 3, scale = 0) == c(0, 0)
+     dinvgamma  (c(42, 0), shape  = 2, scale = 0) == c(0, 0)
+     dinvweibull(c(42, 0), shape = 3, scale = 0) == c(0, 0)
+     dinvexp    (c(42, 0), scale = 0) == c(0, 0)
+ })
> 
> ## Tests on the density
> set.seed(123)                   # reset the seed
> nshpar <- 2                     # (maximum) number of shape parameters
> shpar <- replicate(30, rgamma(nshpar, 5), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; t <- shpar[[c(i, 2)]]
+     Ga <- gamma(a)
+     for (s in scpar)
+     {
+         x <- rtrgamma(100, shape1 = a, shape2 = t, scale = s)
+         y <- x/s
+         u <- y^t
+         stopifnot(exprs = {
+             all.equal(d1 <- dtrgamma(x, shape1 = a, shape2 = t,
+                                     scale = s),
+                       d2 <- dtrgamma(y, shape1 = a, shape2 = t,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       t/(Ga * s^(a * t)) * x^(a * t - 1) * exp(-u),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t/(Ga * x) * u^a * exp(-u),
+                       tolerance = 1e-10)
+         })
+         x <- rinvtrgamma(100, shape1 = a, shape2 = t, scale = s)
+         y <- x/s
+         u <- y^(-t)
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvtrgamma(x, shape1 = a, shape2 = t,
+                                         scale = s),
+                       d2 <- dinvtrgamma(y, shape1 = a, shape2 = t,
+                                         scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       t * s^(a * t)/(Ga * x^(a * t + 1)) * exp(-u),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t/(Ga * x) * u^a * exp(-u),
+                       tolerance = 1e-10)
+         })
+         x <- rinvgamma(100, shape = a, scale = s)
+         y <- x/s
+         u <- y^(-1)
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvgamma(x, shape = a, scale = s),
+                       d2 <- dinvgamma(y, shape = a, scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       s^a/(Ga * x^(a + 1)) * exp(-u),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       1/(Ga * x) * u^a * exp(-u),
+                       tolerance = 1e-10)
+         })
+         x <- rinvweibull(100, shape = t, scale = s)
+         y <- x/s
+         u <- y^(-t)
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvweibull(x, shape = t, scale = s),
+                       d2 <- dinvweibull(y, shape = t, scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       t * s^t/x^(t + 1) * exp(-u),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t/x * u * exp(-u),
+                       tolerance = 1e-10)
+         })
+         x <- rinvexp(100, scale = s)
+         y <- x/s
+         u <- y^(-1)
+         stopifnot(exprs = {
+             all.equal(d1 <- dinvexp(x, scale = s),
+                       d2 <- dinvexp(y, scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d2,
+                       s/x^2 * exp(-u),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       1/x * u * exp(-u),
+                       tolerance = 1e-10)
+         })
+     }
+ }
> 
> ## Tests on the cumulative distribution function.
> scLrg <- c(2, 100, 1e300 * c(0.1, 1, 10, 100), 1e307, xMax, Inf)
> stopifnot(exprs = {
+     ptrgamma(Inf,  2, 3, scale = xMax) == 1
+     ptrgamma(xMax, 2, 3, scale = xMax) == pgamma(1, 2, 1)
+     ptrgamma(xMin, 2, 1, scale = 1)    == pgamma(xMin, 2, 1)
+     all.equal(ptrgamma(1e300, shape1 = 2, shape2 = 1, scale = scLrg, log = TRUE),
+               pgamma(c(5e299, 1e+298, 10, 1, 0.1, 0.01, 1e-7, 1e+300/xMax, 0),
+                      2, 1, log = TRUE))
+ })
> scLrg <- c(2, 100, 1e300 * c(0.1, 1, 10, 100), 1e307, xMax, 0)
> stopifnot(exprs = {
+     pinvtrgamma(Inf,  2, 3, scale = xMax) == 1
+     pinvtrgamma(xMax, 2, 3, scale = xMax) == pgamma(1, 2, 1, lower = FALSE)
+     pinvtrgamma(xMin, 2, 1, scale = 1)    == pgamma(1/xMin, 2, 1, lower = FALSE)
+     all.equal(pinvtrgamma(1e300, shape1 = 2, shape2 = 1, scale = scLrg, log = TRUE),
+               pgamma(c(2e-300, 1e-298, 0.1, 1, 10, 100, 1e+7, xMax/1e+300, 0),
+                      2, 1, lower = FALSE, log = TRUE))
+ })
> stopifnot(exprs = {
+     pinvgamma(Inf,  2, scale = xMax) == 1
+     pinvgamma(xMax, 2, scale = xMax) == pgamma(1, 2, 1, lower = FALSE)
+     pinvgamma(xMin, 2, scale = 1)    == pgamma(1/xMin, 2, 1, lower = FALSE)
+     all.equal(pinvgamma(1e300, shape = 2, scale = scLrg, log = TRUE),
+               pgamma(c(2e-300, 1e-298, 0.1, 1, 10, 100, 1e+7, xMax/1e+300, 0),
+                      2, 1, lower = FALSE, log = TRUE))
+ })
> stopifnot(exprs = {
+     pinvweibull(Inf,  3, scale = xMax) == 1
+     pinvweibull(xMax, 3, scale = xMax) == exp(-1)
+     pinvweibull(xMin, 1, scale = 1)    == exp(-1/xMin)
+     all.equal(pinvweibull(1e300, shape = 1, scale = scLrg, log = TRUE),
+               -c(2e-300, 1e-298, 0.1, 1, 10, 100, 1e+7, xMax/1e+300, 0))
+ })
> stopifnot(exprs = {
+     pinvexp(Inf,  3, scale = xMax) == 1
+     pinvexp(xMax, 3, scale = xMax) == exp(-1)
+     pinvexp(xMin, 1, scale = 1)    == exp(-1/xMin)
+     all.equal(pinvexp(1e300, scale = scLrg, log = TRUE),
+               -c(2e-300, 1e-298, 0.1, 1, 10, 100, 1e+7, xMax/1e+300, 0))
+ })
> 
> ## Tests for first three positive moments and first two negative
> ## moments. (Including for the Gamma, Weibull and Exponential
> ## distributions of base R.)
> ##
> ## Simulation of new parameters ensuring that said moments exist.
> set.seed(123)                   # reset the seed
> nshpar <- 2                     # (maximum) number of shape parameters
> shpar <- replicate(30, c(3, 3) + rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> k <- c(-2, -1, 1, 2, 3)         # orders
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; t <- shpar[[c(i, 2)]]
+     Ga <- gamma(a)
+     for (s in scpar)
+     {
+         stopifnot(exprs = {
+             All.eq(mtrgamma(k, shape1 = a, shape2 = t, scale = s),
+                    s^k * gamma(a + k/t)/Ga)
+             All.eq(mgamma(k, shape = a, scale = s),
+                    s^k * gamma(a + k)/Ga)
+             All.eq(mweibull(k, shape = t, scale = s),
+                    s^k * gamma(1 + k/t))
+             All.eq(mexp(k[k > -1], rate = 1/s),
+                    s^k[k > -1] * gamma(1 + k[k > -1]))
+             All.eq(minvtrgamma(k, shape1 = a, shape2 = t, scale = s),
+                    s^k * gamma(a - k/t)/Ga)
+             All.eq(minvgamma(k, shape = a, scale = s),
+                    s^k * gamma(a - k)/Ga)
+             All.eq(minvweibull(k, shape = t, scale = s),
+                    s^k * gamma(1 - k/t))
+             All.eq(minvexp(k[k < 1], scale = s),
+                    s^k[k < 1] * gamma(1 - k[k < 1]))
+         })
+     }
+ }
> 
> ## Tests for first three positive limited moments and first two
> ## negative limited moments. (Including for the Gamma, Weibull and
> ## Exponential distributions of base R.)
> ##
> ## Limits are taken from quantiles of each distribution.
> order <- c(-2, -1, 1, 2, 3)             # orders
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)     # quantiles
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; t <- shpar[[c(i, 2)]]
+     Ga <- gamma(a)
+     for (s in scpar)
+     {
+         limit <- qtrgamma(q, shape1 = a, shape2 = t, scale = s)
+         y <- limit/s
+         u <- y^t
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levtrgamma(limit, order = k, shape1 = a, shape2 = t, scale = s),
+                        s^k * gamma(a + k/t)/Ga * pgamma(u, a + k/t, scale = 1) +
+                        limit^k * pgamma(u, a, scale = 1, lower = FALSE))
+             })
+         limit <- qgamma(q, shape = a, scale = s)
+         y <- limit/s
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levgamma(limit, order = k, shape = a, scale = s),
+                        s^k * gamma(a + k)/Ga * pgamma(y, a + k, scale = 1) +
+                        limit^k * pgamma(y, a, scale = 1, lower = FALSE))
+             })
+         limit <- qweibull(q, shape = t, scale = s)
+         y <- limit/s
+         u <- y^t
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levweibull(limit, order = k, shape = t, scale = s),
+                        s^k * gamma(1 + k/t) * pgamma(u, 1 + k/t, scale = 1) +
+                        limit^k * pgamma(u, 1, scale = 1, lower = FALSE))
+             })
+         limit <- qexp(q, rate = 1/s)
+         y <- limit/s
+         for (k in order[order > -1])
+             stopifnot(exprs = {
+                 All.eq(levexp(limit, order = k, rate = 1/s),
+                        s^k * gamma(1 + k) * pgamma(y, 1 + k, scale = 1) +
+                        limit^k * pgamma(y, 1, scale = 1, lower = FALSE))
+             })
+         limit <- qinvtrgamma(q, shape1 = a, shape2 = t, scale = s)
+         y <- limit/s
+         u <- y^(-t)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvtrgamma(limit, order = k, shape1 = a, shape2 = t, scale = s),
+                        s^k * (gammainc(a - k/t, u)/Ga) +
+                        limit^k * pgamma(u, a, scale = 1))
+             })
+         limit <- qinvgamma(q, shape = a, scale = s)
+         y <- limit/s
+         u <- y^(-1)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvgamma(limit, order = k, shape = a, scale = s),
+                        s^k * (gammainc(a - k, u)/Ga) +
+                        limit^k * pgamma(u, a, scale = 1))
+             })
+         limit <- qinvweibull(q, shape = t, scale = s)
+         y <- limit/s
+         u <- y^(-t)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvweibull(limit, order = k, shape = t, scale = s),
+                        s^k * gammainc(1 - k/t, u) +
+                        limit^k * (-expm1(-u)))
+             })
+         limit <- qinvexp(q, scale = s)
+         y <- limit/s
+         u <- y^(-1)
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levinvexp(limit, order = k, scale = s),
+                        s^k * gammainc(1 - k, u) +
+                        limit^k * (-expm1(-u)))
+             })
+     }
+ }
> 
> ##
> ## OTHER DISTRIBUTIONS
> ##
> 
> ## Distributions in this category are quite different, so let's treat
> ## them separately.
> 
> ## LOGGAMMA
> 
> ## Tests on the density.
> stopifnot(exprs = {
+     dlgamma(c(42, Inf), shapelog = 2, ratelog = 0) == c(0, 0)
+ })
> assertWarning(stopifnot(exprs = {
+     is.nan(dlgamma(c(0, 42, Inf), shapelog = 2, ratelog = Inf))
+ }))
> x <- rlgamma(100, shapelog = 2, ratelog = 1)
> for(a in round(rlnorm(30), 2))
+ {
+     Ga <- gamma(a)
+     for(r in round(rlnorm(30), 2))
+ 	stopifnot(exprs = {
+             All.eq(dlgamma(x, shapelog = a, ratelog = r),
+                    r^a * (log(x))^(a - 1)/(Ga * x^(r + 1)))
+         })
+ }
> 
> ## Tests on the cumulative distribution function.
> assertWarning(stopifnot(exprs = {
+     is.nan(plgamma(Inf,   1, ratelog = Inf))
+     is.nan(plgamma(Inf, Inf, ratelog = Inf))
+ }))
> scLrg <- log(c(2, 100, 1e300 * c(0.1, 1, 10, 100), 1e307, xMax, Inf))
> stopifnot(exprs = {
+     plgamma(Inf,  2, ratelog = xMax) == 1
+     plgamma(xMax, 2, ratelog = 0) == 0
+     all.equal(plgamma(1e300, 2, ratelog = 1/scLrg, log = TRUE),
+               pgamma(log(1e300), 2, scale = scLrg, log = TRUE))
+ })
> 
> ## Tests for first three positive moments and first two negative
> ## moments.
> k <- c(-2, -1, 1, 2, 3)         # orders
> for(a in round(rlnorm(30), 2))
+ {
+     Ga <- gamma(a)
+     for(r in 3 + round(rlnorm(30), 2))
+ 	stopifnot(exprs = {
+             All.eq(mlgamma(k, shapelog = a, ratelog = r),
+                    (1 - k/r)^(-a))
+         })
+ }
> 
> ## Tests for first three positive limited moments and first two
> ## negative limited moments.
> order <- c(-2, -1, 1, 2, 3)             # orders
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)     # quantiles
> for(a in round(rlnorm(30), 2))
+ {
+     Ga <- gamma(a)
+     for(r in 3 + round(rlnorm(30), 2))
+     {
+         limit <- qlgamma(q, shapelog = a, ratelog = r)
+         for (k in order)
+         {
+             u <- log(limit)
+         stopifnot(exprs = {
+             All.eq(levlgamma(limit, order = k, shapelog = a, ratelog = r),
+                    (1 - k/r)^(-a) * pgamma((r - k) * u, a, scale = 1) +
+                    limit^k * pgamma(r * u, a, scale = 1,lower = FALSE))
+             })
+         }
+     }
+ }
> 
> ## GUMBEL
> 
> ## Tests on the density.
> stopifnot(exprs = {
+     dgumbel(c(1, 3, Inf), alpha = 2,   scale = Inf) == c(0, 0, 0)
+     dgumbel(c(1, 2, 3),   alpha = 2,   scale = 0) == c(0, Inf, 0)
+     dgumbel(c(-Inf, Inf), alpha = 1,   scale = 1) == c(0, 0)
+     dgumbel(1,            alpha = Inf, scale = 1) == 0
+ })
> assertWarning(stopifnot(exprs = {
+     is.nan(dgumbel(Inf,  alpha = Inf,  scale =  1))
+     is.nan(dgumbel(-Inf, alpha = -Inf, scale =  1))
+     is.nan(dgumbel(Inf,  alpha = 1,    scale = -1))
+     is.nan(dgumbel(1,    alpha = 1,    scale = -1))
+     is.nan(dgumbel(1,    alpha = Inf,  scale = -1))
+ }))
> x <- rgumbel(100, alpha = 2, scale = 5)
> for(a in round(rlnorm(30), 2))
+ {
+     Ga <- gamma(a)
+     for(s in round(rlnorm(30), 2))
+     {
+         u <- (x - a)/s
+ 	stopifnot(exprs = {
+             All.eq(dgumbel(x, alpha = a, scale = s),
+                    exp(-(u + exp(-u)))/s)
+         })
+     }
+ }
> 
> ## Tests on the cumulative distribution function.
> assertWarning(stopifnot(exprs = {
+     is.nan(pgumbel(Inf,  alpha = Inf,  scale =  1))
+     is.nan(pgumbel(-Inf, alpha = -Inf, scale =  1))
+     is.nan(pgumbel(Inf,  alpha = 1,    scale = -1))
+     is.nan(pgumbel(1,    alpha = 1,    scale = -1))
+     is.nan(pgumbel(1,    alpha = Inf,  scale = -1))
+ }))
> scLrg <- c(2, 100, 1e300 * c(0.1, 1, 10, 100), 1e307, xMax, Inf)
> stopifnot(exprs = {
+     pgumbel(c(-Inf, Inf),  2, scale = xMax) == c(0, 1)
+     pgumbel(c(xMin, xMax), 2, scale = 0) == c(0, 1)
+     all.equal(pgumbel(1e300, 0, scale = scLrg, log = TRUE),
+               -exp(-c(5e299, 1e+298, 10, 1, 0.1, 0.01, 1e-7, 1e+300/xMax, 0)))
+ })
> 
> ## Test the first two moments, the only ones implemented.
> assertWarning(stopifnot(exprs = {
+     is.nan(mgumbel(c(-2, -1, 3, 4), alpha = 2, scale = 5))
+ }))
> stopifnot(exprs = {
+     All.eq(mgumbel(1, alpha = 2, scale = 5),
+            2 + 5 * 0.577215664901532860606512090082)
+     All.eq(mgumbel(2, alpha = 2, scale = 5),
+            pi^2 * 25/6 + (2 + 5 * 0.577215664901532860606512090082)^2)
+ })
> 
> ## INVERSE GAUSSIAN
> 
> ## Tests on the density.
> stopifnot(exprs = {
+     dinvgauss(c(1, 3, Inf),  mean = 2,   dispersion = Inf) == c(0, 0, 0)
+     dinvgauss(c(0, 42, Inf), mean = 2,   dispersion = 0) == c(Inf, 0, 0)
+     dinvgauss(c(0, Inf),     mean = 1,   dispersion = 1) == c(0, 0)
+     dinvgauss(1,             mean = Inf, dispersion = 2) == dinvgamma(1, 0.5, scale = 0.25)
+ })
> assertWarning(stopifnot(exprs = {
+     is.nan(dinvgauss(-Inf, mean = -1,  dispersion =  1))
+     is.nan(dinvgauss(Inf,  mean = 1,   dispersion = -1))
+     is.nan(dinvgauss(1,    mean = 1,   dispersion = -1))
+     is.nan(dinvgauss(1,    mean = Inf, dispersion = -1))
+ }))
> x <- rinvgauss(100, mean = 2, dispersion = 5)
> for(mu in round(rlnorm(30), 2))
+ {
+     for(phi in round(rlnorm(30), 2))
+ 	stopifnot(exprs = {
+             All.eq(dinvgauss(x, mean = mu, dispersion = phi),
+                    1/sqrt(2*pi*phi*x^3) * exp(-((x/mu - 1)^2)/(2*phi*x)))
+         })
+ }
> 
> ## Tests on the cumulative distribution function.
> assertWarning(stopifnot(exprs = {
+     is.nan(pinvgauss(-Inf, mean = -Inf, dispersion =  1))
+     is.nan(pinvgauss(Inf,  mean = 1,    dispersion = -1))
+     is.nan(pinvgauss(1,    mean = Inf,  dispersion = -1))
+ }))
> x <- c(1:50, 10^c(3:10, 20, 50, 150, 250))
> sqx <- sqrt(x)
> stopifnot(exprs = {
+     pinvgauss(c(0, Inf),  mean = 2, dispersion = xMax) == c(0, 1)
+     pinvgauss(c(0, xMax), mean = xMax, dispersion = 0) == c(0, 1)
+     all.equal(pinvgauss(x, 1, dispersion = 1, log = TRUE),
+               log(pnorm(sqx - 1/sqx) + exp(2) * pnorm(-sqx - 1/sqx)))
+ })
> 
> ## Tests for small value of 'shape'. Added for the patch in 4294e9c.
> q <- runif(100)
> stopifnot(exprs = {
+     all.equal(q,
+               pinvgauss(qinvgauss(q, 0.1, 1e-2), 0.1, 1e-2))
+     all.equal(q,
+               pinvgauss(qinvgauss(q, 0.1, 1e-6), 0.1, 1e-6))
+ })
> 
> ## Tests for first three positive, integer moments.
> k <- 1:3
> for(mu in round(rlnorm(30), 2))
+ {
+     for(phi in round(rlnorm(30), 2))
+ 	stopifnot(exprs = {
+             All.eq(minvgauss(k, mean = mu, dispersion = phi),
+                    c(mu,
+                      mu^2 * (1 + phi * mu),
+                      mu^3 * (1 + 3 * phi * mu + 3 * (phi * mu)^2)))
+         })
+ }
> 
> ## Tests for limited expected value.
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)     # quantiles
> for(mu in round(rlnorm(30), 2))
+ {
+     for(phi in round(rlnorm(30), 2))
+     {
+         limit <- qinvgauss(q, mean = mu, dispersion = phi)
+         stopifnot(exprs = {
+             All.eq(levinvgauss(limit, mean = mu, dispersion = phi),
+                    mu * (pnorm((limit/mu - 1)/sqrt(phi * limit)) -
+                          exp(2/phi/mu) * pnorm(-(limit/mu + 1)/sqrt(phi * limit))) +
+                    limit * pinvgauss(limit, mean = mu, dispersion = phi, lower = FALSE))
+         })
+     }
+ }
> 
> ## GENERALIZED BETA
> stopifnot(exprs = {
+     dgenbeta(c(0, 2.5, 5), shape1 = 0,   shape2 = 0,   shape3 = 3,   scale = 5) == c(Inf, 0, Inf)
+     dgenbeta(c(0, 2.5, 5), shape1 = 0,   shape2 = 0,   shape3 = 0,   scale = 5) == c(Inf, 0, Inf)
+     dgenbeta(c(0, 2.5, 5), shape1 = 0,   shape2 = 2,   shape3 = 0,   scale = 5) == c(Inf, 0, 0)
+     dgenbeta(c(0, 2.5, 5), shape1 = 0,   shape2 = Inf, shape3 = 3,   scale = 5) == c(Inf, 0, 0)
+     dgenbeta(c(0, 2.5, 5), shape1 = 1,   shape2 = Inf, shape3 = 3,   scale = 5) == c(Inf, 0, 0)
+     dgenbeta(c(0, 2.5, 5), shape1 = Inf, shape2 = Inf, shape3 = 3,   scale = 5) == c(0, Inf, 0)
+     dgenbeta(c(0, 2.5, 5), shape1 = Inf, shape2 = Inf, shape3 = Inf, scale = 5) == c(0, 0, Inf)
+ })
> nshpar <- 3                     # number of shape parameters
> shpar <- replicate(30, rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; b <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, b)
+     for (s in scpar)
+     {
+         u <- rbeta(100, a, b)
+         y <- u^(1/t)
+         x <- s * y
+         stopifnot(exprs = {
+             all.equal(d1 <- dgenbeta(x, shape1 = a, shape2 = b, shape3 = t,
+                                      scale = s),
+                       d2 <- dgenbeta(y, shape1 = a, shape2 = b, shape3 = t,
+                                      scale = 1)/s,
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * y^(a*t - 1) * (1 - y^t)^(b - 1)/(s * Be),
+                       tolerance = 1e-10)
+             all.equal(d1,
+                       t * u^a * (1 - u)^(b - 1)/(x * Be),
+                       tolerance = 1e-10)
+         })
+     }
+ }
> 
> ## Tests on the cumulative distribution function.
> scLrg <- 1e300 * c(0.5, 1, 2, 4)
> stopifnot(exprs = {
+     all.equal(pgenbeta(1e300,
+                        shape1 = 3, shape2 = 1,
+                        shape3 = rep(c(1, 2), each = length(scLrg)),
+                        scale = scLrg, log = TRUE),
+               c(0, pbeta(c(1, 1/2, 1/4), 3, 1, log = TRUE),
+                 0, pbeta(c(1, 1/4, 1/16), 3, 1, log = TRUE)))
+ })
> 
> ## Tests for first three positive moments and first two negative
> ## moments.
> ##
> ## Simulation of new parameters ensuring that said moments exist.
> set.seed(123)                   # reset the seed
> nshpar <- 3                     # number of shape parameters
> shpar <- replicate(30, sqrt(c(3, 0, 3)) + rlnorm(nshpar, 2), simplify = FALSE)
> scpar <- rlnorm(30, 2)          # scale parameters
> k <- c(-2, -1, 1, 2, 3)         # orders
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; b <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, b)
+     for (s in scpar)
+         stopifnot(exprs = {
+             All.eq(mgenbeta(k, shape1 = a, shape2 = b, shape3 = t, scale = s),
+                    s^k * beta(a + k/t, b)/Be)
+         })
+ }
> 
> ## Tests for first three positive limited moments and first two
> ## negative limited moments.
> ##
> ## Simulation of new parameters ensuring that said moments exist.
> order <- c(-2, -1, 1, 2, 3)             # orders
> q <- c(0.25, 0.50, 0.75, 0.9, 0.95)     # quantiles
> for (i in seq_along(shpar))
+ {
+     a <- shpar[[c(i, 1)]]; g <- shpar[[c(i, 2)]]; t <- shpar[[c(i, 3)]]
+     Be <- beta(a, b)
+     for (s in scpar)
+     {
+         limit <- qgenbeta(q, shape1 = a, shape2 = b, shape3 = t, scale = s)
+         u <- (limit/s)^t
+         for (k in order)
+             stopifnot(exprs = {
+                 All.eq(levgenbeta(limit, order = k, shape1 = a, shape2 = b, shape3 = t, scale = s),
+                        s^k * beta(a + k/t, b)/Be * pbeta(u, a + k/t, b) +
+                        limit^k * pbeta(u, a, b, lower = FALSE))
+             })
+     }
+ }
> 
> ##
> ## RANDOM NUMBERS (all continuous distributions)
> ##
> set.seed(123)
> n <- 20
> m <- rnorm(1)
> 
> ## Generate variates
> Rfpareto <- rfpareto(n, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Rpareto4 <- rpareto4(n, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Rpareto3 <- rpareto3(n, min = m,               shape  = 1.5,             scale = 2)
> Rpareto2 <- rpareto2(n, min = m, shape  = 0.8,                           scale = 2)
> Rtrbeta        <- rtrbeta      (n, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Rburr          <- rburr        (n, shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Rllogis        <- rllogis      (n, shape  = 1.5, scale = 2)
> Rparalogis     <- rparalogis   (n, shape  = 0.8, scale = 2)
> Rgenpareto     <- rgenpareto   (n, shape1 = 0.8, shape2 = 2, scale = 2)
> Rpareto        <- rpareto      (n, shape  = 0.8, scale = 2)
> Rpareto1       <- rpareto1     (n, shape  = 0.8, min = 2)
> Rinvburr       <- rinvburr     (n, shape1 = 1.5, shape2 = 2, scale = 2)
> Rinvpareto     <- rinvpareto   (n, shape  = 2,   scale = 2)
> Rinvparalogis  <- rinvparalogis(n, shape  = 2,   scale = 2)
> Rtrgamma       <- rtrgamma     (n, shape1 = 2, shape2 = 3, scale = 5)
> Rinvtrgamma    <- rinvtrgamma  (n, shape1 = 2, shape2 = 3, scale = 5)
> Rinvgamma      <- rinvgamma    (n, shape  = 2,             scale = 5)
> Rinvweibull    <- rinvweibull  (n,             shape  = 3, scale = 5)
> Rinvexp        <- rinvexp      (n,                         scale = 5)
> Rlgamma <- rlgamma(n, shapelog = 1.5, ratelog = 5)
> Rgumbel <- rgumbel(n, alpha = 2, scale = 5)
> Rinvgauss <- rinvgauss(n, mean = 2, dispersion = 5)
> Rgenbeta <- rgenbeta(n, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> 
> ## Compute quantiles
> Pfpareto <- pfpareto(Rfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Ppareto4 <- ppareto4(Rpareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Ppareto3 <- ppareto3(Rpareto3, min = m,               shape  = 1.5,             scale = 2)
> Ppareto2 <- ppareto2(Rpareto2, min = m, shape  = 0.8,                           scale = 2)
> Ptrbeta        <- ptrbeta      (Rtrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Pburr          <- pburr        (Rburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Pllogis        <- pllogis      (Rllogis,       shape  = 1.5, scale = 2)
> Pparalogis     <- pparalogis   (Rparalogis,    shape  = 0.8, scale = 2)
> Pgenpareto     <- pgenpareto   (Rgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2)
> Ppareto        <- ppareto      (Rpareto,       shape  = 0.8, scale = 2)
> Ppareto1       <- ppareto1     (Rpareto1,      shape  = 0.8, min = 2)
> Pinvburr       <- pinvburr     (Rinvburr,      shape1 = 1.5, shape2 = 2, scale = 2)
> Pinvpareto     <- pinvpareto   (Rinvpareto,    shape  = 2,   scale = 2)
> Pinvparalogis  <- pinvparalogis(Rinvparalogis, shape  = 2,   scale = 2)
> Ptrgamma       <- ptrgamma     (Rtrgamma,      shape1 = 2, shape2 = 3, scale = 5)
> Pinvtrgamma    <- pinvtrgamma  (Rinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5)
> Pinvgamma      <- pinvgamma    (Rinvgamma,     shape  = 2,             scale = 5)
> Pinvweibull    <- pinvweibull  (Rinvweibull,               shape  = 3, scale = 5)
> Pinvexp        <- pinvexp      (Rinvexp,                               scale = 5)
> Plgamma <- plgamma(Rlgamma, shapelog = 1.5, ratelog = 5)
> Pgumbel <- pgumbel(Rgumbel, alpha = 2, scale = 5)
> Pinvgauss <- pinvgauss(Rinvgauss, mean = 2, dispersion = 5)
> Pgenbeta <- pgenbeta(Rgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> 
> ## Just compute pdf
> Dfpareto <- dfpareto(Rfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Dpareto4 <- dpareto4(Rpareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Dpareto3 <- dpareto3(Rpareto3, min = m,               shape  = 1.5,             scale = 2)
> Dpareto2 <- dpareto2(Rpareto2, min = m, shape  = 0.8,                           scale = 2)
> Dtrbeta        <- dtrbeta      (Rtrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> Dburr          <- dburr        (Rburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2)
> Dllogis        <- dllogis      (Rllogis,       shape  = 1.5, scale = 2)
> Dparalogis     <- dparalogis   (Rparalogis,    shape  = 0.8, scale = 2)
> Dgenpareto     <- dgenpareto   (Rgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2)
> Dpareto        <- dpareto      (Rpareto,       shape  = 0.8, scale = 2)
> Dpareto1       <- dpareto1     (Rpareto1,      shape  = 0.8, min = 2)
> Dinvburr       <- dinvburr     (Rinvburr,      shape1 = 1.5, shape2 = 2, scale = 2)
> Dinvpareto     <- dinvpareto   (Rinvpareto,    shape  = 2,   scale = 2)
> Dinvparalogis  <- dinvparalogis(Rinvparalogis, shape  = 2,   scale = 2)
> Dtrgamma       <- dtrgamma     (Rtrgamma,      shape1 = 2, shape2 = 3, scale = 5)
> Dinvtrgamma    <- dinvtrgamma  (Rinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5)
> Dinvgamma      <- dinvgamma    (Rinvtrgamma,   shape  = 2,             scale = 5)
> Dinvweibull    <- dinvweibull  (Rinvweibull,               shape  = 3, scale = 5)
> Dinvexp        <- dinvexp      (Rinvexp,                               scale = 5)
> Dlgamma <- dlgamma(Rlgamma, shapelog = 1.5, ratelog = 5)
> Dgumbel <- dgumbel(Rgumbel, alpha = 2, scale = 5)
> Dinvgauss <- dinvgauss(Rinvgauss, mean = 2, dispersion = 5)
> Dgenbeta <- dgenbeta(Rgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
> 
> ## Check q<dist>(p<dist>(.)) identity
> stopifnot(exprs = {
+     All.eq(Rfpareto, qfpareto(Pfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
+     All.eq(Rpareto4, qpareto4(Ppareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2))
+     All.eq(Rpareto3, qpareto3(Ppareto3, min = m,               shape  = 1.5,             scale = 2))
+     All.eq(Rpareto2, qpareto2(Ppareto2, min = m, shape  = 0.8,                           scale = 2))
+     All.eq(Rtrbeta,       qtrbeta      (Ptrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
+     All.eq(Rburr,         qburr        (Pburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2))
+     All.eq(Rllogis,       qllogis      (Pllogis,       shape  = 1.5, scale = 2))
+     All.eq(Rparalogis,    qparalogis   (Pparalogis,    shape  = 0.8, scale = 2))
+     All.eq(Rgenpareto,    qgenpareto   (Pgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2))
+     All.eq(Rpareto,       qpareto      (Ppareto,       shape  = 0.8, scale = 2))
+     All.eq(Rpareto1,      qpareto1     (Ppareto1,      shape  = 0.8, min = 2))
+     All.eq(Rinvburr,      qinvburr     (Pinvburr,      shape1 = 1.5, shape2 = 2, scale = 2))
+     All.eq(Rinvpareto,    qinvpareto   (Pinvpareto,    shape  = 2,   scale = 2))
+     All.eq(Rinvparalogis, qinvparalogis(Pinvparalogis, shape  = 2,   scale = 2))
+     All.eq(Rtrgamma,      qtrgamma     (Ptrgamma,      shape1 = 2, shape2 = 3, scale = 5))
+     All.eq(Rinvtrgamma,   qinvtrgamma  (Pinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5))
+     All.eq(Rinvgamma,     qinvgamma    (Pinvgamma,     shape  = 2,             scale = 5))
+     All.eq(Rinvweibull,   qinvweibull  (Pinvweibull,               shape  = 3, scale = 5))
+     All.eq(Rinvexp,       qinvexp      (Pinvexp,                               scale = 5))
+     All.eq(Rlgamma, qlgamma(Plgamma, shapelog = 1.5, ratelog = 5))
+     All.eq(Rgumbel, qgumbel(Pgumbel, alpha = 2, scale = 5))
+     All.eq(Rinvgauss, qinvgauss(Pinvgauss, mean = 2, dispersion = 5))
+     All.eq(Rgenbeta, qgenbeta(Pgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
+ })
> 
> ## Check q<dist>(p<dist>(.)) identity for special cases
> stopifnot(exprs = {
+     All.eq(Rfpareto - m, qtrbeta(Pfpareto, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
+     All.eq(Rpareto4 - m, qburr  (Ppareto4, shape1 = 0.8, shape2 = 1.5,             scale = 2))
+     All.eq(Rpareto3 - m, qllogis(Ppareto3,               shape  = 1.5,             scale = 2))
+     All.eq(Rpareto2 - m, qpareto(Ppareto2, shape  = 0.8,                           scale = 2))
+ })
> 
> ## Check q<dist>(p<dist>(.)) identity with upper tail
> stopifnot(exprs = {
+     All.eq(Rfpareto, qfpareto(1 - Pfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE))
+     All.eq(Rpareto4, qpareto4(1 - Ppareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2, lower = FALSE))
+     All.eq(Rpareto3, qpareto3(1 - Ppareto3, min = m,               shape  = 1.5,             scale = 2, lower = FALSE))
+     All.eq(Rpareto2, qpareto2(1 - Ppareto2, min = m, shape  = 0.8,                           scale = 2, lower = FALSE))
+     All.eq(Rtrbeta,       qtrbeta      (1 - Ptrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE))
+     All.eq(Rburr,         qburr        (1 - Pburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2, lower = FALSE))
+     All.eq(Rllogis,       qllogis      (1 - Pllogis,       shape  = 1.5, scale = 2, lower = FALSE))
+     All.eq(Rparalogis,    qparalogis   (1 - Pparalogis,    shape  = 0.8, scale = 2, lower = FALSE))
+     All.eq(Rgenpareto,    qgenpareto   (1 - Pgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2, lower = FALSE))
+     All.eq(Rpareto,       qpareto      (1 - Ppareto,       shape  = 0.8, scale = 2, lower = FALSE))
+     All.eq(Rpareto1,      qpareto1     (1 - Ppareto1,      shape  = 0.8, min = 2, lower = FALSE))
+     All.eq(Rinvburr,      qinvburr     (1 - Pinvburr,      shape1 = 1.5, shape2 = 2, scale = 2, lower = FALSE))
+     All.eq(Rinvpareto,    qinvpareto   (1 - Pinvpareto,    shape  = 2,   scale = 2, lower = FALSE))
+     All.eq(Rinvparalogis, qinvparalogis(1 - Pinvparalogis, shape  = 2,   scale = 2, lower = FALSE))
+     All.eq(Rtrgamma,      qtrgamma     (1 - Ptrgamma,      shape1 = 2, shape2 = 3, scale = 5, lower = FALSE))
+     All.eq(Rinvtrgamma,   qinvtrgamma  (1 - Pinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5, lower = FALSE))
+     All.eq(Rinvgamma,     qinvgamma    (1 - Pinvgamma,     shape  = 2,             scale = 5, lower = FALSE))
+     All.eq(Rinvweibull,   qinvweibull  (1 - Pinvweibull,               shape  = 3, scale = 5, lower = FALSE))
+     All.eq(Rinvexp,       qinvexp      (1 - Pinvexp,                               scale = 5, lower = FALSE))
+     All.eq(Rlgamma, qlgamma(1 - Plgamma, shapelog = 1.5, ratelog = 5, lower = FALSE))
+     All.eq(Rgumbel, qgumbel(1 - Pgumbel, alpha = 2, scale = 5, lower = FALSE))
+     All.eq(Rinvgauss, qinvgauss(1 - Pinvgauss, mean = 2, dispersion = 5, lower = FALSE))
+     All.eq(Rgenbeta, qgenbeta(1 - Pgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE))
+ })
> 
> ## Check q<dist>(p<dist>(., log), log) identity
> stopifnot(exprs = {
+     All.eq(Rfpareto, qfpareto(log(Pfpareto), min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, log = TRUE))
+     All.eq(Rpareto4, qpareto4(log(Ppareto4), min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2, log = TRUE))
+     All.eq(Rpareto3, qpareto3(log(Ppareto3), min = m,               shape  = 1.5,             scale = 2, log = TRUE))
+     All.eq(Rpareto2, qpareto2(log(Ppareto2), min = m, shape  = 0.8,                           scale = 2, log = TRUE))
+     All.eq(Rtrbeta,       qtrbeta      (log(Ptrbeta),       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, log = TRUE))
+     All.eq(Rburr,         qburr        (log(Pburr),         shape1 = 0.8, shape2 = 1.5,             scale = 2, log = TRUE))
+     All.eq(Rllogis,       qllogis      (log(Pllogis),       shape  = 1.5, scale = 2, log = TRUE))
+     All.eq(Rparalogis,    qparalogis   (log(Pparalogis),    shape  = 0.8, scale = 2, log = TRUE))
+     All.eq(Rgenpareto,    qgenpareto   (log(Pgenpareto),    shape1 = 0.8, shape2 = 2, scale = 2, log = TRUE))
+     All.eq(Rpareto,       qpareto      (log(Ppareto),       shape  = 0.8, scale = 2, log = TRUE))
+     All.eq(Rpareto1,      qpareto1     (log(Ppareto1),      shape  = 0.8, min = 2, log = TRUE))
+     All.eq(Rinvburr,      qinvburr     (log(Pinvburr),      shape1 = 1.5, shape2 = 2, scale = 2, log = TRUE))
+     All.eq(Rinvpareto,    qinvpareto   (log(Pinvpareto),    shape  = 2,   scale = 2, log = TRUE))
+     All.eq(Rinvparalogis, qinvparalogis(log(Pinvparalogis), shape  = 2,   scale = 2, log = TRUE))
+     All.eq(Rtrgamma,      qtrgamma     (log(Ptrgamma),      shape1 = 2, shape2 = 3, scale = 5, log = TRUE))
+     All.eq(Rinvtrgamma,   qinvtrgamma  (log(Pinvtrgamma),   shape1 = 2, shape2 = 3, scale = 5, log = TRUE))
+     All.eq(Rinvgamma,     qinvgamma    (log(Pinvgamma),     shape  = 2,             scale = 5, log = TRUE))
+     All.eq(Rinvweibull,   qinvweibull  (log(Pinvweibull),               shape  = 3, scale = 5, log = TRUE))
+     All.eq(Rinvexp,       qinvexp      (log(Pinvexp),                               scale = 5, log = TRUE))
+     All.eq(Rlgamma, qlgamma(log(Plgamma), shapelog = 1.5, ratelog = 5, log = TRUE))
+     All.eq(Rgumbel, qgumbel(log(Pgumbel), alpha = 2, scale = 5, log = TRUE))
+     All.eq(Rinvgauss, qinvgauss(log(Pinvgauss), mean = 2, dispersion = 5, log = TRUE))
+     All.eq(Rgenbeta, qgenbeta(log(Pgenbeta), shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, log = TRUE))
+ })
> 
> ## Check q<dist>(p<dist>(., log), log) identity with upper tail
> stopifnot(exprs = {
+     All.eq(Rfpareto, qfpareto(log1p(-Pfpareto), min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rpareto4, qpareto4(log1p(-Ppareto4), min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rpareto3, qpareto3(log1p(-Ppareto3), min = m,               shape  = 1.5,             scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rpareto2, qpareto2(log1p(-Ppareto2), min = m, shape  = 0.8,                           scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rtrbeta,       qtrbeta      (log1p(-Ptrbeta),       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rburr,         qburr        (log1p(-Pburr),         shape1 = 0.8, shape2 = 1.5,             scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rllogis,       qllogis      (log1p(-Pllogis),       shape  = 1.5, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rparalogis,    qparalogis   (log1p(-Pparalogis),    shape  = 0.8, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rgenpareto,    qgenpareto   (log1p(-Pgenpareto),    shape1 = 0.8, shape2 = 2, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rpareto,       qpareto      (log1p(-Ppareto),       shape  = 0.8, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rpareto1,      qpareto1     (log1p(-Ppareto1),      shape  = 0.8, min = 2, lower = FALSE, log = TRUE))
+     All.eq(Rinvburr,      qinvburr     (log1p(-Pinvburr),      shape1 = 1.5, shape2 = 2, scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rinvpareto,    qinvpareto   (log1p(-Pinvpareto),    shape  = 2,   scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rinvparalogis, qinvparalogis(log1p(-Pinvparalogis), shape  = 2,   scale = 2, lower = FALSE, log = TRUE))
+     All.eq(Rtrgamma,      qtrgamma     (log1p(-Ptrgamma),      shape1 = 2, shape2 = 3, scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rinvtrgamma,   qinvtrgamma  (log1p(-Pinvtrgamma),   shape1 = 2, shape2 = 3, scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rinvgamma,     qinvgamma    (log1p(-Pinvgamma),     shape  = 2,             scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rinvweibull,   qinvweibull  (log1p(-Pinvweibull),               shape  = 3, scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rinvexp,       qinvexp      (log1p(-Pinvexp),                               scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rlgamma, qlgamma(log1p(-Plgamma), shapelog = 1.5, ratelog = 5, lower = FALSE, log = TRUE))
+     All.eq(Rgumbel, qgumbel(log1p(-Pgumbel), alpha = 2, scale = 5, lower = FALSE, log = TRUE))
+     All.eq(Rinvgauss, qinvgauss(log1p(-Pinvgauss), mean = 2, dispersion = 5, lower = FALSE, log = TRUE))
+     All.eq(Rgenbeta, qgenbeta(log1p(-Pgenbeta), shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2, lower = FALSE, log = TRUE))
+ })
> 
> 
> ###
> ### DISCRETE DISTRIBUTIONS
> ###
> 
> ## Reset seed
> set.seed(123)
> 
> ## Define a small function to compute probabilities for the (a, b, 1)
> ## family of discrete distributions using the recursive relation
> ##
> ##   p[k] = (a + b/k)p[k - 1], k = 2, 3, ...
> ##
> ## for a, b and p[1] given.
> dab1 <- function(x, a, b, p1)
+ {
+     x <- floor(x)
+     if (x < 1)
+         stop("recursive computations possible for x >= 2 only")
+     for (k in seq(2, length.out = x - 1))
+     {
+         p2 <- (a + b/k) * p1
+         p1 <- p2
+     }
+     p1
+ }
> 
> ## ZERO-TRUNCATED (a, b, 1) CLASS
> 
> ## Tests on the probability mass function:
> ##
> ## 1. probability is 0 at x = 0;
> ## 2. pmf satisfies the recursive relation
> lambda <- rlnorm(30, 2)                 # Poisson parameters
> r <- lambda                             # size for negative binomial
> prob <- runif(30)                       # probs
> size <- round(lambda)                   # size for binomial
> stopifnot(exprs = {
+     dztpois(0, lambda) == 0
+     dztnbinom(0, r, prob) == 0
+     dztgeom(0, prob) == 0
+     dztbinom(0, size, prob) == 0
+     dlogarithmic(0, prob) == 0
+ })
> 
> x <- sapply(size, sample, size = 1)
> stopifnot(exprs = {
+     All.eq(dztpois(x, lambda),
+            mapply(dab1, x,
+                   a = 0,
+                   b = lambda,
+                   p1 = lambda/(exp(lambda) - 1)))
+     All.eq(dztnbinom(x, r, prob),
+            mapply(dab1, x,
+                   a = 1 - prob,
+                   b = (r - 1) * (1 - prob),
+                   p1 = r * prob^r * (1 - prob)/(1 - prob^r)))
+     All.eq(dztgeom(x, prob),
+            mapply(dab1, x,
+                   a = 1 - prob,
+                   b = 0,
+                   p1 = prob))
+     All.eq(dztbinom(x, size, prob),
+            mapply(dab1, x,
+                   a = -prob/(1 - prob),
+                   b = (size + 1) * prob/(1 - prob),
+                   p1 = size * prob * (1 - prob)^(size - 1)/(1 - (1 - prob)^size)))
+     All.eq(dlogarithmic(x, prob),
+            mapply(dab1, x,
+                   a = prob,
+                   b = -prob,
+                   p1 = -prob/log1p(-prob)))
+ })
> 
> ## Tests on cumulative distribution function.
> for (l in lambda)
+     stopifnot(exprs = {
+         all.equal(cumsum(dztpois(0:20, l)),
+                   pztpois(0:20, l),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(r))
+     stopifnot(exprs = {
+         all.equal(cumsum(dztnbinom(0:20, r[i], prob[i])),
+                   pztnbinom(0:20, r[i], prob[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(r))
+     stopifnot(exprs = {
+         all.equal(cumsum(dztgeom(0:20, prob[i])),
+                   pztgeom(0:20, prob[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(size))
+     stopifnot(exprs = {
+         all.equal(cumsum(dztbinom(0:20, size[i], prob[i])),
+                   pztbinom(0:20, size[i], prob[i]),
+                   tolerance = 1e-8)
+     })
> for (p in prob)
+     stopifnot(exprs = {
+         all.equal(cumsum(dlogarithmic(0:20, p)),
+                   plogarithmic(0:20, p),
+                   tolerance = 1e-8)
+     })
> 
> ## ZERO-MODIFIED (a, b, 1) CLASS
> 
> ## Tests on the probability mass function:
> ##
> ## 1. probability is p0 at x = 0 (trivial, but...);
> ## 2. pmf satisfies the recursive relation
> lambda <- rlnorm(30, 2)                 # Poisson parameters
> r <- lambda                             # size for negative binomial
> prob <- runif(30)                       # probs
> size <- round(lambda)                   # size for binomial
> p0 <- runif(30)                         # probs at 0
> stopifnot(exprs = {
+     dzmpois(0, lambda, p0) == p0
+     dzmnbinom(0, r, prob, p0) == p0
+     dzmgeom(0, prob, p0) == p0
+     dzmbinom(0, size, prob, p0) == p0
+     dzmlogarithmic(0, prob, p0) == p0
+ })
> 
> x <- sapply(size, sample, size = 1)
> stopifnot(exprs = {
+     All.eq(dzmpois(x, lambda, p0),
+            mapply(dab1, x,
+                   a = 0,
+                   b = lambda,
+                   p1 = (1 - p0) *lambda/(exp(lambda) - 1)))
+     All.eq(dzmnbinom(x, r, prob, p0),
+            mapply(dab1, x,
+                   a = 1 - prob,
+                   b = (r - 1) * (1 - prob),
+                   p1 = (1 - p0) * r * prob^r * (1 - prob)/(1 - prob^r)))
+     All.eq(dzmgeom(x, prob, p0),
+            mapply(dab1, x,
+                   a = 1 - prob,
+                   b = 0,
+                   p1 = (1 - p0) * prob))
+     All.eq(dzmbinom(x, size, prob, p0),
+            mapply(dab1, x,
+                   a = -prob/(1 - prob),
+                   b = (size + 1) * prob/(1 - prob),
+                   p1 = (1 - p0) * size * prob * (1 - prob)^(size - 1)/(1 - (1 - prob)^size)))
+     All.eq(dzmlogarithmic(x, prob, p0),
+            mapply(dab1, x,
+                   a = prob,
+                   b = -prob,
+                   p1 = -(1 - p0) * prob/log1p(-prob)))
+ })
> 
> ## Tests on cumulative distribution function.
> for (i in seq_along(lambda))
+     stopifnot(exprs = {
+         all.equal(cumsum(dzmpois(0:20, lambda[i], p0 = p0[i])),
+                   pzmpois(0:20, lambda[i], p0 = p0[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(r))
+     stopifnot(exprs = {
+         all.equal(cumsum(dzmnbinom(0:20, r[i], prob[i], p0[i])),
+                   pzmnbinom(0:20, r[i], prob[i], p0[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(r))
+     stopifnot(exprs = {
+         all.equal(cumsum(dzmgeom(0:20, prob[i], p0[i])),
+                   pzmgeom(0:20, prob[i], p0[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(size))
+     stopifnot(exprs = {
+         all.equal(cumsum(dzmbinom(0:20, size[i], prob[i], p0[i])),
+                   pzmbinom(0:20, size[i], prob[i], p0[i]),
+                   tolerance = 1e-8)
+     })
> for (i in seq_along(prob))
+     stopifnot(exprs = {
+         all.equal(cumsum(dzmlogarithmic(0:20, prob[i], p0[i])),
+                   pzmlogarithmic(0:20, prob[i], p0[i]),
+                   tolerance = 1e-8)
+     })
> 
> ## POISSON-INVERSE GAUSSIAN
> 
> ## Reset seed
> set.seed(123)
> 
> ## Define a small function to compute probabilities for the PIG
> ## directly using the Bessel function.
> dpigBK <- function(x, mu, phi)
+ {
+     M_LN2 <- 0.693147180559945309417232121458
+     M_SQRT_2dPI <- 0.225791352644727432363097614947
+ 
+     phimu  <- phi * mu
+     lphi <- log(phi)
+     y <- x - 0.5
+ 
+     logA = -lphi/2 - M_SQRT_2dPI
+     logB = (M_LN2 + lphi + log1p(1/(2 * phimu * mu)))/2;
+ 
+     exp(logA + 1/phimu - lfactorial(x) - y * logB) *
+         besselK(exp(logB - lphi), y)
+ }
> 
> ## Tests on the probability mass function.
> mu <- rlnorm(30, 2)
> phi <- rlnorm(30, 2)
> x <- 0:100
> for (i in seq_along(phi))
+ {
+     stopifnot(exprs = {
+         all.equal(dpoisinvgauss(x, mean = mu[i], dispersion = phi[i]),
+                   dpigBK(x, mu[i], phi[i]))
+         all.equal(dpoisinvgauss(x, mean = Inf, dispersion = phi[i]),
+                   dpigBK(x, Inf, phi[i]))
+     })
+ }
> 
> ## Tests on cumulative distribution function.
> for (i in seq_along(phi))
+     stopifnot(exprs = {
+         all.equal(cumsum(dpoisinvgauss(0:20, mu[i], phi[i])),
+                   ppoisinvgauss(0:20, mu[i], phi[i]),
+                   tolerance = 1e-8)
+         all.equal(cumsum(dpoisinvgauss(0:20, Inf, phi[i])),
+                   ppoisinvgauss(0:20, Inf, phi[i]),
+                   tolerance = 1e-8)
+     })
> 
> ##
> ## RANDOM NUMBERS (all discrete distributions)
> ##
> set.seed(123)
> n <- 20
> 
> ## Generate variates.
> ##
> ## For zero-modified distributions, we simulate two sets of values:
> ## one with p0m < p0 (suffix 'p0lt') and one with p0m > p0 (suffix
> ## 'p0gt').
> Rztpois      <- rztpois     (n, lambda = 12)
> Rztnbinom    <- rztnbinom   (n, size = 7, prob = 0.01)
> Rztgeom      <- rztgeom     (n, prob = pi/16)
> Rztbinom     <- rztbinom    (n, size = 55, prob = pi/16)
> Rlogarithmic <- rlogarithmic(n, prob = 0.99)
> Rzmpoisp0lt        <- rzmpois       (n, lambda = 6, p0 = 0.001)
> Rzmpoisp0gt        <- rzmpois       (n, lambda = 6, p0 = 0.010)
> Rzmnbinomp0lt      <- rzmnbinom     (n, size = 7, prob = 0.8, p0 = 0.01)
> Rzmnbinomp0gt      <- rzmnbinom     (n, size = 7, prob = 0.8, p0 = 0.40)
> Rzmgeomp0lt        <- rzmgeom       (n, prob = pi/16, p0 = 0.01)
> Rzmgeomp0gt        <- rzmgeom       (n, prob = pi/16, p0 = 0.40)
> Rzmbinomp0lt       <- rzmbinom      (n, size = 12, prob = pi/16, p0 = 0.01)
> Rzmbinomp0gt       <- rzmbinom      (n, size = 12, prob = pi/16, p0 = 0.12)
> Rzmlogarithmicp0lt <- rzmlogarithmic(n, prob = 0.99, p0 = 0.05)
> Rzmlogarithmicp0gt <- rzmlogarithmic(n, prob = 0.99, p0 = 0.55)
> Rpoisinvgauss    <- rpoisinvgauss(n, mean = 12,  dispersion = 0.1)
> RpoisinvgaussInf <- rpoisinvgauss(n, mean = Inf, dispersion = 1.1)
> 
> ## Compute quantiles
> Pztpois      <- pztpois     (Rztpois,      lambda = 12)
> Pztnbinom    <- pztnbinom   (Rztnbinom,    size = 7, prob = 0.01)
> Pztgeom      <- pztgeom     (Rztgeom,      prob = pi/16)
> Pztbinom     <- pztbinom    (Rztbinom,     size = 55, prob = pi/16)
> Plogarithmic <- plogarithmic(Rlogarithmic, prob = 0.99)
> Pzmpoisp0lt        <- pzmpois       (Rzmpoisp0lt,        lambda = 6, p0 = 0.001)
> Pzmpoisp0gt        <- pzmpois       (Rzmpoisp0gt,        lambda = 6, p0 = 0.010)
> Pzmnbinomp0lt      <- pzmnbinom     (Rzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01)
> Pzmnbinomp0gt      <- pzmnbinom     (Rzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40)
> Pzmgeomp0lt        <- pzmgeom       (Rzmgeomp0lt,        prob = pi/16, p0 = 0.01)
> Pzmgeomp0gt        <- pzmgeom       (Rzmgeomp0gt,        prob = pi/16, p0 = 0.40)
> Pzmbinomp0lt       <- pzmbinom      (Rzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01)
> Pzmbinomp0gt       <- pzmbinom      (Rzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12)
> Pzmlogarithmicp0lt <- pzmlogarithmic(Rzmlogarithmicp0lt, prob = 0.99, p0 = 0.05)
> Pzmlogarithmicp0gt <- pzmlogarithmic(Rzmlogarithmicp0gt, prob = 0.99, p0 = 0.55)
> Ppoisinvgauss    <- ppoisinvgauss(Rpoisinvgauss,    mean = 12,  dispersion = 0.1)
> PpoisinvgaussInf <- ppoisinvgauss(RpoisinvgaussInf, mean = Inf, dispersion = 1.1)
> 
> ## Just compute pmf
> Dztpois      <- dztpois     (Rztpois,      lambda = 12)
> Dztnbinom    <- dztnbinom   (Rztnbinom,    size = 7, prob = 0.01)
> Dztgeom      <- dztgeom     (Rztgeom,      prob = pi/16)
> Dztbinom     <- dztbinom    (Rztbinom,     size = 55, prob = pi/16)
> Dlogarithmic <- dlogarithmic(Rlogarithmic, prob = pi/16)
> Dzmpoisp0lt        <- dzmpois       (Rzmpoisp0lt,        lambda = 6, p0 = 0.001)
> Dzmpoisp0gt        <- dzmpois       (Rzmpoisp0gt,        lambda = 6, p0 = 0.010)
> Dzmnbinomp0lt      <- dzmnbinom     (Rzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01)
> Dzmnbinomp0gt      <- dzmnbinom     (Rzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40)
> Dzmgeomp0lt        <- dzmgeom       (Rzmgeomp0lt,        prob = pi/16, p0 = 0.01)
> Dzmgeomp0gt        <- dzmgeom       (Rzmgeomp0gt,        prob = pi/16, p0 = 0.40)
> Dzmbinomp0lt       <- dzmbinom      (Rzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01)
> Dzmbinomp0gt       <- dzmbinom      (Rzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12)
> Dzmlogarithmicp0lt <- dzmlogarithmic(Rzmlogarithmicp0lt, prob = 0.99, p0 = 0.05)
> Dzmlogarithmicp0gt <- dzmlogarithmic(Rzmlogarithmicp0gt, prob = 0.99, p0 = 0.55)
> Dpoisinvgauss    <- dpoisinvgauss(Rpoisinvgauss,    mean = 12,  dispersion = 0.1)
> DpoisinvgaussInf <- dpoisinvgauss(RpoisinvgaussInf, mean = Inf, dispersion = 1.1)
> 
> ## Check q<dist>(p<dist>(.)) identity
> stopifnot(exprs = {
+     Rztpois      == qztpois     (Pztpois,      lambda = 12)
+     Rztnbinom    == qztnbinom   (Pztnbinom,    size = 7, prob = 0.01)
+     Rztgeom      == qztgeom     (Pztgeom,      prob = pi/16)
+     Rztbinom     == qztbinom    (Pztbinom,     size = 55, prob = pi/16)
+     Rlogarithmic == qlogarithmic(Plogarithmic, prob = 0.99)
+     Rzmpoisp0lt        == qzmpois       (Pzmpoisp0lt,        lambda = 6, p0 = 0.001)
+     Rzmpoisp0gt        == qzmpois       (Pzmpoisp0gt,        lambda = 6, p0 = 0.010)
+     Rzmnbinomp0lt      == qzmnbinom     (Pzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01)
+     Rzmnbinomp0gt      == qzmnbinom     (Pzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40)
+     Rzmgeomp0lt        == qzmgeom       (Pzmgeomp0lt,        prob = pi/16, p0 = 0.01)
+     Rzmgeomp0gt        == qzmgeom       (Pzmgeomp0gt,        prob = pi/16, p0 = 0.40)
+     Rzmbinomp0lt       == qzmbinom      (Pzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01)
+     Rzmbinomp0gt       == qzmbinom      (Pzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12)
+     Rzmlogarithmicp0lt == qzmlogarithmic(Pzmlogarithmicp0lt, prob = 0.99, p0 = 0.05)
+     Rzmlogarithmicp0gt == qzmlogarithmic(Pzmlogarithmicp0gt, prob = 0.99, p0 = 0.55)
+     Rpoisinvgauss    == qpoisinvgauss(Ppoisinvgauss,    mean = 12,  dispersion = 0.1)
+     RpoisinvgaussInf == qpoisinvgauss(PpoisinvgaussInf, mean = Inf, dispersion = 1.1)
+ })
> 
> ## Check q<dist>(p<dist>(.)) identity with upper tail
> stopifnot(exprs = {
+     Rztpois      == qztpois     (1 - Pztpois,      lambda = 12, lower = FALSE)
+     Rztnbinom    == qztnbinom   (1 - Pztnbinom,    size = 7, prob = 0.01, lower = FALSE)
+     Rztgeom      == qztgeom     (1 - Pztgeom,      prob = pi/16, lower = FALSE)
+     Rztbinom     == qztbinom    (1 - Pztbinom,     size = 55, prob = pi/16, lower = FALSE)
+     Rlogarithmic == qlogarithmic(1 - Plogarithmic, prob = 0.99, lower = FALSE)
+     Rzmpoisp0lt        == qzmpois       (1 - Pzmpoisp0lt,        lambda = 6, p0 = 0.001, lower = FALSE)
+     Rzmpoisp0gt        == qzmpois       (1 - Pzmpoisp0gt,        lambda = 6, p0 = 0.010, lower = FALSE)
+     Rzmnbinomp0lt      == qzmnbinom     (1 - Pzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01, lower = FALSE)
+     Rzmnbinomp0gt      == qzmnbinom     (1 - Pzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40, lower = FALSE)
+     Rzmgeomp0lt        == qzmgeom       (1 - Pzmgeomp0lt,        prob = pi/16, p0 = 0.01, lower = FALSE)
+     Rzmgeomp0gt        == qzmgeom       (1 - Pzmgeomp0gt,        prob = pi/16, p0 = 0.40, lower = FALSE)
+     Rzmbinomp0lt       == qzmbinom      (1 - Pzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01, lower = FALSE)
+     Rzmbinomp0gt       == qzmbinom      (1 - Pzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12, lower = FALSE)
+     Rzmlogarithmicp0lt == qzmlogarithmic(1 - Pzmlogarithmicp0lt, prob = 0.99, p0 = 0.05, lower = FALSE)
+     Rzmlogarithmicp0gt == qzmlogarithmic(1 - Pzmlogarithmicp0gt, prob = 0.99, p0 = 0.55, lower = FALSE)
+     Rpoisinvgauss    == qpoisinvgauss(1 - Ppoisinvgauss,    mean = 12,  dispersion = 0.1, lower = FALSE)
+     RpoisinvgaussInf == qpoisinvgauss(1 - PpoisinvgaussInf, mean = Inf, dispersion = 1.1, lower = FALSE)
+ })
> 
> ## Check q<dist>(p<dist>(., log), log) identity
> stopifnot(exprs = {
+     Rztpois      == qztpois     (log(Pztpois),      lambda = 12, log = TRUE)
+     Rztnbinom    == qztnbinom   (log(Pztnbinom),    size = 7, prob = 0.01, log = TRUE)
+     Rztgeom      == qztgeom     (log(Pztgeom),      prob = pi/16, log = TRUE)
+     Rztbinom     == qztbinom    (log(Pztbinom),     size = 55, prob = pi/16, log = TRUE)
+     Rlogarithmic == qlogarithmic(log(Plogarithmic), prob = 0.99, log = TRUE)
+     Rzmpoisp0lt        == qzmpois       (log(Pzmpoisp0lt),        lambda = 6, p0 = 0.001, log = TRUE)
+     Rzmpoisp0gt        == qzmpois       (log(Pzmpoisp0gt),        lambda = 6, p0 = 0.010, log = TRUE)
+     Rzmnbinomp0lt      == qzmnbinom     (log(Pzmnbinomp0lt),      size = 7, prob = 0.8, p0 = 0.01, log = TRUE)
+     Rzmnbinomp0gt      == qzmnbinom     (log(Pzmnbinomp0gt),      size = 7, prob = 0.8, p0 = 0.40, log = TRUE)
+     Rzmgeomp0lt        == qzmgeom       (log(Pzmgeomp0lt),        prob = pi/16, p0 = 0.01, log = TRUE)
+     Rzmgeomp0gt        == qzmgeom       (log(Pzmgeomp0gt),        prob = pi/16, p0 = 0.40, log = TRUE)
+     Rzmbinomp0lt       == qzmbinom      (log(Pzmbinomp0lt),       size = 12, prob = pi/16, p0 = 0.01, log = TRUE)
+     Rzmbinomp0gt       == qzmbinom      (log(Pzmbinomp0gt),       size = 12, prob = pi/16, p0 = 0.12, log = TRUE)
+     Rzmlogarithmicp0lt == qzmlogarithmic(log(Pzmlogarithmicp0lt), prob = 0.99, p0 = 0.05, log = TRUE)
+     Rzmlogarithmicp0gt == qzmlogarithmic(log(Pzmlogarithmicp0gt), prob = 0.99, p0 = 0.55, log = TRUE)
+     Rpoisinvgauss    == qpoisinvgauss(log(Ppoisinvgauss),    mean = 12,  dispersion = 0.1, log = TRUE)
+     RpoisinvgaussInf == qpoisinvgauss(log(PpoisinvgaussInf), mean = Inf, dispersion = 1.1, log = TRUE)
+ })
> 
> ## Check q<dist>(p<dist>(., log), log) identity with upper tail
> stopifnot(exprs = {
+     Rztpois      == qztpois     (log1p(-Pztpois),      lambda = 12, lower = FALSE, log = TRUE)
+     Rztnbinom    == qztnbinom   (log1p(-Pztnbinom),    size = 7, prob = 0.01, lower = FALSE, log = TRUE)
+     Rztgeom      == qztgeom     (log1p(-Pztgeom),      prob = pi/16, lower = FALSE, log = TRUE)
+     Rztbinom     == qztbinom    (log1p(-Pztbinom),     size = 55, prob = pi/16, lower = FALSE, log = TRUE)
+     Rlogarithmic == qlogarithmic(log1p(-Plogarithmic), prob = 0.99, lower = FALSE, log = TRUE)
+     Rzmpoisp0lt        == qzmpois       (log1p(-Pzmpoisp0lt),        lambda = 6, p0 = 0.001, lower = FALSE, log = TRUE)
+     Rzmpoisp0gt        == qzmpois       (log1p(-Pzmpoisp0gt),        lambda = 6, p0 = 0.010, lower = FALSE, log = TRUE)
+     Rzmnbinomp0lt      == qzmnbinom     (log1p(-Pzmnbinomp0lt),      size = 7, prob = 0.8, p0 = 0.01, lower = FALSE, log = TRUE)
+     Rzmnbinomp0gt      == qzmnbinom     (log1p(-Pzmnbinomp0gt),      size = 7, prob = 0.8, p0 = 0.40, lower = FALSE, log = TRUE)
+     Rzmgeomp0lt        == qzmgeom       (log1p(-Pzmgeomp0lt),        prob = pi/16, p0 = 0.01, lower = FALSE, log = TRUE)
+     Rzmgeomp0gt        == qzmgeom       (log1p(-Pzmgeomp0gt),        prob = pi/16, p0 = 0.40, lower = FALSE, log = TRUE)
+     Rzmbinomp0lt       == qzmbinom      (log1p(-Pzmbinomp0lt),       size = 12, prob = pi/16, p0 = 0.01, lower = FALSE, log = TRUE)
+     Rzmbinomp0gt       == qzmbinom      (log1p(-Pzmbinomp0gt),       size = 12, prob = pi/16, p0 = 0.12, lower = FALSE, log = TRUE)
+     Rzmlogarithmicp0lt == qzmlogarithmic(log1p(-Pzmlogarithmicp0lt), prob = 0.99, p0 = 0.05, lower = FALSE, log = TRUE)
+     Rzmlogarithmicp0gt == qzmlogarithmic(log1p(-Pzmlogarithmicp0gt), prob = 0.99, p0 = 0.55, lower = FALSE, log = TRUE)
+     Rpoisinvgauss    == qpoisinvgauss(log1p(-Ppoisinvgauss),    mean = 12,  dispersion = 0.1, lower = FALSE, log = TRUE)
+     RpoisinvgaussInf == qpoisinvgauss(log1p(-PpoisinvgaussInf), mean = Inf, dispersion = 1.1, lower = FALSE, log = TRUE)
+ })
> 
> proc.time()
   user  system elapsed 
  29.06    0.51   29.56