Package 'RTMBdist'

Description

Smooth approximation to the absolute value function

Usage

abs_smooth(x, epsilon = 1e-06)

Arguments

Details

We approximate the absolute value here as

$\vert x \vert \approx \sqrt{x^2 + \epsilon}$

Value

Smooth absolute value of x.

Examples

abs(0)
abs_smooth(0, 1e-4)

Description

Density, distribution function, quantile function, and random generation for the Box–Cox Cole and Green distribution.

Usage

dbccg(x, mu = 1, sigma = 0.1, nu = 1, log = FALSE)

pbccg(q, mu = 1, sigma = 0.1, nu = 1, lower.tail = TRUE, log.p = FALSE)

qbccg(p, mu = 1, sigma = 0.1, nu = 1, lower.tail = TRUE, log.p = FALSE)

rbccg(n, mu = 1, sigma = 0.1, nu = 1)

Arguments

Details

dbccg and pbccg allow for automatic differentiation with RTMB. The parameterisation follows the BCCG family of the gamlss.dist package.

The density is

$f(x; \mu, \sigma, \nu) = \frac{x^{\nu-1}}{\mu^{\nu} \sigma \sqrt{2\pi}} \exp\!\left(-\tfrac{z^2}{2}\right) \Bigl[\Phi\!\left(\tfrac{1}{\sigma|\nu|}\right)\Bigr]^{-1}, \quad x > 0,$

where $z = [(x/\mu)^\nu - 1]/(\nu\sigma)$ for $\nu \neq 0$ and $z = \log(x/\mu)/\sigma$ for $\nu = 0$, and $\Phi$ is the standard normal CDF.

Value

dbccg gives the density, pbccg gives the distribution function, qbccg gives the quantile function, and rbccg generates random deviates.

References

Cole, T. J. and Green, P. J. (1992) Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in Medicine, 11, 1305-1319.

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC, doi:10.1201/9780429298547. An older version can be found in https://www.gamlss.com/.

See Also

bct, bcpe, gengamma

Examples

x <- rbccg(5, mu = 10, sigma = 0.2, nu = 0.5)
d <- dbccg(x, mu = 10, sigma = 0.2, nu = 0.5)
p <- pbccg(x, mu = 10, sigma = 0.2, nu = 0.5)
q <- qbccg(p, mu = 10, sigma = 0.2, nu = 0.5)

Description

Density, distribution function, quantile function, and random generation for the Box-Cox Power Exponential distribution.

Usage

dbcpe(x, mu = 5, sigma = 0.1, nu = 1, tau = 2, log = FALSE)

pbcpe(q, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE, log.p = FALSE)

qbcpe(p, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE, log.p = FALSE)

rbcpe(n, mu = 5, sigma = 0.1, nu = 1, tau = 2)

Arguments

Details

dbcpe and pbcpe allow for automatic differentiation with RTMB. The parameterisation follows the BCPE family of the gamlss.dist package.

The density is

$f(x; \mu, \sigma, \nu, \tau) = \frac{x^{\nu-1}}{\mu^{\nu} \sigma} \frac{f_T(z;\tau)}{F_T\!\left(1/(\sigma|\nu|);\tau\right)}, \quad x > 0,$

where $z = [(x/\mu)^\nu - 1]/(\nu\sigma)$ for $\nu \neq 0$ and $z = \log(x/\mu)/\sigma$ for $\nu = 0$, and $f_T(\cdot;\tau)$ and $F_T(\cdot;\tau)$ are the PDF and CDF of the power exponential (PE) distribution with shape $\tau$.

Value

dbcpe gives the density, pbcpe gives the distribution function, qbcpe gives the quantile function, and rbcpe generates random deviates.

References

Rigby, R. A. and Stasinopoulos, D. M. (2004) Smooth centile curves for skew and kurtotic data modelled using the Box-Cox Power Exponential distribution. Statistics in Medicine, 23, 3053-3076.

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC, doi:10.1201/9780429298547. An older version can be found in https://www.gamlss.com/.

See Also

bccg, bct, powerexp

Examples

x <- rbcpe(1, mu = 5, sigma = 0.1, nu = 1, tau = 1)
d <- dbcpe(x, mu = 5, sigma = 0.1, nu = 1, tau = 1)
p <- pbcpe(x, mu = 5, sigma = 0.1, nu = 1, tau = 1)
q <- qbcpe(p, mu = 5, sigma = 0.1, nu = 1, tau = 1)

Description

Density, distribution function, quantile function, and random generation for the Box–Cox t distribution.

Usage

dbct(x, mu = 5, sigma = 0.1, nu = 1, tau = 2, log = FALSE)

pbct(q, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE, log.p = FALSE)

qbct(p, mu = 5, sigma = 0.1, nu = 1, tau = 2, lower.tail = TRUE, log.p = FALSE)

rbct(n, mu = 5, sigma = 0.1, nu = 1, tau = 2)

Arguments

Details

dbct and pbct allow for automatic differentiation with RTMB. The parameterisation follows the BCT family of the gamlss.dist package.

The density is

$f(x; \mu, \sigma, \nu, \tau) = \frac{x^{\nu-1}}{\mu^{\nu} \sigma} \frac{f_t(z;\tau)}{F_t\!\left(1/(\sigma|\nu|);\tau\right)}, \quad x > 0,$

where $z = [(x/\mu)^\nu - 1]/(\nu\sigma)$ for $\nu \neq 0$ and $z = \log(x/\mu)/\sigma$ for $\nu = 0$, and $f_t(\cdot;\tau)$ and $F_t(\cdot;\tau)$ are the PDF and CDF of Student's $t$ distribution with $\tau$ degrees of freedom.

Value

dbct gives the density, pbct gives the distribution function, qbct gives the quantile function, and rbct generates random deviates.

References

Rigby, R. A. and Stasinopoulos, D. M. (2006) Using the Box-Cox t distribution in GAMLSS to model skewness and kurtosis. Statistical Modelling, 6(3), 209. doi:10.1191/1471082X06st122oa

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC, doi:10.1201/9780429298547. An older version can be found in https://www.gamlss.com/.

See Also

bccg, bcpe, skewt

Examples

x <- rbct(1, mu = 10, sigma = 0.2, nu = 0.5, tau = 4)
d <- dbct(x, mu = 10, sigma = 0.2, nu = 0.5, tau = 4)
p <- pbct(x, mu = 10, sigma = 0.2, nu = 0.5, tau = 4)
q <- qbct(p, mu = 10, sigma = 0.2, nu = 0.5, tau = 4)

Description

Probability mass function, distribution function, quantile function, and random generation for the Bell distribution.

Usage

dbell(x, theta, log = FALSE)

pbell(q, theta, lower.tail = TRUE, log.p = FALSE)

qbell(p, theta, lower.tail = TRUE, log.p = FALSE)

rbell(n, theta)

Arguments

Details

This implementation of dbell and pbell allows for automatic differentiation with RTMB with respect to theta.

The Bell distribution (Castellares et al. 2018) is a one-parameter distribution for overdispersed counts with probability mass function

$P(X = x;\, \theta) = \frac{e^{1 - e^{\theta}}\, \theta^{x}\, B_x}{x!}, \quad x = 0, 1, 2, \ldots,$

for $\theta > 0$, where $B_x$ is the $x$-th Bell number, i.e. the number of ways a set of $x$ elements can be partitioned into non-empty subsets.

Its mean and variance are

$E(X) = \theta e^{\theta}, \qquad \mathrm{Var}(X) = \theta e^{\theta} (1 + \theta),$

so the dispersion index is $\mathrm{Var}(X) / E(X) = 1 + \theta > 1$ and the distribution is always overdispersed. Note that, unlike the negative binomial or generalised Poisson distributions, the Bell distribution has no separate dispersion parameter: the amount of overdispersion is tied to the mean. As $\theta \to 0$ it approaches the Poisson distribution.

The distribution arises as a compound Poisson sum

$X = \sum_{i=1}^{N} Y_i, \qquad N \sim \mathrm{Pois}(e^{\theta} - 1), \quad Y_i \sim \mathrm{ztPois}(\theta),$

with the $Y_i$ independent of $N$, and rbell uses exactly this representation. It is therefore infinitely divisible, and a member of the one-parameter exponential family with natural parameter $\log\theta$ and sufficient statistic $x$.

The Bell numbers grow faster than any exponential and overflow double precision at $x = 219$, so $\log B_x$ is used throughout rather than $B_x$. As x is data, $\log B_x$ is constant with respect to the parameters and is cached across calls.

Neither the distribution function nor the quantile function has a closed form; both are obtained by summing the probability mass function, with qbell choosing its summation range automatically from the mean and variance.

See bell2 for the parameterisation by the mean.

Value

dbell gives the probability mass function, pbell gives the distribution function, qbell gives the quantile function, and rbell generates random deviates.

References

Castellares, F., Ferrari, S. L. P., and Lemonte, A. J. (2018). On the Bell distribution and its associated regression model for count data. Applied Mathematical Modelling 56, 172-185. doi:10.1016/j.apm.2017.12.014

Examples

set.seed(123)
x <- rbell(1, 1)
d <- dbell(x, 1)
p <- pbell(x, 1)
q <- qbell(p, 1)

# mean and variance
theta <- 0.8
xs <- 0:100
sum(xs * dbell(xs, theta)) # theta * exp(theta)

Description

Probability mass function, distribution function, quantile function, and random generation for the Bell distribution reparameterised in terms of its mean.

Usage

dbell2(x, mu, log = FALSE)

pbell2(q, mu, lower.tail = TRUE, log.p = FALSE)

qbell2(p, mu, lower.tail = TRUE, log.p = FALSE)

rbell2(n, mu)

Arguments

Details

This implementation of dbell2 and pbell2 allows for automatic differentiation with RTMB with respect to mu.

The Bell distribution has mean $\mu = \theta e^{\theta}$, which is a bijection from $\theta > 0$ to $\mu > 0$ and is inverted by the principal branch of the Lambert W function,

$\theta = W(\mu).$

Every positive mean therefore corresponds to exactly one $\theta$. All four functions simply apply this transformation and hand over to their bell counterparts.

In this parameterisation the variance is $\mu (1 + W(\mu))$, so the distribution is always overdispersed relative to the Poisson distribution, but the degree of overdispersion is determined by the mean rather than by a free parameter.

lambertW is AD-compatible to arbitrary order, so mu may be a parameter of a model fitted by Laplace approximation.

Value

dbell2 gives the probability mass function, pbell2 gives the distribution function, qbell2 gives the quantile function, and rbell2 generates random deviates.

References

Castellares, F., Ferrari, S. L. P., and Lemonte, A. J. (2018). On the Bell distribution and its associated regression model for count data. Applied Mathematical Modelling 56, 172-185. doi:10.1016/j.apm.2017.12.014

See Also

bell for the natural parameterisation.

Examples

set.seed(123)
x <- rbell2(1, 3)
d <- dbell2(x, 3)
p <- pbell2(x, 3)
q <- qbell2(p, 3)

# the two parameterisations agree
all.equal(dbell2(0:5, 3), dbell(0:5, lambertW(3)))

Description

Density, distribution function, quantile function, and random generation for the beta distribution reparameterised in terms of mean and concentration.

Usage

dbeta(x, shape1, shape2, log = FALSE, eps = 0)

dbeta2(x, mu, phi, log = FALSE, eps = 0)

pbeta2(q, mu, phi, lower.tail = TRUE, log.p = FALSE)

qbeta2(p, mu, phi, lower.tail = TRUE, log.p = FALSE)

rbeta2(n, mu, phi)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

Currently, dbeta masks RTMB::dbeta because the latter has a numerically unstable gradient.

dbeta and dbeta2 compute the beta density. dbeta2 reparameterises by mean $\mu$ and concentration $\phi$:

$f(x; \mu, \phi) = \frac{x^{\mu\phi - 1}(1-x)^{(1-\mu)\phi - 1}}{B(\mu\phi,\, (1-\mu)\phi)}, \quad x \in (0,1).$

Value

dbeta2 gives the density, pbeta2 gives the distribution function, qbeta2 gives the quantile function, and rbeta2 generates random deviates.

Examples

set.seed(123)
x <- rbeta2(1, 0.5, 1)
d <- dbeta2(x, 0.5, 1)
p <- pbeta2(x, 0.5, 1)
q <- qbeta2(p, 0.5, 1)

Description

Density, distribution function and random generation for the beta-binomial distribution.

Usage

dbetabinom(x, size, shape1, shape2, log = FALSE)

pbetabinom(q, size, shape1, shape2, lower.tail = TRUE, log.p = FALSE)

rbetabinom(n, size, shape1, shape2)

Arguments

Details

This implementation of dbetabinom allows for automatic differentiation with RTMB.

$P(X = k;\, n, a, b) = \binom{n}{k} \frac{B(k+a,\, n-k+b)}{B(a,\, b)}, \quad k = 0, 1, \ldots, n.$

The distribution function has no closed form and is computed by summing the probability mass function over $0, \ldots, q$. It is AD-compatible in the parameters, while q and size must be numeric data. This is also what one-step-ahead (OSA) residuals via RTMB::oneStepPredict need, so these are supported, e.g. with method = "cdf" and discrete = TRUE.

Value

dbetabinom gives the density, pbetabinom gives the distribution function, and rbetabinom generates random samples.

Examples

set.seed(123)
x <- rbetabinom(1, 10, 2, 5)
d <- dbetabinom(x, 10, 2, 5)
p <- pbetabinom(x, 10, 2, 5)

Description

Density, distribution function, quantile function, and random generation for the Beta prime distribution.

Usage

dbetaprime(x, shape1, shape2, log = FALSE)

pbetaprime(q, shape1, shape2, lower.tail = TRUE, log.p = FALSE)

qbetaprime(p, shape1, shape2, lower.tail = TRUE, log.p = FALSE)

rbetaprime(n, shape1, shape2)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

If $X \sim \text{Beta}(\alpha, \beta)$, then $\frac{X}{1-X} \sim \text{Betaprime}(\alpha, \beta)$

$f(x;\, a, b) = \frac{x^{a-1}}{(1+x)^{a+b}\, B(a,b)}, \quad x > 0.$

Value

dbetaprime gives the density, pbetaprime gives the distribution function, qbetaprime gives the quantile function, and rbetaprime generates random deviates.

Examples

x <- rbetaprime(1, 2, 1)
d <- dbetaprime(x, 2, 1)
p <- pbetaprime(x, 2, 1)
q <- qbetaprime(p, 2, 1)

Description

Probability mass function, distribution function and random generation for the beta-negative binomial distribution.

Usage

dbnbinom(x, size, shape1, shape2, log = FALSE)

pbnbinom(q, size, shape1, shape2, lower.tail = TRUE, log.p = FALSE)

rbnbinom(n, size, shape1, shape2)

Arguments

Details

dbnbinom allows for automatic differentiation with RTMB.

The beta-negative binomial arises by giving the success probability of a negative binomial a beta prior, in the same way that the beta-binomial does for the binomial:

$P(X = k;\, r, a, b) = \frac{\Gamma(k + r)}{k!\, \Gamma(r)} \frac{B(a + r,\, b + k)}{B(a,\, b)}, \quad k = 0, 1, 2, \ldots$

The extra beta layer gives a much heavier tail than the negative binomial: the mean $rb / (a - 1)$ exists only for $a > 1$ and the variance

$\frac{r b (r + a - 1)(b + a - 1)}{(a - 2)(a - 1)^2}$

only for $a > 2$. As $a \to \infty$ with $b/(a+b)$ held fixed the distribution collapses to the negative binomial.

Note that the three parameters are only weakly identified: the likelihood has a long ridge along which size and shape2 trade off against each other, so fitting all three at once needs either a lot of data or a restriction. The mean parameterisation is usually the more stable one to estimate in.

The distribution function has no closed form and is computed by summing the probability mass function over $0, \ldots, q$, so its cost grows with the largest q. It is AD-compatible in the parameters, while q must be numeric data. This is also what one-step-ahead (OSA) residuals via RTMB::oneStepPredict need, so these are supported, e.g. with method = "cdf" and discrete = TRUE.

Value

dbnbinom gives the probability mass function, pbnbinom gives the distribution function, and rbnbinom generates random deviates.

References

Johnson, N. L., Kemp, A. W. and Kotz, S. (2005) Univariate Discrete Distributions, 3rd edition, Wiley, doi:10.1002/0471715816.

See Also

bnbinom2, betabinom, nbinom2

Examples

set.seed(123)
x <- rbnbinom(5, size = 3, shape1 = 4, shape2 = 2)
d <- dbnbinom(x, size = 3, shape1 = 4, shape2 = 2)
p <- pbnbinom(x, size = 3, shape1 = 4, shape2 = 2)

Description

Probability mass function, distribution function and random generation for the beta-negative binomial distribution reparameterised in terms of its mean.

Usage

dbnbinom2(x, mu, sigma, nu, log = FALSE)

pbnbinom2(q, mu, sigma, nu, lower.tail = TRUE, log.p = FALSE)

rbnbinom2(n, mu, sigma, nu)

Arguments

Details

dbnbinom2 allows for automatic differentiation with RTMB. The parameterisation follows the BNB family of the gamlss.dist package, in which $\mu$ is exactly the mean.

Writing $r$, $a$ and $b$ for the arguments size, shape1 and shape2 of dbnbinom, the reparameterisation is

$a = \frac{1}{\sigma} + 1, \qquad b = \frac{\mu\nu}{\sigma}, \qquad r = \frac{1}{\nu}.$

This gives $E(X) = \mu$ for every admissible $\sigma$ and $\nu$, where the original parameterisation needs $a > 1$ for a mean to exist at all. The variance is

$\mathrm{Var}(X) = \frac{\mu (\sigma + \nu)(\mu\nu + 1)}{\nu (1 - \sigma)},$

which is finite for $\sigma < 1$ and increases without bound as $\sigma \to 1$. Both $\sigma$ and $\nu$ add dispersion beyond the negative binomial, which is recovered as $\sigma \to 0$.

Because the mean is pinned to a single parameter, this parameterisation is the more stable of the two to estimate in; see bnbinom for the identifiability problem it avoids.

The distribution function has no closed form and is computed by summing the probability mass function over $0, \ldots, q$, so its cost grows with the largest q. It is AD-compatible in the parameters, while q must be numeric data. This is also what one-step-ahead (OSA) residuals via RTMB::oneStepPredict need, so these are supported, e.g. with method = "cdf" and discrete = TRUE.

Value

dbnbinom2 gives the probability mass function, pbnbinom2 gives the distribution function, and rbnbinom2 generates random deviates.

References

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC, doi:10.1201/9780429298547. An older version can be found in https://www.gamlss.com/.

See Also

bnbinom, nbinom2, betabinom

Examples

set.seed(123)
x <- rbnbinom2(5, mu = 4, sigma = 0.4, nu = 0.5)
d <- dbnbinom2(x, mu = 4, sigma = 0.4, nu = 0.5)
p <- pbnbinom2(x, mu = 4, sigma = 0.4, nu = 0.5)

Description

Construct a function that computes the log density or CDF of the bivariate Clayton copula, intended to be used with dcopula.

Usage

cclayton(theta)

Cclayton(theta)

Arguments

Details

The Clayton copula density is

$c(u,v;\theta) = (1+\theta) (uv)^{-(1+\theta)} \left( u^{-\theta} + v^{-\theta} - 1 \right)^{-(2\theta+1)/\theta}, \quad \theta > 0.$

Value

A function of two arguments (u,v) returning log copula density (cclayton) or copula CDF (Cclayton).

See Also

cgaussian(), cgumbel(), cfrank()

Examples

x <- c(0.5, 1); y <- c(0.2, 0.8)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dbeta(y, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pbeta(y, 2, 1)
dcopula(d1, d2, p1, p2, copula = cclayton(2), log = TRUE)

# CDF version (for discrete copulas)
Cclayton(1.5)(0.5, 0.4)

Description

Turns a copula for two linear variables into a circular-linear copula that is symmetric in the sign of the angle, intended to be used with dcopula to join a turning angle and a step length.

Usage

cfold(copula)

Arguments

Details

The copula density is

$c(u, v) = c_0\bigl(1 - |2u - 1|, \, v\bigr),$

where $u$ is the distribution function of the circular margin, $v$ that of the linear margin and $c_0$ is the density of the copula passed as copula, e.g. cgaussian, cclayton, cgumbel or cfrank. This is a copula for any $c_0$. It is the rectangular patchwork copula of Hodel and Fieberg (2022), with the linear copula in the rectangle $u \le 1/2$ and its mirror image in the rectangle $u > 1/2$.

The map $u \mapsto 1 - |2u - 1|$ folds the circle at $u = 1/2$. For a circular margin that is symmetric about its mean direction and whose distribution function is cut open at the antipode of that mean direction, as by default in pwrpcauchy and pvm, $u = 1/2$ is the mean direction. $1 - |2u - 1|$ is then the distribution function of the angular distance to the mean direction, counted from the antipode: it is 1 for angles at the mean direction and 0 for angles opposite to it. In words, $c_0$ links the straightness of a step to its length. A copula with positive dependence, such as cgaussian(rho) with rho > 0, makes long steps straight and lets short steps turn in either direction, the pattern most common in movement data. Unlike cjw, the resulting copula is symmetric, $c(u, v) = c(1 - u, v)$.

In dcopula, the circular margin must be the first one, i.e. d1 and p1 belong to the angle.

Value

Function of two arguments (u, v) returning the log copula density, with u for the circular and v for the linear margin.

References

Hodel, F. H. and Fieberg, J. R. (2022) Circular-linear copulae for animal movement data. Methods in Ecology and Evolution, 13, doi:10.1111/2041-210X.13821.

Durante, F., Saminger-Platz, S. and Sarkoci, P. (2009) Rectangular patchwork for bivariate copulas and tail dependence. Communications in Statistics - Theory and Methods, 38, 2515-2527, doi:10.1080/03610920802571203.

See Also

dcopula(), cjw(), vm, wrpcauchy

Examples

# turning angles with von Mises margin, step lengths with Weibull margin
angle <- c(-2, -0.3, 0.1, 1.5); step <- c(0.2, 1.1, 2.4, 0.6)
d1 <- dvm(angle, 0, 2, log = TRUE); p1 <- pvm(angle, 0, 2)
d2 <- dweibull(step, 2, 1, log = TRUE); p2 <- pweibull(step, 2, 1)
dcopula(d1, d2, p1, p2, copula = cfold(cgaussian(0.5)), log = TRUE)

# wrapped Cauchy margin, which also allows for automatic differentiation
d1 <- dwrpcauchy(angle, 0, 0.5, log = TRUE); p1 <- pwrpcauchy(angle, 0, 0.5)
dcopula(d1, d2, p1, p2, copula = cfold(cclayton(2)), log = TRUE)

# the copula is symmetric in the sign of the angle
cop <- cfold(cclayton(2))
cop(c(0.2, 0.8), 0.9)

Description

Returns a function computing the log density of the bivariate Frank copula, intended to be used with dcopula.

Usage

cfrank(theta)

Cfrank(theta)

Arguments

Details

The Frank copula density is

$c(u,v;\theta) = \frac{\theta (1-e^{-\theta}) e^{-\theta(u+v)}} {\left[(1 - e^{-\theta}) - (1 - e^{-\theta u})(1 - e^{-\theta v}) \right]^2}, \quad \theta \ne 0.$

Value

A function of two arguments (u, v) returning either the log copula density (cfrank) or the copula CDF (Cfrank).

See Also

cgaussian(), cclayton(), cgumbel()

Examples

x <- c(0.5, 1); y <- c(1, 2)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2)
dcopula(d1, d2, p1, p2, copula = cfrank(2), log = TRUE)

Description

Returns a function computing the log density of the bivariate Gaussian copula, intended to be used with dcopula.

Usage

cgaussian(rho = 0)

Arguments

Value

Function of two arguments (u,v) returning log copula density.

The Gaussian copula density is

$c(u,v;\rho) = \frac{1}{\sqrt{1-\rho^2}} \exp\left\{-\frac{1}{2(1-\rho^2)} (z_1^2 - 2 \rho z_1 z_2 + z_2^2) + \frac{1}{2}(z_1^2 + z_2^2) \right\},$

where $z_1 = \Phi^{-1}(u)$, $z_2 = \Phi^{-1}(v)$, and $-1 < \rho < 1$.

See Also

cclayton(), cgumbel(), cfrank()

Examples

x <- c(0.5, 1); y <- c(1, 2)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2)
dcopula(d1, d2, p1, p2, copula = cgaussian(0.5), log = TRUE)

Description

Returns a function computing the log density of the multivariate Gaussian copula, parameterised by the inverse correlation matrix.

Usage

cgmrf(Q)

Arguments

Details

Caution: Parameterising the inverse correlation directly is difficult, as inverting it needs to yield a positive definite matrix with unit diagonal. Hence we still advise parameterising the correaltion matrix R and computing its inverse. This function is useful when you need access to the precision (i.e. inverse correlation) in your likelihood function.

Value

Function with matrix argument U returning log copula density.

See Also

cmvgauss()

Examples

x <- c(0.5, 1); y <- c(1, 2); z <- c(0.2, 0.8)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE); d3 <- dbeta(z, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2); p3 <- pbeta(z, 2, 1)
R <- matrix(c(1,0.5,0.3,0.5,1,0.4,0.3,0.4,1), nrow = 3)

## Based on correlation matrix
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cmvgauss(R), log = TRUE)

## Based on precision matrix
Q <- solve(R)
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cgmrf(Q), log = TRUE)

## Parameterisation inside a model
# using RTMB::unstructured to get a valid correlation matrix
library(RTMB)
d <- 5 # dimension
cor_func <- unstructured(d)
npar <- length(cor_func$parms())
R <- cor_func$corr(rep(0.1, npar))

Description

Construct functions that compute either the log density or the CDF of the bivariate Gumbel copula, intended for use with dcopula.

Usage

cgumbel(theta)

Cgumbel(theta)

Arguments

Details

The Gumbel copula density

$c(u,v;\theta) = \frac{e^{-A}}{uv} \, (\tilde u \tilde v)^{\theta - 1} S^{1/\theta - 2} (A + \theta - 1),$

where $\tilde u = -\log u$, $\tilde v = -\log v$, $S = \tilde u^\theta + \tilde v^\theta$ and $A = S^{1/\theta}$, so that $e^{-A}$ is the copula CDF.

Value

A function of two arguments (u, v) returning either the log copula density (cgumbel) or the copula CDF (Cgumbel).

See Also

cgaussian(), cclayton(), cfrank()

Examples

x <- c(0.5, 1); y <- c(0.2, 0.4)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dbeta(y, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pbeta(y, 2, 1)
dcopula(d1, d2, p1, p2, copula = cgumbel(1.5), log = TRUE)

# CDF version (for discrete copulas)
Cgumbel(1.5)(0.5, 0.4)

Description

Returns a function computing the log density of the circular-linear copula of Johnson and Wehrly (1978), intended to be used with dcopula to join a circular and a linear margin, such as the turning angles and step lengths of an animal track.

Usage

cjw(g, ..., q = 1)

Arguments

Details

The copula density is

$c(u, v) = 2\pi \, g\bigl(2\pi(u - q v)\bigr),$

where $u$ is the distribution function of the circular margin, $v$ that of the linear margin, $g$ is a density on the circle, the binding density, and $q = \pm 1$. The joint density of an angle $\theta$ and a step length $s$ is then

$f(\theta, s) = 2\pi \, g\bigl(2\pi(F_1(\theta) - q F_2(s))\bigr) f_1(\theta) f_2(s).$

Because $g$ is periodic, $c$ is a copula for any circular density and any integer $q \neq 0$; only $q = \pm 1$ is allowed here. Any circular density of this package that allows for automatic differentiation can be used for $g$, in particular dvm and dwrpcauchy. Its concentration controls the strength of the dependence, and q = -1 reverses its direction.

In dcopula, the circular margin must be the first one, i.e. d1 and p1 belong to the angle. Its distribution function can be pwrpcauchy, which allows for automatic differentiation. Where the circle is cut open for this distribution function only shifts $u$ by a constant, which the location of $g$ absorbs. The likelihood is therefore the same for any origin, but the meaning of the location of $g$ is not.

The dependence of this copula is a helix: as $v$ goes from 0 to 1, the angle it favours turns once around the circle. The copula is hence not symmetric in the sign of the angle, $c(u, v) \neq c(1 - u, v)$. With a turning angle margin centred at 0 and the location of $g$ at 0, for example, medium steps are straight, short steps turn one way and long steps the other. It therefore cannot capture the pattern most common in movement data, where long steps are straight and short steps turn in either direction. For this, see cfold.

Random pairs from the copula are obtained by drawing $v$ uniformly and $z$ from $g$, and setting $u = (z / (2\pi) + q v) \bmod 1$; see the examples.

Value

Function of two arguments (u, v) returning the log copula density, with u for the circular and v for the linear margin.

References

Johnson, R. A. and Wehrly, T. E. (1978) Some angular-linear distributions and related regression models. Journal of the American Statistical Association, 73, 602-606, doi:10.1080/01621459.1978.10480062.

Hodel, F. H. and Fieberg, J. R. (2022) Circular-linear copulae for animal movement data. Methods in Ecology and Evolution, 13, doi:10.1111/2041-210X.13821.

See Also

dcopula(), cfold(), vm, wrpcauchy

Examples

# turning angles with wrapped Cauchy margin, step lengths with Weibull margin
angle <- c(-2, -0.3, 0.1, 1.5); step <- c(0.2, 1.1, 2.4, 0.6)
d1 <- dwrpcauchy(angle, 0, 0.5, log = TRUE); p1 <- pwrpcauchy(angle, 0, 0.5)
d2 <- dweibull(step, 2, 1, log = TRUE); p2 <- pweibull(step, 2, 1)

# von Mises binding density
dcopula(d1, d2, p1, p2, copula = cjw(dvm, mu = 0, kappa = 2), log = TRUE)

# wrapped Cauchy binding density, dependence in the other direction
dcopula(d1, d2, p1, p2, copula = cjw(dwrpcauchy, mu = 0, rho = 0.7, q = -1), log = TRUE)

# simulation from the joint distribution with von Mises binding density
n <- 1000
v <- runif(n)
z <- rvm(n, mu = 0, kappa = 2)
u <- (z / (2 * pi) + v) %% 1
angle <- qwrpcauchy(u, 0, 0.5)
step <- qweibull(v, 2, 1)

Description

Returns a function computing the log density of the multivariate Gaussian copula, intended to be used with dmvcopula.

Usage

cmvgauss(R)

Arguments

Value

Function with matrix argument U returning log copula density.

See Also

cgmrf()

Examples

x <- c(0.5, 1); y <- c(1, 2); z <- c(0.2, 0.8)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE); d3 <- dbeta(z, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2); p3 <- pbeta(z, 2, 1)
R <- matrix(c(1,0.5,0.3,0.5,1,0.4,0.3,0.4,1), nrow = 3)

## Based on correlation matrix
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cmvgauss(R), log = TRUE)

## Based on precision matrix
Q <- solve(R)
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cgmrf(Q), log = TRUE)

## Parameterisation inside a model
# using RTMB::unstructured to get a valid correlation matrix
library(RTMB)
d <- 5 # dimension
cor_func <- unstructured(d)
npar <- length(cor_func$parms())
R <- cor_func$corr(rep(0.1, npar))

Description

Probability mass function, distribution function, quantile function, and random generation for the Conway-Maxwell-binomial (CMB) distribution.

Usage

dcombinom(x, size, prob, nu = 1, log = FALSE)

pcombinom(q, size, prob, nu = 1, lower.tail = TRUE, log.p = FALSE)

qcombinom(p, size, prob, nu = 1, lower.tail = TRUE, log.p = FALSE)

rcombinom(n, size, prob, nu = 1)

Arguments

Details

This implementation of dcombinom and pcombinom allows for automatic differentiation with RTMB, including differentiation with respect to x, such that one-step-ahead (OSA) residuals are supported.

The CMB distribution generalises the binomial distribution by an additional dispersion parameter $\nu$ in the same way that the Conway-Maxwell-Poisson distribution generalises the Poisson distribution. Its probability mass function is

$P(X = x;\, n, p, \nu) = \frac{1}{Z(n, p, \nu)} \binom{n}{x}^{\nu} p^x (1-p)^{n-x}, \quad x = 0, 1, \ldots, n,$

with normalising constant

$Z(n, p, \nu) = \sum_{k=0}^{n} \binom{n}{k}^{\nu} p^k (1-p)^{n-k}.$

For $\nu = 1$ this reduces to the binomial distribution. Values $\nu > 1$ give under-dispersion and $\nu < 1$ over-dispersion relative to a binomial distribution with the same mean. As the support is finite, $Z$ converges for every real $\nu$, so $\nu$ is not restricted to be positive; negative values yield strongly over-dispersed, U-shaped distributions.

The distribution arises as the sum of $n$ exchangeable, possibly associated Bernoulli variables, where $\nu$ controls the association. Note that prob is not the mean divided by size unless $\nu = 1$; the mean has no closed form for general $\nu$ and must be obtained by summation over the support.

Value

dcombinom gives the probability mass function, pcombinom gives the distribution function, qcombinom gives the quantile function, and rcombinom generates random deviates.

References

Shmueli, G., Minka, T. P., Kadane, J. B., Borle, S., and Boatwright, P. (2005). A useful distribution for fitting discrete data: revival of the Conway-Maxwell-Poisson distribution. Journal of the Royal Statistical Society: Series C 54(1), 127-142.

Kadane, J. B. (2016). Sums of possibly associated Bernoulli variables: the Conway-Maxwell-binomial distribution. Bayesian Analysis 11(2), 403-420.

https://en.wikipedia.org/wiki/Conway-Maxwell-binomial_distribution

Examples

set.seed(123)
x <- rcombinom(1, size = 10, prob = 0.4, nu = 1.5)
d <- dcombinom(x, 10, 0.4, 1.5)
p <- pcombinom(x, 10, 0.4, 1.5)
q <- qcombinom(p, 10, 0.4, 1.5)

# nu = 1 recovers the binomial distribution
all.equal(dcombinom(0:10, 10, 0.3, 1), dbinom(0:10, 10, 0.3))

Description

Computes the joint density (or log-density) of a bivariate distribution constructed from two arbitrary margins combined with a specified copula.

Usage

dcopula(d1, d2, p1, p2, copula = cgaussian(0), log = FALSE)

Arguments

Details

The joint density is

$f(x,y) = c(F_1(x), F_2(y)) \, f_1(x) f_2(y),$

where $F_i$ are the marginal CDFs, $f_i$ are the marginal densities, and $c$ is the copula density.

The marginal densities d1, d2 and CDFs p1, p2 must be differentiable for automatic differentiation (AD) to work.

Available copula constructors are:

• cgaussian (Gaussian copula)

• cclayton (Clayton copula)

• cgumbel (Gumbel copula)

• cfrank (Frank copula)

Circular-linear copulas, e.g. for turning angles and step lengths, with the circular margin as the first one:

• cjw (Johnson-Wehrly copula with a circular binding density such as dvm or dwrpcauchy)

• cfold (folds any of the copulas above into one that is symmetric in the sign of the angle)

Value

Joint density (or log-density) under the bivariate copula.

See Also

ddcopula(), dmvcopula()

Examples

# Normal + Exponential margins with Gaussian copula
x <- c(0.5, 1); y <- c(1, 2)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2)
dcopula(d1, d2, p1, p2, copula = cgaussian(0.5), log = TRUE)

# Normal + Beta margins with Clayton copula
x <- c(0.5, 1); y <- c(0.2, 0.8)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dbeta(y, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pbeta(y, 2, 1)
dcopula(d1, d2, p1, p2, copula = cclayton(2), log = TRUE)

# Normal + Beta margins with Gumbel copula
x <- c(0.5, 1); y <- c(0.2, 0.4)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dbeta(y, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pbeta(y, 2, 1)
dcopula(d1, d2, p1, p2, copula = cgumbel(1.5), log = TRUE)

# Normal + Exponential margins with Frank copula
x <- c(0.5, 1); y <- c(1, 2)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2)
dcopula(d1, d2, p1, p2, copula = cfrank(2), log = TRUE)

Description

Computes the joint probability mass function of two discrete margins combined with a copula CDF.

Usage

ddcopula(d1, d2, p1, p2, Copula, log = FALSE)

Arguments

Details

The joint probability mass function for two discrete margins is

$\Pr(Y_1 = y_1, Y_2 = y_2) = C(F_1(y_1), F_2(y_2)) - C(F_1(y_1-1), F_2(y_2)) - C(F_1(y_1), F_2(y_2-1)) + C(F_1(y_1-1), F_2(y_2-1)),$

where $F_i$ are the marginal CDFs, and $C$ is the copula CDF.

Available copula CDF constructors are:

• Cclayton (Clayton copula)

• Cgumbel (Gumbel copula)

• Cfrank (Frank copula)

Value

Joint probability (or log-probability) under chosen copula

See Also

dcopula(), dmvcopula()

Examples

x <- c(3,5); y <- c(2,4)
d1 <- dpois(x, 4); d2 <- dpois(y, 3)
p1 <- ppois(x, 4); p2 <- ppois(y, 3)
ddcopula(d1, d2, p1, p2, Copula = Cclayton(2), log = FALSE)

Description

Density and and random generation for the Dirichlet distribution.

Usage

ddirichlet(x, alpha, log = FALSE)

rdirichlet(n, alpha)

Arguments

Details

This implementation of ddirichlet allows for automatic differentiation with RTMB.

$f(\mathbf{x};\,\boldsymbol{\alpha}) = \frac{\Gamma\!\left(\sum_i \alpha_i\right)}{\prod_i \Gamma(\alpha_i)} \prod_i x_i^{\alpha_i - 1},$

where $\mathbf{x}$ lies on the unit simplex and $\alpha_i > 0$.

Value

ddirichlet gives the density, rdirichlet generates random deviates.

Examples

# single alpha
alpha <- c(1,2,3)
x <- rdirichlet(1, alpha)
d <- ddirichlet(x, alpha)
# vectorised over alpha
alpha <- rbind(alpha, 2*alpha)
x <- rdirichlet(2, alpha)

Description

Density and and random generation for the Dirichlet-multinomial distribution.

Usage

ddirmult(x, size, alpha, log = FALSE)

rdirmult(n, size, alpha)

Arguments

Details

This implementation of ddirmult allows for automatic differentiation with RTMB.

$P(\mathbf{x};\,\boldsymbol{\alpha}, n) = \frac{\Gamma(\alpha_0)\,\Gamma(n+1)}{\Gamma(n+\alpha_0)} \prod_i \frac{\Gamma(x_i + \alpha_i)}{\Gamma(\alpha_i)\,\Gamma(x_i+1)},$

where $\alpha_0 = \sum_i \alpha_i$ and $n = \sum_i x_i$.

Value

ddirmult gives the density and rdirmult generates random samples.

Examples

# single alpha
alpha <- c(1,2,3)
size <- 10
x <- rdirmult(1, size, alpha)
d <- ddirmult(x, size, alpha)
# vectorised over alpha and size
alpha <- rbind(alpha, 2*alpha)
size <- c(size, 3*size)
x <- rdirmult(2, size, alpha)

Description

Computes the joint density (or log-density) of a distribution constructed from any number of arbitrary margins combined with a specified copula.

Usage

dmvcopula(D, P, copula = cmvgauss(diag(ncol(D))), log = FALSE)

Arguments

Details

The joint density is

$f(x_1, \dots, x_d) = c(F_1(x_1), \dots, F_d(x_d)) \, f_1(x_1) \dots f_d(x_d),$

where $F_i$ are the marginal CDFs, $f_i$ are the marginal densities, and $c$ is the copula density.

The marginal densities d_1, ..., d_d and CDFs p_1, ..., p_d must be differentiable for automatic differentiation (AD) to work.

Available multivariate copula constructors are:

• cmvgauss (Multivariate Gaussian copula)

• cgmrf (Multivariate Gaussian copula parameterised by precision (inverse correlation) matrix)

Value

Joint density (or log-density) under the chosen copula.

See Also

dcopula(), ddcopula()

Examples

x <- c(0.5, 1); y <- c(1, 2); z <- c(0.2, 0.8)
d1 <- dnorm(x, 1, log = TRUE); d2 <- dexp(y, 2, log = TRUE); d3 <- dbeta(z, 2, 1, log = TRUE)
p1 <- pnorm(x, 1); p2 <- pexp(y, 2); p3 <- pbeta(z, 2, 1)
R <- matrix(c(1,0.5,0.3,0.5,1,0.4,0.3,0.4,1), nrow = 3)

### Multivariate Gaussian copula
## Based on correlation matrix
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cmvgauss(R), log = TRUE)

## Based on precision matrix
Q <- solve(R)
dmvcopula(cbind(d1, d2, d3), cbind(p1, p2, p3), copula = cgmrf(Q), log = TRUE)

## Parameterisation inside a model
# using RTMB::unstructured to get a valid correlation matrix
library(RTMB)
d <- 5 # dimension
cor_func <- unstructured(d)
npar <- length(cor_func$parms())
R <- cor_func$corr(rep(0.1, npar))

Description

AD-compatible error function and complementary error function

Usage

erf(x)

erfc(x)

Arguments

Value

erf(x) returns the error function and erfc(x) returns the complementary error function.

Examples

erf(1)
erfc(1)