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