Package 'dgumbel'

The Gumbel Distribution and Derivatives

Description

Density function, distribution function, quantile function and random generation, and their gradient functions for the Gumbel distribution with location and scale parameters.

Usage

dgumbel(x, location=0, scale=1, log = FALSE, grad=FALSE) 
pgumbel(q, location=0, scale=1, lower.tail = TRUE, log.p = FALSE, grad=FALSE) 
qgumbel(p, location=0, scale=1, lower.tail = TRUE, grad=FALSE)
rgumbel(n, location=0, scale=1)

Arguments

x, q	Vector of quantiles.
p	Vector of probabilities.
n	Number of observations.
location, scale	Location and scale parameters.
log, log.p	Logical; if TRUE, probabilities p are given as log(p).
lower.tail	Logical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
grad	Logical; if TRUE, the gradient w.r.t. parameters location and scale is given instead of function value.

Details

The Gumbel distribution function with parameters $\code{location} = a$ and $\code{scale} = b$ is

$G(z) = \exp\left\{-\exp\left[-\left(\frac{z-a}{b}\right) \right]\right\}$

for all real $z$, where $b > 0$. Gradients are exact numerical derivatives implemented using automatic differentiation. dgumbel builds on the Eigen linear algebra library, Adept for automatic differentiation and RcppEigen for bindings to R and loading Eigen.

Value

dgumbel gives the density function, pgumbel gives the distribution function, qgumbel gives the quantile function, and rgumbel generates random deviates. If grad=TRUE is supplied, then the gradient is returned instead of the objective function.

Examples

dgumbel(-1:2, -1, 0.5)
pgumbel(-1:2, -1, 0.5)
qgumbel(seq(0.9, 0.6, -0.1), 2, 0.5)
rgumbel(6, -1, 0.5)
p <- (1:9)/10
pgumbel(qgumbel(p, -1, 2), -1, 2)
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

## Random number generation
loc = .5
scale = 3.2
n <- 1000
x <- rgumbel(n, loc, scale)

## The density
hist(x, freq=FALSE)
xs <- sort(x)
fx <- dgumbel(xs, loc, scale)
points(xs,fx, type="l", col=2, lwd=2)

## The distribution
edf <- sapply(xs, function(x){sum(xs<=x)/n})
plot(xs, edf)
Fx <- pgumbel(xs, loc, scale)
points(xs, Fx, type="l", col=2, lwd=2) 

## The quantile function
q <- qgumbel(0.6, loc, scale)
polygon(c(xs[xs <= q], q), c(Fx[xs<=q], 0), col=3)

## Negative log likelihood: Objective and gradient
nll <- function(par, data) -sum(dgumbel(data, par[1], par[2], log=TRUE))
dnll <- function(par, data) -rowSums(dgumbel(data, par[1], par[2], log=TRUE, grad=TRUE))

## Parameter estimation
par_start <- c(3,1)
opt <- nlminb(par_start, objective=nll, gradient=dnll, data=x, control = list(trace=5))
opt$convergence
opt$par

---