* using log directory 'd:/Rcompile/CRANpkg/local/3.6/actuar.Rcheck'
* using R version 3.6.3 (2020-02-29)
* using platform: x86_64-w64-mingw32 (64-bit)
* using session charset: ISO8859-1
* checking for file 'actuar/DESCRIPTION' ... OK
* checking extension type ... Package
* this is package 'actuar' version '3.1-2'
* package encoding: UTF-8
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking whether package 'actuar' can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking 'build' directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* loading checks for arch 'i386'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
** checking loading without being on the library search path ... OK
** checking use of S3 registration ... OK
* loading checks for arch 'x64'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
** checking loading without being on the library search path ... OK
** checking use of S3 registration ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... [16s] OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking contents of 'data' directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking pragmas in C/C++ headers and code ... OK
* checking compiled code ... OK
* checking sizes of PDF files under 'inst/doc' ... OK
* checking installed files from 'inst/doc' ... OK
* checking files in 'vignettes' ... OK
* checking examples ...
** running examples for arch 'i386' ... [5s] OK
** running examples for arch 'x64' ... [5s] OK
* checking for unstated dependencies in 'tests' ... OK
* checking tests ...
** running tests for arch 'i386' ... [36s] ERROR
  Running 'betaint-tests.R' [1s]
  Running 'dpqr-tests.R' [34s]
  Running 'rmixture-tests.R' [0s]
Running the tests in 'tests/dpqr-tests.R' failed.
Complete output:
  > ### 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 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(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),
  +               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(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),
  +               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(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),
  +               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(ppareto2(1e300 + m, min = m,
  +                        shape = 3,
  +                        scale = scLrg, log = TRUE),
  +               ppareto (1e300,
  +                        shape = 3,
  +                        scale = scLrg, log = TRUE),
  +               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))
   [1]  2.1501443 36.1801601  3.4857782  3.8152237 48.1620649  6.2800581
   [7]  0.8887332 20.3945333  5.7414426  3.7080842  1.1805086  6.6355089
  [13]  8.9991768  1.3915895  2.1947954  0.3214738  0.8047136 11.0704139
  [19]  3.9580371 21.4194474
  > (Rpareto4 <- rpareto4(n, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2))
   [1]   0.9630583  13.1761661   2.9329200   4.5090762   0.3199809   2.3155668
   [7]   0.3382328   0.3289525   0.4013289   2.3900373   0.5650040   1.1406143
  [13]   0.7557142 935.3464541   2.0994818   5.7709037   2.9784554   1.2264002
  [19]   3.3087978  11.3780756
  > (Rpareto3 <- rpareto3(n, min = m,               shape  = 1.5,             scale = 2))
   [1]  3.6959876  0.6942011  1.9357369  0.2723682  7.9135783  1.8210896
   [7] -0.4373613 -0.0745573 -0.0523196  5.0612333  6.5130680  0.7512528
  [13]  2.5195016  0.7371812  2.7414920  4.7509974  0.2920283  8.5263000
  [19]  1.6250735  1.3790832
  > (Rpareto2 <- rpareto2(n, min = m, shape  = 0.8,                           scale = 2))
   [1]  1.2270094  5.3510903  2.3353246 -0.4405288  2.4078013 -0.2479907
   [7] -0.3238822  1.1591109  3.5228019 19.3945399 -0.3860656  6.4016123
  [13] 63.7924661 -0.4216471  0.4519406 20.3242459  1.6689672 -0.4394662
  [19]  1.3446661  3.6391867
  > (Rtrbeta        <- rtrbeta      (n, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
   [1]  6.1237100  1.8744979  0.8919564  5.5252895  6.7227135  3.9531080
   [7] 13.1689817 93.6582516  1.9038216 14.9230278  1.4039983 15.3786102
  [13]  2.5191359  4.6901629  3.0622399  3.7123791  6.1707434  2.2434564
  [19]  8.8206621  2.7460235
  > (Rburr          <- rburr        (n, shape1 = 0.8, shape2 = 1.5,             scale = 2))
   [1] 1.6972393 7.5315819 0.6865106 1.1587441 1.5084275 1.7591859 3.6241865
   [8] 2.2737625 0.6427548 1.9574343 0.7808315 4.4158951 1.3245975 5.2549780
  [15] 1.8867210 2.6143859 5.2546735 2.0573305 0.4871449 0.5335860
  > (Rllogis        <- rllogis      (n, shape  = 1.5, scale = 2))
   [1] 3.8274680 3.2907796 0.1192988 1.4431102 0.3295236 2.1869917 2.5716056
   [8] 1.2881878 6.2799134 1.6445086 4.3330211 0.2304517 1.5205031 1.9214033
  [15] 2.6020944 0.5290367 2.9005850 3.6538487 5.7384436 5.6974725
  > (Rparalogis     <- rparalogis   (n, shape  = 0.8, scale = 2))
   [1]   3.29151302  13.37863602  17.92802167   1.13667610 226.01431848
   [6]   0.96763930   6.88912228   4.92951979   0.40678730   0.12897409
  [11]  10.80274500   0.04828211   0.81173250   1.05757854 191.61011753
  [16]   5.33184273   3.36726751   2.15896583   0.98331351   0.13597285
  > (Rgenpareto     <- rgenpareto   (n, shape1 = 0.8, shape2 = 2, scale = 2))
   [1]   9.1397476   6.1739488   1.3600727   6.3233986   1.8592948   1.3282090
   [7]   8.0056767   5.8415159   0.6771558   4.3369265   5.4286585   2.1379114
  [13]  54.4459206  22.6959404 418.9189235  20.9191590   2.9588252   1.7106009
  [19]  20.3113896 748.7666620
  > (Rpareto        <- rpareto      (n, shape  = 0.8, scale = 2))
   [1]  33.9962019   0.6482620   0.7085492 678.6061093   0.7325132   0.9670834
   [7]   1.5623838   2.9940088  18.2960986 806.6045582   3.3896018   2.8500929
  [13]   4.4979186   3.2131947   1.0510948  72.5769592   5.3042620   0.6318582
  [19]   0.5030023  10.0470496
  > (Rpareto1       <- rpareto1     (n, shape  = 0.8, min = 2))
   [1]   7.324825   2.425920   2.437014   9.164465  21.963669   3.101489
   [7]  33.942852 138.371106   2.001490 132.705913   7.749185   2.234682
  [13]   3.655667   9.550610   2.924884   2.509210   8.499629   4.846730
  [19]   3.138326   3.483903
  > (Rinvburr       <- rinvburr     (n, shape1 = 1.5, shape2 = 2, scale = 2))
   [1]  3.424701 16.306924  2.238696  1.131552  2.735385  1.532752  2.546524
   [8]  2.063354 18.850984  2.133496  1.547773  3.285891  1.201398 13.395263
  [15]  2.681530  1.304574  3.276841 20.546751  3.605976  2.255694
  > (Rinvpareto     <- rinvpareto   (n, shape  = 2,   scale = 2))
   [1]   2.6364375  21.2350209   1.2218165   1.5657496  35.7094332   2.4952425
   [7]   3.0345392  15.4531524   1.5698249   0.3075728   3.5196489   4.5595843
  [13]   3.7060677   2.8250761  26.9229897   4.1466020   5.4236884 107.6204565
  [19]  14.6817948   1.6834442
  > (Rinvparalogis  <- rinvparalogis(n, shape  = 2,   scale = 2))
   [1]  2.236658 16.527069  3.606370  5.226055  2.503909  5.341359  3.316043
   [8]  9.325952  1.738741  2.129130  1.334796  1.841420 18.528845  2.185957
  [15]  4.799531  5.586437  1.386768  1.958369  2.082247  1.365363
  > (Rtrgamma       <- rtrgamma     (n, shape1 = 2, shape2 = 3, scale = 5))
   [1] 3.682281 3.858648 5.302705 3.276745 4.955305 5.715560 8.069265 5.639096
   [9] 6.286054 4.978647 4.888418 7.318054 5.970505 9.402681 6.161277 5.672815
  [17] 4.790313 4.085980 6.778603 5.142050
  > (Rinvtrgamma    <- rinvtrgamma  (n, shape1 = 2, shape2 = 3, scale = 5))
   [1] 4.896678 5.502744 3.720271 6.670609 4.478685 5.287671 5.219636 3.529218
   [9] 3.699375 4.542030 5.643312 3.550759 3.868373 6.828480 4.594627 8.059268
  [17] 7.495471 9.208972 4.880526 2.914834
  > (Rinvgamma      <- rinvgamma    (n, shape  = 2,             scale = 5))
   [1]  1.2205324  3.5881100  4.1895661  3.6655585  1.7861110  4.5201158
   [7]  3.2086082  1.4224586  0.8892067  1.2518328  3.7275703  2.0948572
  [13]  1.7293930  5.8413666  1.7179378  2.0365153 17.9958054  4.2963595
  [19]  2.9040597  3.1393318
  > (Rinvweibull    <- rinvweibull  (n,             shape  = 3, scale = 5))
   [1]  7.224733 12.787122 18.517641  3.293167  3.358632  9.643576  7.464673
   [8]  7.407516  4.612062  5.849253 10.994832  6.881836  8.940351  4.725462
  [15]  4.528093  3.250087  9.783872  4.633481  4.752201 10.601778
  > (Rinvexp        <- rinvexp      (n,                         scale = 5))
   [1]  14.996537 118.677380  81.634288 298.663312   4.793279 225.621141
   [7]  13.405519   3.318091   8.293811   3.212230   7.746716   4.915751
  [13]   4.202300  19.268502   1.904340   1.290878  12.353135   9.809097
  [19] 163.396473   9.077455
  > (Rlgamma <- rlgamma(n, shapelog = 1.5, ratelog = 5))
   [1] 1.152755 1.102730 1.117140 1.196528 1.078281 1.527771 1.284417 1.678350
   [9] 1.022102 1.020058 1.433508 1.206858 1.601457 4.158848 1.507774 1.197388
  [17] 1.390355 1.673158 1.253893 1.056516
  > (Rgumbel <- rgumbel(n, alpha = 2, scale = 5))
   [1] -3.56843891  8.55406059 14.68383243 -0.05101188  2.94374512 10.97699471
   [7] -0.67222684 -2.81636644  2.31074792  3.54843150 -4.86475796  9.30684196
  [13] 13.88371439  1.09287275 -1.13845047  1.73624280  7.70901865  5.89655608
  [19]  3.81417732 14.96592464
  > (Rinvgauss <- rinvgauss(n, mean = 2, dispersion = 5))
   [1] 0.13481571 0.16104316 0.04773399 0.64341790 8.10365104 7.73392520
   [7] 0.15611086 0.51332331 0.08342178 0.37838020 0.08553184 1.91153584
  [13] 0.22865342 1.10550916 0.15152142 8.03631756 0.20389745 0.13487839
  [19] 4.13259400 0.33219371
  > (Rgenbeta <- rgenbeta(n, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
   [1] 0.8800450 1.2536617 1.0975603 1.1998557 1.3560796 1.5922056 1.6304331
   [8] 1.4349133 0.1705367 1.2221339 0.8540434 0.7998316 0.4017762 1.8316415
  [15] 1.0605475 1.5362786 0.7223314 0.6452797 0.7486176 1.0404133
  > 
  > ## Compute quantiles
  > (Pfpareto <- pfpareto(Rfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
   [1] 0.30163096 0.94610852 0.46104755 0.49235515 0.96138464 0.65626515
   [7] 0.11131189 0.89641355 0.62873489 0.48248468 0.15566788 0.67258567
  [13] 0.75394437 0.18846984 0.30785062 0.03818643 0.09905751 0.80010545
  [19] 0.50506625 0.90194222
  > (Ppareto4 <- ppareto4(Rpareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2))
   [1] 0.3348848 0.9051593 0.6160304 0.7256164 0.1853600 0.5514837 0.1899356
   [8] 0.1876105 0.2056577 0.5601683 0.2455248 0.3707789 0.2898176 0.9993752
  [15] 0.5246834 0.7798811 0.6201835 0.3872290 0.6482021 0.8888646
  > (Ppareto3 <- ppareto3(Rpareto3, min = m,               shape  = 1.5,             scale = 2))
   [1] 0.75638053 0.33194441 0.58235322 0.21180417 0.89713536 0.56510726
   [7] 0.01504302 0.10694889 0.11353094 0.82494735 0.86930431 0.34689807
  [13] 0.65648353 0.34324187 0.67962676 0.81230888 0.21770570 0.90640501
  [19] 0.53322096 0.48849454
  > (Ppareto2 <- ppareto2(Rpareto2, min = m, shape  = 0.8,                           scale = 2))
   [1] 0.40001104 0.66717646 0.51138697 0.04552617 0.51709760 0.10964978
   [7] 0.08556181 0.39126502 0.58931022 0.85290531 0.06470020 0.69877110
  [13] 0.93927943 0.05227306 0.27940373 0.85770570 0.45071534 0.04590876
  [19] 0.41451665 0.59548972
  > (Ptrbeta        <- ptrbeta      (Rtrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
   [1] 0.61879377 0.17640768 0.03921634 0.58191841 0.65054577 0.45166676
   [7] 0.82535895 0.98225934 0.18096926 0.84790455 0.10467761 0.85292133
  [13] 0.27430685 0.51969709 0.34895912 0.42628087 0.62146063 0.23331197
  [19] 0.73266434 0.30656737
  > (Pburr          <- pburr        (Rburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2))
   [1] 0.37002695 0.81617151 0.13635589 0.25343200 0.33171535 0.38198213
   [7] 0.62776194 0.47016431 0.12531766 0.41824990 0.16023224 0.68755184
  [13] 0.29170968 0.73498219 0.40565681 0.51871020 0.73496727 0.43540957
  [19] 0.08681178 0.09812561
  > (Pllogis        <- pllogis      (Rllogis,       shape  = 1.5, scale = 2))
   [1] 0.72583338 0.67851724 0.01435912 0.38000669 0.06268591 0.53346730
   [7] 0.59316741 0.34076968 0.84765338 0.42713294 0.76127397 0.03764106
  [13] 0.39863427 0.48497027 0.59742666 0.11975346 0.63590814 0.71176072
  [19] 0.82935476 0.82782825
  > (Pparalogis     <- pparalogis   (Rparalogis,    shape  = 0.8, scale = 2))
   [1] 0.51795739 0.74703507 0.78374521 0.32562361 0.95233637 0.29914691
   [7] 0.64811136 0.59105600 0.17904868 0.08114265 0.71747167 0.03889521
  [13] 0.27160557 0.31362492 0.94715606 0.60477987 0.52215462 0.43974674
  [19] 0.30173841 0.08431646
  > (Pgenpareto     <- pgenpareto   (Rgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2))
   [1] 0.58074223 0.47982175 0.12587420 0.48618404 0.18116201 0.12221888
   [7] 0.54762905 0.46503498 0.04782608 0.38488687 0.44534844 0.20997091
  [13] 0.87756744 0.76768562 0.97512176 0.75409990 0.28549993 0.16514790
  [19] 0.74903177 0.98432607
  > (Ppareto        <- ppareto      (Rpareto,       shape  = 0.8, scale = 2))
   [1] 0.9009582 0.2011684 0.2154247 0.9905701 0.2209341 0.2706093 0.3698681
   [8] 0.5190892 0.8433631 0.9917845 0.5475416 0.5077067 0.6104129 0.5353341
  [15] 0.2867210 0.9446981 0.6452169 0.1971877 0.1642912 0.7622506
  > (Ppareto1       <- ppareto1     (Rpareto1,      shape  = 0.8, min = 2))
   [1] 0.6460138974 0.1431145798 0.1462366262 0.7041045446 0.8529516754
   [6] 0.2960079387 0.8961933118 0.9662722270 0.0005954721 0.9651251959
  [11] 0.6616087162 0.0849362384 0.3827647304 0.7137146513 0.2622025970
  [16] 0.1659456908 0.6857292196 0.5074334522 0.3026262338 0.3585376460
  > (Pinvburr       <- pinvburr     (Rinvburr,      shape1 = 1.5, shape2 = 2, scale = 2))
   [1] 0.6439229 0.9778534 0.4147353 0.1194048 0.5260297 0.2250734 0.4864118
   [8] 0.3702148 0.9833502 0.3883191 0.2292448 0.6232975 0.1365402 0.9674695
  [15] 0.5150718 0.1630703 0.6219023 0.9859542 0.6687715 0.4189159
  > (Pinvpareto     <- pinvpareto   (Rinvpareto,    shape  = 2,   scale = 2))
   [1] 0.32334499 0.83525532 0.14381704 0.19281595 0.89673868 0.30811955
   [7] 0.36330054 0.78394648 0.19337868 0.01776581 0.40660787 0.48316767
  [13] 0.42184495 0.34280880 0.86648331 0.45510805 0.53376487 0.96384333
  [19] 0.77459154 0.20887635
  > (Pinvparalogis  <- pinvparalogis(Rinvparalogis, shape  = 2,   scale = 2))
   [1] 0.30878683 0.97134245 0.58490009 0.76082363 0.37270939 0.76919391
   [7] 0.53767718 0.91399545 0.18529644 0.28221842 0.09496241 0.21048708
  [13] 0.97709900 0.29630218 0.72598303 0.78568783 0.10541775 0.23959463
  [19] 0.27054487 0.10105849
  > (Ptrgamma       <- ptrgamma     (Rtrgamma,      shape1 = 2, shape2 = 3, scale = 5))
   [1] 0.06139920 0.07821445 0.33478350 0.03290653 0.25446472 0.44006711
   [7] 0.92223148 0.42003842 0.59049704 0.25954807 0.24017528 0.82016714
  [13] 0.50757438 0.99010382 0.55798855 0.42885047 0.21998559 0.10437553
  [19] 0.71101263 0.29645511
  > (Pinvtrgamma    <- pinvtrgamma  (Rinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5))
   [1] 0.71199284 0.82657373 0.30246923 0.93269513 0.59479812 0.79234963
   [7] 0.78015650 0.22375134 0.29371407 0.61482121 0.84575184 0.23240962
  [13] 0.36458699 0.94042412 0.63082385 0.97564582 0.96376476 0.98847893
  [19] 0.70810050 0.03886046
  > (Pinvgamma      <- pinvgamma    (Rinvgamma,     shape  = 2,             scale = 5))
   [1] 0.08475380 0.59408184 0.66499919 0.60430698 0.23118414 0.69677206
   [7] 0.53850316 0.13430706 0.02393430 0.09200818 0.61224305 0.31132662
  [13] 0.21600077 0.78854974 0.21292621 0.29662060 0.96785848 0.67575740
  [19] 0.48652799 0.52729373
  > (Pinvweibull    <- pinvweibull  (Rinvweibull,               shape  = 3, scale = 5))
   [1] 0.71786749 0.94196711 0.98050670 0.03019718 0.03690869 0.86989874
   [7] 0.74043087 0.73525834 0.27966517 0.53547040 0.91024013 0.68145298
  [13] 0.83952188 0.30586397 0.26018332 0.02622513 0.87505533 0.28462757
  [19] 0.31200659 0.90041490
  > (Pinvexp        <- pinvexp      (Rinvexp,                               scale = 5))
   [1] 0.71647615 0.95874415 0.94058922 0.98339810 0.35235110 0.97808271
   [7] 0.68867852 0.22159791 0.54724473 0.21086212 0.52443448 0.36162826
  [13] 0.30427459 0.77144427 0.07239765 0.02078895 0.66713986 0.60065719
  [19] 0.96986304 0.57647968
  > (Plgamma <- plgamma(Rlgamma, shapelog = 1.5, ratelog = 5))
   [1] 0.29950795 0.19339871 0.22478912 0.38381094 0.13949021 0.76312504
   [7] 0.52525935 0.84079002 0.02546943 0.02218361 0.69213357 0.40236132
  [13] 0.80562288 0.99741901 0.74979268 0.38537480 0.65174313 0.83866479
  [19] 0.48026319 0.09217148
  > (Pgumbel <- pgumbel(Rgumbel, alpha = 2, scale = 5))
   [1] 0.04756923 0.76368413 0.92392736 0.22154650 0.43692475 0.84699440
   [7] 0.18149966 0.07278399 0.39072851 0.48014023 0.01931287 0.79301059
  [13] 0.91132790 0.30151812 0.15361920 0.34848239 0.72669918 0.63209103
  [19] 0.49872477 0.92794623
  > (Pinvgauss <- pinvgauss(Rinvgauss, mean = 2, dispersion = 5))
   [1] 0.24621364 0.29232542 0.04490232 0.63430088 0.94520442 0.94255315
   [7] 0.28416016 0.58564314 0.13412747 0.51427731 0.13930376 0.81612195
  [13] 0.38532532 0.73539198 0.27635988 0.94473692 0.35489802 0.24633213
  [19] 0.89805579 0.48205524
  > (Pgenbeta <- pgenbeta(Rgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2))
   [1] 0.3631563 0.6066533 0.5029626 0.5708196 0.6745405 0.8241147 0.8466035
   [8] 0.7259740 0.0274693 0.5856566 0.3471103 0.3142575 0.1074184 0.9501596
  [15] 0.4786289 0.7901209 0.2688807 0.2258892 0.2840466 0.4654637
  > 
  > ## Just compute pdf
  > dfpareto(Rfpareto, min = m, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
   [1] 0.1399564518 0.0017259176 0.0994225091 0.0907703847 0.0009388851
   [6] 0.0478642379 0.1476999952 0.0056691695 0.0545545221 0.0934980864
  [11] 0.1548802412 0.0440370395 0.0267662682 0.1553789105 0.1386298469
  [16] 0.1023132847 0.1438539291 0.0185175481 0.0872640267 0.0051322165
  > dpareto4(Rpareto4, min = m, shape1 = 0.8, shape2 = 1.5,             scale = 2)
   [1] 2.092090e-01 7.849000e-03 9.202973e-02 5.205075e-02 2.509939e-01
   [6] 1.184500e-01 2.503886e-01 2.507041e-01 2.478546e-01 1.148100e-01
  [11] 2.387845e-01 1.951225e-01 2.253605e-01 8.009920e-07 1.298032e-01
  [16] 3.542945e-02 9.038857e-02 1.884102e-01 7.954462e-02 1.045395e-02
  > dpareto3(Rpareto3, min = m,               shape  = 1.5,             scale = 2)
   [1] 0.06493737 0.26511688 0.14615219 0.30067431 0.01633519 0.15478962
   [7] 0.18052404 0.29483602 0.29707902 0.03853167 0.02409280 0.25907778
  [13] 0.10982853 0.26057762 0.09891113 0.04305675 0.29966416 0.01400414
  [19] 0.17082414 0.19324056
  > dpareto2(Rpareto2, min = m, shape  = 0.8,                           scale = 2)
   [1] 0.126730843 0.033654378 0.079841989 0.360187839 0.077757727 0.308015060
   [7] 0.327082484 0.130925300 0.054009014 0.005359857 0.344111660 0.026889172
  [13] 0.000732092 0.354484500 0.191366986 0.004974299 0.103897308 0.359863073
  [19] 0.119941018 0.052197717
  > dtrbeta      (Rtrbeta,       shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
   [1] 0.0570369577 0.1555673503 0.1035237196 0.0665148627 0.0492306884
   [6] 0.1020089739 0.0145463164 0.0002262037 0.1555426982 0.0113490458
  [11] 0.1457415128 0.0106841021 0.1453431293 0.0832438442 0.1291222322
  [16] 0.1089618705 0.0563677377 0.1517203077 0.0309678768 0.1389063061
  > dburr        (Rburr,         shape1 = 0.8, shape2 = 1.5,             scale = 2)
   [1] 0.19542595 0.02576370 0.25276283 0.23661164 0.21040779 0.19056558
   [7] 0.08741546 0.15322390 0.25166466 0.17544334 0.25308891 0.06507229
  [13] 0.22472617 0.04901066 0.18074821 0.13235370 0.04901545 0.16814684
  [19] 0.24139375 0.24565018
  > dllogis      (Rllogis,       shape  = 1.5, scale = 2)
   [1] 0.07798862 0.09942853 0.17795151 0.24488942 0.26746060 0.17070019
   [7] 0.14076022 0.26158341 0.03084528 0.22318861 0.06291312 0.23578185
  [13] 0.23649244 0.19499350 0.13864296 0.29888067 0.11973221 0.08422245
  [19] 0.03699403 0.03752418
  > dparalogis   (Rparalogis,    shape  = 0.8, scale = 2)
   [1] 0.0560814287 0.0099302302 0.0065814624 0.1476586632 0.0001319623
   [6] 0.1662930150 0.0238306012 0.0357304669 0.2822907628 0.4576818499
  [11] 0.0132904186 0.6163577975 0.1878444095 0.1558679582 0.0001720329
  [16] 0.0325737981 0.0547390283 0.0855798878 0.1643835775 0.4493582320
  > dgenpareto   (Rgenpareto,    shape1 = 0.8, shape2 = 2, scale = 2)
   [1] 2.684539e-02 4.314596e-02 1.145440e-01 4.200443e-02 1.062460e-01
   [6] 1.148850e-01 3.176102e-02 4.585556e-02 1.077430e-01 6.181705e-02
  [11] 4.958267e-02 1.005053e-01 1.700516e-03 7.174359e-03 4.715859e-05
  [16] 8.150315e-03 8.380022e-02 1.091167e-01 8.532015e-03 1.667704e-05
  > dpareto      (Rpareto,       shape  = 0.8, scale = 2)
   [1] 2.201162e-03 2.413150e-01 2.317330e-01 1.108413e-05 2.280877e-01
   [6] 1.966620e-01 1.415079e-01 7.703804e-02 6.174067e-03 8.128096e-06
  [11] 6.716020e-02 8.120147e-02 4.796454e-02 7.130613e-02 1.870224e-01
  [16] 5.932334e-04 3.885765e-02 2.440290e-01 2.671061e-01 1.578806e-02
  > dpareto1     (Rpareto1,      shape  = 0.8, min = 2)
   [1] 0.0386615239 0.2825765942 0.2802653518 0.0258298066 0.0053560568
   [6] 0.1815881296 0.0024466227 0.0001949989 0.3994642746 0.0002102381
  [11] 0.0349343864 0.3275861591 0.1350747410 0.0239804880 0.2017987707
  [16] 0.2659176934 0.0295797160 0.0813029020 0.1777696449 0.1472974173
  > dinvburr     (Rinvburr,      shape1 = 1.5, shape2 = 2, scale = 2)
   [1] 0.143450771 0.002665966 0.246687663 0.239806555 0.200975476 0.277527173
   [7] 0.218614021 0.260745348 0.001741912 0.255399728 0.277902414 0.153832313
  [13] 0.250546302 0.004724867 0.205974181 0.263067432 0.154532055 0.001351179
  [19] 0.130890753 0.245217978
  > dinvpareto   (Rinvpareto,    shape  = 2,   scale = 2)
   [1] 0.1058093978 0.0067714782 0.1461381396 0.1381432821 0.0026637468
   [6] 0.1098786664 0.0951203778 0.0116266724 0.1380290387 0.1001248075
  [11] 0.0837191821 0.0646184123 0.0797925937 0.1005952892 0.0044509604
  [16] 0.0714244768 0.0530268126 0.0003267985 0.0126505910 0.1347399648
  > dinvparalogis(Rinvparalogis, shape  = 2,   scale = 2)
   [1] 0.245363072 0.003393056 0.152591923 0.074391606 0.231910341 0.070830207
   [7] 0.172998301 0.017236796 0.242781625 0.248537305 0.196880438 0.247456776
  [13] 0.002429309 0.247057222 0.089518595 0.063913177 0.205342456 0.249834140
  [19] 0.249391861 0.201945679
  > dtrgamma     (Rtrgamma,      shape1 = 2, shape2 = 3, scale = 5)
   [1] 0.08717975 0.10372132 0.24419934 0.05473657 0.21672116 0.26295719
   [7] 0.09817362 0.26081199 0.25834262 0.21882842 0.21051358 0.17524394
  [13] 0.26540482 0.01825347 0.26244407 0.26183651 0.20100476 0.12669934
  [19] 0.22742261 0.23260109
  > dinvtrgamma  (Rinvtrgamma,   shape1 = 2, shape2 = 3, scale = 5)
   [1] 0.239475145 0.144904524 0.419376163 0.052346715 0.322553232 0.174137245
   [7] 0.184380287 0.400142684 0.418562405 0.309625583 0.128275432 0.403697507
  [13] 0.417296637 0.045727103 0.298885435 0.016717260 0.026208108 0.007111335
  [19] 0.242492896 0.168493379
  > dinvgamma    (Rinvtrgamma,   shape  = 2,             scale = 5)
   [1] 0.07669678 0.06047637 0.12662797 0.03980294 0.09112605 0.06568703
   [7] 0.06745272 0.13791778 0.12781189 0.08873723 0.05735231 0.13659188
  [13] 0.11858096 0.03775422 0.08681166 0.02568120 0.03046763 0.01859974
  [19] 0.07719945 0.18160405
  > dinvweibull  (Rinvweibull,               shape  = 3, scale = 5)
   [1] 0.098807142 0.013212221 0.003127082 0.096281495 0.108770321 0.037717901
   [7] 0.089427837 0.091575829 0.231787100 0.171539833 0.023357847 0.113933111
  [13] 0.049277078 0.230028926 0.232086259 0.088139184 0.035811597 0.231568162
  [19] 0.229411851 0.026727552
  > dinvexp      (Rinvexp,                               scale = 5)
   [1] 1.592905e-02 3.403587e-04 7.057076e-04 5.512335e-05 7.667965e-02
   [6] 9.606960e-05 1.916105e-02 1.006374e-01 3.977803e-02 1.021775e-01
  [11] 4.369442e-02 7.482600e-02 8.615127e-02 1.038912e-02 9.981722e-02
  [16] 6.237813e-02 2.185911e-02 3.121322e-02 1.816331e-04 3.498048e-02
  > dlgamma(Rlgamma, shapelog = 1.5, ratelog = 5)
   [1] 2.027070049 2.194015355 2.160125982 1.821013385 2.203495968 0.645870840
   [7] 1.405752448 0.406163236 1.635988340 1.578126848 0.872440662 1.770414531
  [13] 0.513202765 0.002910821 0.688049737 1.816807835 1.002588406 0.412544321
  [19] 1.543989633 2.126882358
  > dgumbel(Rgumbel, alpha = 2, scale = 5)
   [1] 0.02897508 0.04117800 0.01462056 0.06677956 0.07235424 0.02813058
   [7] 0.06194589 0.03814258 0.07343682 0.07045358 0.01524552 0.03678280
  [13] 0.01692382 0.07229953 0.05755431 0.07347177 0.04639868 0.05799079
  [19] 0.06939265 0.01387864
  > dinvgauss(Rinvgauss, mean = 2, dispersion = 5)
   [1] 1.890778452 1.633114908 2.324193919 0.321833272 0.006894296 0.007458802
   [7] 1.678083963 0.435602357 2.462797824 0.644282631 2.443298153 0.067500493
  [13] 1.157896367 0.150738266 1.721376411 0.006992163 1.304747571 1.890109068
  [19] 0.020660571 0.755832129
  > dgenbeta(Rgenbeta, shape1 = 0.8, shape2 = 1.5, shape3 = 2, scale = 2)
   [1] 0.6204901 0.6657241 0.6594897 0.6658230 0.6583174 0.5968054 0.5793481
   [8] 0.6454206 0.2571983 0.6660817 0.6136692 0.5980234 0.4228707 0.4307708
  [15] 0.6552199 0.6180389 0.5723392 0.5429270 0.5814873 0.6524751
  > 
  > ## 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), tol = 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]), tol = 1e-8)
  +     })
  > for (i in seq_along(r))
  +     stopifnot(exprs = {
  +         all.equal(cumsum(dztgeom(0:20, prob[i])),
  +                   pztgeom(0:20, prob[i]), tol = 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]), tol = 1e-8)
  +     })
  > for (p in prob)
  +     stopifnot(exprs = {
  +         all.equal(cumsum(dlogarithmic(0:20, p)),
  +                   plogarithmic(0:20, p), tol = 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]), tol = 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]), tol = 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]), tol = 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]), tol = 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]), tol = 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]), tol = 1e-8)
  +         all.equal(cumsum(dpoisinvgauss(0:20, Inf, phi[i])),
  +                   ppoisinvgauss(0:20, Inf, phi[i]), tol = 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))
   [1] 10 15 11 16 18  6 12 16 12 11 18 11 13 12  8 17 10  6 10 18
  > (Rztnbinom    <- rztnbinom   (n, size = 7, prob = 0.01))
   [1] 1025  798  757 1534  768  811  689  723  528  427 1228 1049  796  894  276
  [16]  646  856  479  546  490
  > (Rztgeom      <- rztgeom     (n, prob = pi/16))
   [1]  7  6  5 10 13 12  1  5  3  5  3  3  1  5  3  5  9  5  2  2
  > (Rztbinom     <- rztbinom    (n, size = 55, prob = pi/16))
   [1]  8 13  7 10 11 11  9 11 16 11 14 15 12 10  8 15  9  6 16 12
  > (Rlogarithmic <- rlogarithmic(n, prob = 0.99))
   [1] 24  1 18  5  8  1  4  7  3  4 19 22  4  1 82  4  7 16  4  2
  > (Rzmpoisp0lt        <- rzmpois       (n, lambda = 6, p0 = 0.001))
   [1] 5 6 4 5 4 6 5 7 5 5 6 7 4 5 4 7 4 9 8 7
  > (Rzmpoisp0gt        <- rzmpois       (n, lambda = 6, p0 = 0.010))
   [1]  5  9  8  7  6  4  9  4 12 10  5  4  4  6  5  5  4  6  4  2
  > (Rzmnbinomp0lt      <- rzmnbinom     (n, size = 7, prob = 0.8, p0 = 0.01))
   [1] 1 1 4 1 2 1 4 2 1 2 3 2 1 3 2 1 2 1 2 1
  > (Rzmnbinomp0gt      <- rzmnbinom     (n, size = 7, prob = 0.8, p0 = 0.40))
   [1] 1 0 0 1 0 0 1 2 5 1 2 0 3 1 0 0 3 3 2 4
  > (Rzmgeomp0lt        <- rzmgeom       (n, prob = pi/16, p0 = 0.01))
   [1]  1  2  1  2 18  2  6  8  1  2  2  1  1 22 20  1 11  4  3  3
  > (Rzmgeomp0gt        <- rzmgeom       (n, prob = pi/16, p0 = 0.40))
   [1]  8  6  3  0  7  3  1  0  0  7 10  1  1  9  0  0  9  9  2  0
  > (Rzmbinomp0lt       <- rzmbinom      (n, size = 12, prob = pi/16, p0 = 0.01))
   [1] 3 2 3 3 1 4 2 3 5 1 1 3 2 4 5 3 1 6 1 3
  > (Rzmbinomp0gt       <- rzmbinom      (n, size = 12, prob = pi/16, p0 = 0.12))
   [1] 2 1 1 2 2 1 2 1 3 4 2 0 3 0 0 5 2 0 2 2
  > (Rzmlogarithmicp0lt <- rzmlogarithmic(n, prob = 0.99, p0 = 0.05))
   [1]   1  41 127   4  34  12   1   0   1   3  36   2   2  53   6  56  27   6   1
  [20]   2
  > (Rzmlogarithmicp0gt <- rzmlogarithmic(n, prob = 0.99, p0 = 0.55))
   [1]  0  9  6  0 84  0  0  1 20 26  0  6  0  0  0  0  4  2  6 13
  > (Rpoisinvgauss    <- rpoisinvgauss(n, mean = 12,  dispersion = 0.1))
   [1] 19  1 27 19 15 11  1 10  4 29 10  1  3 12  7  5  7 97  8 23
  > (RpoisinvgaussInf <- rpoisinvgauss(n, mean = Inf, dispersion = 1.1))
   [1]  15   0  13   1  14   0   3   0 962   1  28   0 169   2 128  12   1  61   0
  [20]   8
  > 
  > ## Compute quantiles
  > (Pztpois      <- pztpois     (Rztpois,      lambda = 12))
   [1] 0.34722541 0.84441470 0.46159402 0.89870837 0.96258328 0.04581644
   [7] 0.57596264 0.89870837 0.57596264 0.46159402 0.96258328 0.46159402
  [13] 0.68153368 0.57596264 0.15502259 0.93703332 0.34722541 0.04581644
  [19] 0.34722541 0.96258328
  > (Pztnbinom    <- pztnbinom   (Rztnbinom,    size = 7, prob = 0.01))
   [1] 0.88997750 0.69393730 0.64180911 0.99430218 0.65628911 0.70939020
   [7] 0.54503524 0.59489596 0.29006618 0.14800061 0.96302942 0.90254438
  [13] 0.69151318 0.79550085 0.02498486 0.47861352 0.75876748 0.21682950
  [19] 0.31830940 0.23268782
  > (Pztgeom      <- pztgeom     (Rztgeom,      prob = pi/16))
   [1] 0.7834938 0.7305965 0.6647753 0.8876244 0.9416725 0.9274218 0.1963495
   [8] 0.6647753 0.4809591 0.6647753 0.4809591 0.4809591 0.1963495 0.6647753
  [15] 0.4809591 0.6647753 0.8601686 0.6647753 0.3541459 0.3541459
  > (Pztbinom     <- pztbinom    (Rztbinom,     size = 55, prob = pi/16))
   [1] 0.2216795 0.8220816 0.1288058 0.4733550 0.6064659 0.6064659 0.3401774
   [8] 0.6064659 0.9686998 0.6064659 0.8927166 0.9398875 0.7257130 0.4733550
  [15] 0.2216795 0.9398875 0.3401774 0.0654511 0.9686998 0.7257130
  > (Plogarithmic <- plogarithmic(Rlogarithmic, prob = 0.99))
   [1] 0.7706793 0.2149758 0.7214703 0.4850700 0.5731034 0.2149758 0.4437691
   [8] 0.5480570 0.3916214 0.4437691 0.7309124 0.7560780 0.4437691 0.2149758
  [15] 0.9359981 0.4437691 0.5480570 0.7006357 0.4437691 0.3213888
  > (Pzmpoisp0lt        <- pzmpois       (Rzmpoisp0lt,        lambda = 6, p0 = 0.001))
   [1] 0.4448579 0.6057192 0.2839966 0.4448579 0.2839966 0.6057192 0.4448579
   [8] 0.7436002 0.4448579 0.4448579 0.6057192 0.7436002 0.2839966 0.4448579
  [15] 0.2839966 0.7436002 0.2839966 0.9159516 0.8470110 0.7436002
  > (Pzmpoisp0gt        <- pzmpois       (Rzmpoisp0gt,        lambda = 6, p0 = 0.010))
   [1] 0.4498592 0.9167088 0.8483893 0.7459101 0.6092712 0.2904471 0.9167088
   [8] 0.2904471 0.9912391 0.9577004 0.4498592 0.2904471 0.2904471 0.6092712
  [15] 0.4498592 0.4498592 0.2904471 0.6092712 0.2904471 0.0690415
  > (Pzmnbinomp0lt      <- pzmnbinom     (Rzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01))
   [1] 0.3777981 0.3777981 0.9368513 0.3777981 0.6720366 0.3777981 0.9368513
   [8] 0.6720366 0.3777981 0.6720366 0.8485797 0.6720366 0.3777981 0.8485797
  [15] 0.6720366 0.3777981 0.6720366 0.3777981 0.6720366 0.3777981
  > (Pzmnbinomp0gt      <- pzmnbinom     (Rzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40))
   [1] 0.6229080 0.4000000 0.4000000 0.6229080 0.4000000 0.4000000 0.6229080
   [8] 0.8012343 0.9852671 0.6229080 0.8012343 0.4000000 0.9082301 0.6229080
  [15] 0.4000000 0.4000000 0.9082301 0.9082301 0.8012343 0.9617280
  > (Pzmgeomp0lt        <- pzmgeom       (Rzmgeomp0lt,        prob = pi/16, p0 = 0.01))
   [1] 0.2043860 0.3606045 0.2043860 0.3606045 0.9806427 0.3606045 0.7332906
   [8] 0.8277446 0.2043860 0.3606045 0.3606045 0.2043860 0.2043860 0.9919255
  [15] 0.9874980 0.2043860 0.9105924 0.5870438 0.4861495 0.4861495
  > (Pzmgeomp0gt        <- pzmgeom       (Rzmgeomp0gt,        prob = pi/16, p0 = 0.40))
   [1] 0.8956028 0.8383579 0.6885755 0.4000000 0.8700963 0.6885755 0.5178097
   [8] 0.4000000 0.4000000 0.8700963 0.9325746 0.5178097 0.5178097 0.9161011
  [15] 0.4000000 0.4000000 0.9161011 0.9161011 0.6124876 0.4000000
  > (Pzmbinomp0lt       <- pzmbinom      (Rzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01))
   [1] 0.7909672 0.5423821 0.7909672 0.7909672 0.2371476 0.9276205 0.5423821
   [8] 0.7909672 0.9810404 0.2371476 0.2371476 0.7909672 0.5423821 0.9276205
  [15] 0.9810404 0.7909672 0.2371476 0.9962673 0.2371476 0.7909672
  > (Pzmbinomp0gt       <- pzmbinom      (Rzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12))
   [1] 0.5932285 0.3219090 0.3219090 0.5932285 0.5932285 0.3219090 0.5932285
   [8] 0.3219090 0.8141930 0.9356627 0.5932285 0.1200000 0.8141930 0.1200000
  [15] 0.1200000 0.9831470 0.5932285 0.1200000 0.5932285 0.5932285
  > (Pzmlogarithmicp0lt <- pzmlogarithmic(Rzmlogarithmicp0lt, prob = 0.99, p0 = 0.05))
   [1] 0.2542270 0.8608469 0.9712412 0.4715806 0.8348042 0.6660721 0.2542270
   [8] 0.0500000 0.2542270 0.4220403 0.8429410 0.3553193 0.3553193 0.8934536
  [15] 0.5431862 0.8999300 0.8004979 0.5431862 0.2542270 0.3553193
  > (Pzmlogarithmicp0gt <- pzmlogarithmic(Rzmlogarithmicp0gt, prob = 0.99, p0 = 0.55))
   [1] 0.5500000 0.8178149 0.7836145 0.5500000 0.9722104 0.5500000 0.5500000
   [8] 0.6467391 0.8829067 0.9027400 0.5500000 0.7836145 0.5500000 0.5500000
  [15] 0.5500000 0.5500000 0.7496961 0.6946249 0.7836145 0.8484196
  > (Ppoisinvgauss    <- ppoisinvgauss(Rpoisinvgauss,    mean = 12,  dispersion = 0.1))
   [1] 0.82602138 0.07783346 0.90272595 0.82602138 0.75958956 0.65734419
   [7] 0.07783346 0.62339905 0.30314405 0.91502328 0.62339905 0.07783346
  [13] 0.22730218 0.68745450 0.49217525 0.37323559 0.49217525 0.99775681
  [19] 0.54152275 0.87109211
  > (PpoisinvgaussInf <- ppoisinvgauss(RpoisinvgaussInf, mean = Inf, dispersion = 1.1))
   [1] 0.8072193 0.2596554 0.7935270 0.4347151 0.8007243 0.2596554 0.6021435
   [8] 0.2596554 0.9754794 0.4347151 0.8576561 0.2596554 0.9415763 0.5374926
  [15] 0.9329037 0.7854913 0.4347151 0.9030378 0.2596554 0.7404399
  > 
  > ## Just compute pmf
  > dztpois     (Rztpois,      lambda = 12)
   [1] 0.10483790 0.07239156 0.11436862 0.05429367 0.02554996 0.02548143
   [7] 0.11436862 0.05429367 0.11436862 0.11436862 0.02554996 0.11436862
  [13] 0.10557103 0.11436862 0.06552369 0.03832495 0.10483790 0.02548143
  [19] 0.10483790 0.02554996
  > dztnbinom   (Rztnbinom,    size = 7, prob = 0.01)
   [1] 5.520319e-04 1.210509e-03 1.333824e-03 3.697815e-05 1.301686e-03
   [6] 1.169906e-03 1.506008e-03 1.426657e-03 1.552589e-03 1.209801e-03
  [11] 2.114903e-04 4.980992e-04 1.216711e-03 9.093551e-04 4.132601e-04
  [16] 1.579316e-03 1.027695e-03 1.421979e-03 1.582258e-03 1.457539e-03
  > dztgeom     (Rztgeom,      prob = pi/16)
   [1] 0.05289725 0.06582121 0.08190279 0.02745584 0.01425070 0.01773246
   [7] 0.19634954 0.08190279 0.12681315 0.08190279 0.12681315 0.12681315
  [13] 0.19634954 0.08190279 0.12681315 0.08190279 0.03416390 0.08190279
  [19] 0.15779640 0.15779640
  > dztbinom    (Rztbinom,     size = 55, prob = pi/16)
   [1] 0.09287368 0.09636862 0.06335468 0.13317760 0.13311093 0.13311093
   [7] 0.11849791 0.13311093 0.02881228 0.13311093 0.07063494 0.04717098
  [13] 0.11924710 0.13317760 0.09287368 0.04717098 0.11849791 0.03704401
  [19] 0.02881228 0.11924710
  > dlogarithmic(Rlogarithmic, prob = pi/16)
   [1] 2.055295e-18 8.982514e-01 4.782264e-14 2.670223e-04 1.263331e-06
   [6] 8.982514e-01 1.699916e-03 7.353248e-06 1.154347e-02 1.699916e-03
  [11] 8.895745e-15 5.815712e-17 1.699916e-03 8.982514e-01 5.956692e-60
  [16] 1.699916e-03 7.353248e-06 1.395489e-12 1.699916e-03 8.818562e-02
  > dzmpois       (Rzmpoisp0lt,        lambda = 6, p0 = 0.001)
   [1] 0.16086125 0.16086125 0.13405104 0.16086125 0.13405104 0.16086125
   [7] 0.16086125 0.13788107 0.16086125 0.16086125 0.16086125 0.13788107
  [13] 0.13405104 0.16086125 0.13405104 0.13788107 0.13405104 0.06894054
  [19] 0.10341081 0.13788107
  > dzmpois       (Rzmpoisp0gt,        lambda = 6, p0 = 0.010)
   [1] 0.15941205 0.06831945 0.10247918 0.13663890 0.15941205 0.13284338
   [7] 0.06831945 0.13284338 0.01117955 0.04099167 0.15941205 0.13284338
  [13] 0.13284338 0.15941205 0.15941205 0.15941205 0.13284338 0.15941205
  [19] 0.13284338 0.04428113
  > dzmnbinom     (Rzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01)
   [1] 0.36779812 0.36779812 0.08827155 0.36779812 0.29423850 0.36779812
   [7] 0.08827155 0.29423850 0.36779812 0.29423850 0.17654310 0.29423850
  [13] 0.36779812 0.17654310 0.29423850 0.36779812 0.29423850 0.36779812
  [19] 0.29423850 0.36779812
  > dzmnbinom     (Rzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40)
   [1] 0.22290795 0.40000000 0.40000000 0.22290795 0.40000000 0.40000000
   [7] 0.22290795 0.17832636 0.02353908 0.22290795 0.17832636 0.40000000
  [13] 0.10699582 0.22290795 0.40000000 0.40000000 0.10699582 0.10699582
  [19] 0.17832636 0.05349791
  > dzmgeom       (Rzmgeomp0lt,        prob = pi/16, p0 = 0.01)
   [1] 0.194386045 0.156218435 0.194386045 0.156218435 0.004729415 0.156218435
   [7] 0.065163000 0.042085788 0.194386045 0.156218435 0.156218435 0.194386045
  [13] 0.194386045 0.001972769 0.003054512 0.194386045 0.021844246 0.100894310
  [19] 0.125545017 0.125545017
  > dzmgeom       (Rzmgeomp0gt,        prob = pi/16, p0 = 0.40)
   [1] 0.02550654 0.03949273 0.07608789 0.40000000 0.03173835 0.07608789
   [7] 0.11780972 0.40000000 0.40000000 0.03173835 0.01647350 0.11780972
  [13] 0.11780972 0.02049834 0.40000000 0.40000000 0.02049834 0.02049834
  [19] 0.09467784 0.40000000
  > dzmbinom      (Rzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01)
   [1] 0.24858507 0.30523448 0.24858507 0.24858507 0.22714763 0.13665334
   [7] 0.30523448 0.24858507 0.05341988 0.22714763 0.22714763 0.24858507
  [13] 0.30523448 0.13665334 0.05341988 0.24858507 0.22714763 0.01522693
  [19] 0.22714763 0.24858507
  > dzmbinom      (Rzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12)
   [1] 0.27131954 0.20190901 0.20190901 0.27131954 0.27131954 0.20190901
   [7] 0.27131954 0.20190901 0.22096450 0.12146963 0.27131954 0.12000000
  [13] 0.22096450 0.12000000 0.12000000 0.04748434 0.27131954 0.12000000
  [19] 0.27131954 0.27131954
  > dzmlogarithmic(Rzmlogarithmicp0lt, prob = 0.99, p0 = 0.05)
   [1] 0.2042269801 0.0033322459 0.0004532564 0.0495403086 0.0043111747
   [6] 0.0152376857 0.2042269801 0.0500000000 0.2042269801 0.0667209544
  [11] 0.0039906388 0.1010923552 0.1010923552 0.0022849009 0.0323696377
  [16] 0.0020982672 0.0058245773 0.0323696377 0.2042269801 0.1010923552
  > dzmlogarithmic(Rzmlogarithmicp0gt, prob = 0.99, p0 = 0.55)
   [1] 0.5500000000 0.0099183875 0.0153329863 0.5500000000 0.0005000851
   [6] 0.5500000000 0.5500000000 0.0967390958 0.0039961403 0.0028940667
  [11] 0.5500000000 0.0153329863 0.5500000000 0.5500000000 0.5500000000
  [16] 0.5500000000 0.0234664620 0.0478858524 0.0153329863 0.0065960055
  > dpoisinvgauss(Rpoisinvgauss,    mean = 12,  dispersion = 0.1)
   [1] 0.0141875134 0.0534970428 0.0069390804 0.0141875134 0.0214217730
   [6] 0.0339451407 0.0534970428 0.0383803950 0.0758418721 0.0059021685
  [11] 0.0383803950 0.0534970428 0.0777438594 0.0301103042 0.0559189041
  [16] 0.0700915412 0.0559189041 0.0001072344 0.0493475070 0.0097718885
  > dpoisinvgauss(RpoisinvgaussInf, mean = Inf, dispersion = 1.1)
   [1] 0.0064949909 0.2596554485 0.0080356959 0.1750596677 0.0071973307
   [6] 0.2596554485 0.0646508555 0.2596554485 0.0000127472 0.1750596677
  [11] 0.0025580940 0.2596554485 0.0001730496 0.1027775188 0.0002624907
  [16] 0.0090499168 0.1750596677 0.0007972343 0.2596554485 0.0164681562
  > 
  > ## Check q<dist>(p<dist>(.)) identity
  > stopifnot(exprs = {
  +     All.eq(Rztpois,      qztpois     (Pztpois,      lambda = 12))
  +     All.eq(Rztnbinom,    qztnbinom   (Pztnbinom,    size = 7, prob = 0.01))
  +     All.eq(Rztgeom,      qztgeom     (Pztgeom,      prob = pi/16))
  +     All.eq(Rztbinom,     qztbinom    (Pztbinom,     size = 55, prob = pi/16))
  +     All.eq(Rlogarithmic, qlogarithmic(Plogarithmic, prob = 0.99))
  +     All.eq(Rzmpoisp0lt,        qzmpois       (Pzmpoisp0lt,        lambda = 6, p0 = 0.001))
  +     All.eq(Rzmpoisp0gt,        qzmpois       (Pzmpoisp0gt,        lambda = 6, p0 = 0.010))
  +     All.eq(Rzmnbinomp0lt,      qzmnbinom     (Pzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01))
  +     All.eq(Rzmnbinomp0gt,      qzmnbinom     (Pzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40))
  +     All.eq(Rzmgeomp0lt,        qzmgeom       (Pzmgeomp0lt,        prob = pi/16, p0 = 0.01))
  +     All.eq(Rzmgeomp0gt,        qzmgeom       (Pzmgeomp0gt,        prob = pi/16, p0 = 0.40))
  +     All.eq(Rzmbinomp0lt,       qzmbinom      (Pzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01))
  +     All.eq(Rzmbinomp0gt,       qzmbinom      (Pzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12))
  +     All.eq(Rzmlogarithmicp0lt, qzmlogarithmic(Pzmlogarithmicp0lt, prob = 0.99, p0 = 0.05))
  +     All.eq(Rzmlogarithmicp0gt, qzmlogarithmic(Pzmlogarithmicp0gt, prob = 0.99, p0 = 0.55))
  +     All.eq(Rpoisinvgauss,    qpoisinvgauss(Ppoisinvgauss,    mean = 12,  dispersion = 0.1))
  +     All.eq(RpoisinvgaussInf, qpoisinvgauss(PpoisinvgaussInf, mean = Inf, dispersion = 1.1))
  + })
  > 
  > ## Check q<dist>(p<dist>(.)) identity with upper tail
  > stopifnot(exprs = {
  +     All.eq(Rztpois,      qztpois     (1 - Pztpois,      lambda = 12, lower = FALSE))
  +     All.eq(Rztnbinom,    qztnbinom   (1 - Pztnbinom,    size = 7, prob = 0.01, lower = FALSE))
  +     All.eq(Rztgeom,      qztgeom     (1 - Pztgeom,      prob = pi/16, lower = FALSE))
  +     All.eq(Rztbinom,     qztbinom    (1 - Pztbinom,     size = 55, prob = pi/16, lower = FALSE))
  +     All.eq(Rlogarithmic, qlogarithmic(1 - Plogarithmic, prob = 0.99, lower = FALSE))
  +     All.eq(Rzmpoisp0lt,        qzmpois       (1 - Pzmpoisp0lt,        lambda = 6, p0 = 0.001, lower = FALSE))
  +     All.eq(Rzmpoisp0gt,        qzmpois       (1 - Pzmpoisp0gt,        lambda = 6, p0 = 0.010, lower = FALSE))
  +     All.eq(Rzmnbinomp0lt,      qzmnbinom     (1 - Pzmnbinomp0lt,      size = 7, prob = 0.8, p0 = 0.01, lower = FALSE))
  +     All.eq(Rzmnbinomp0gt,      qzmnbinom     (1 - Pzmnbinomp0gt,      size = 7, prob = 0.8, p0 = 0.40, lower = FALSE))
  +     All.eq(Rzmgeomp0lt,        qzmgeom       (1 - Pzmgeomp0lt,        prob = pi/16, p0 = 0.01, lower = FALSE))
  +     All.eq(Rzmgeomp0gt,        qzmgeom       (1 - Pzmgeomp0gt,        prob = pi/16, p0 = 0.40, lower = FALSE))
  +     All.eq(Rzmbinomp0lt,       qzmbinom      (1 - Pzmbinomp0lt,       size = 12, prob = pi/16, p0 = 0.01, lower = FALSE))
  +     All.eq(Rzmbinomp0gt,       qzmbinom      (1 - Pzmbinomp0gt,       size = 12, prob = pi/16, p0 = 0.12, lower = FALSE))
  +     All.eq(Rzmlogarithmicp0lt, qzmlogarithmic(1 - Pzmlogarithmicp0lt, prob = 0.99, p0 = 0.05, lower = FALSE))
  +     All.eq(Rzmlogarithmicp0gt, qzmlogarithmic(1 - Pzmlogarithmicp0gt, prob = 0.99, p0 = 0.55, lower = FALSE))
  +     All.eq(Rpoisinvgauss,    qpoisinvgauss(1 - Ppoisinvgauss,    mean = 12,  dispersion = 0.1, lower = FALSE))
  +     All.eq(RpoisinvgaussInf, qpoisinvgauss(1 - PpoisinvgaussInf, mean = Inf, dispersion = 1.1, lower = FALSE))
  + })
  Error: All.eq(Rzmlogarithmicp0lt, qzmlogarithmic(1 - Pzmlogarithmicp0lt,  .... is not TRUE
  Execution halted
** running tests for arch 'x64' ... [40s] OK
  Running 'betaint-tests.R' [1s]
  Running 'dpqr-tests.R' [39s]
  Running 'rmixture-tests.R' [1s]
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in 'inst/doc' ... OK
* checking re-building of vignette outputs ... [99s] OK
* checking PDF version of manual ... OK
* DONE
Status: 1 ERROR