Package 'bgev'

Bimodal GEV (generalized extreme value) distribution

Description

Functions to compute the density, distribution function, quantile function, and to generate random variates for the BGEV (bimodal generalized extreme value)

Usage

dbgev(x, mu = 1, sigma = 1, xi = 0.3, delta = 2)

pbgev(q, mu = 1, sigma = 1, xi = 0.3, delta = 2)

qbgev(p, mu = 1, sigma = 1, xi = 0.3, delta = 2)

rbgev(n, mu = 1, sigma = 1, xi = 0.3, delta = 2)

Arguments

x	Numeric vector of values for calculating density.
mu	location parameter
sigma	scale parameter (sigma > 0)
xi	shape parameter in R
delta	shape parameter (delta > -1)
q	Numeric vector of quantiles.
p	Numeric vector of probabilities.
n	Number of observations for random generation.

Details

This distribution was proposed by Cira EG Otiniano, Bianca S Paiva, Roberto Vila and Marcelo Bourguignon (2021)

Value

dbgev	density values
pbgev	distribution function values
qbgev	quantile function values
rbgev	random variates

Note

BGEV distribution is equivalent to the GEV distribution when delta = 0. When comparing BGEV with GEV from package EnvStats, the shape parameter of GEV is changed to -xi due to reparametrization

Author(s)

Thiago do Rego Sousa and Yasmin Lirio

References

Otiniano, Cira E. G., et al. (2023). A bimodal model for extremes data. Environmental and Ecological Statistics, 1–28. doi:10.1007/s10651-023-00566-7

Examples

par(mfrow = c(2, 2))
set.seed(1000)
r <- rbgev(n = 1000)
plot(r, type = "l", main = "BGEV Random Values")

hist(r, probability = TRUE, border = "white", ylim = c(0,1))
x <- seq(min(r), max(r), length = 201)
lines(x, dbgev(x), lwd = 2)

plot(sort(r), (1:1000)/1000, main = "Probability", ylab = "Probability")
lines(x, pbgev(x), lwd = 2)

round(qbgev(pbgev(q = seq(0, 3, by = 0.1)), 6),2)

n <- 10000
y <- rbgev(n)
qqplot(qbgev(ppoints(n)), y,
      main = "QQ-plot for bgev")
qqline(y, distribution = function(p) qbgev(p),
      probs = c(0.1, 0.6), col = 2)

---

Log-likelihood function for the BGEV distribution

Description

Log-likelihood function for the BGEV distribution

Usage

bgev_log_likelihood(x, pars)

Arguments

x	Numeric vector of observations.
pars	Vector of parameters (mu, sigma, xi, delta). See bgev.

Value

The log-likelihood value.

Author(s)

Thiago do Rego Sousa and Yasmin Lirio

---

Maximum Likelihood Estimation for the BGEV distribution

Description

Fits the BGEV distribution by maximum likelihood with a multistart local search (Nelder-Mead). The primary start comes from bgev_start_using_quantiles; additional starts are drawn inside a loose data-driven box (see bgev_start_box). The run with the highest log-likelihood is returned. Multistart guards against the multiple local maxima that are known to occur in bimodal likelihoods.

Usage

bgev_mle(
  x,
  likelihood = c("continuous_density", "grouped_likelihood"),
  h = 1,
  n_starts = 10,
  k = 3,
  control = list(maxit = 2000),
  ...
)

Arguments

x	Numeric vector of observations.
likelihood	Which likelihood to maximise: "continuous_density" (the default, using the density) or "grouped_likelihood" for discrete or rounded data, which uses the interval probability F(x + h/2) - F(x - h/2) of each observation (see h).
h	Rounding resolution for "grouped_likelihood" (default 1, for integer data). Ignored when likelihood = "continuous_density".
n_starts	Number of starting points for the multistart search (the quantile start plus n_starts - 1 perturbed starts).
k	Half-width multiplier for the data-driven box of perturbed starts.
control	List of control parameters passed to optim (defaults to a higher maxit than optim's own, which the bounded search otherwise hits on this likelihood).
...	Additional arguments passed to optim.

Details

The fit carries diagnostics (convergence, agreement across starts, and a support-boundary check) because the BGEV support depends on the parameters, so the usual regularity conditions can fail near the boundary. Standard errors from the inverse observed-information Hessian are returned, but only for an admissible (regular) optimum; near the boundary they are not reliable and are returned as NA.

Although the BGEV distribution is defined for delta > -1, estimation is restricted to delta > 0. This is deliberate: bimodality – the purpose of the model – occurs only for delta > 0, and for delta < 0 the density is unbounded at x = mu (the transform derivative (delta + 1)|x - mu|^delta diverges), giving a singular likelihood whose global maximum is a spurious spike on a data point. The revised BGEV reference estimates on delta >= 0 for the same reason.

The search is run on a reparametrised scale – log(sigma) and log(delta) – so that sigma > 0 and delta > 0 hold automatically. The only remaining penalty guards the data-dependent support (observations beyond the fitted endpoint), which is not a box constraint and cannot be transformed away. Among the multistart results, the fit returned is the highest-likelihood one whose Hessian at the optimum is positive-definite (a genuine interior maximum, rejecting spurious spikes); admissible is FALSE if none qualified. Estimates are reported on the natural scale.

Value

A list with par (named estimate c(mu, sigma, xi, delta)), se (standard errors from the inverse observed-information Hessian, NA unless admissible), loglik (maximised log-likelihood, positive, on the chosen scale), likelihood (which likelihood was used), convergence (optim code, 0 = success), start (the quantile start), n_starts, loglik_starts, agree (TRUE when several starts reach the selected maximum), admissible (TRUE when the returned optimum has a positive-definite Hessian), boundary (support-boundary diagnostic), and optimum (gradient norm, Hessian positive-definiteness, eigenvalue ratio and se).

Author(s)

Thiago do Rego Sousa and Yasmin Lirio

Examples

set.seed(1)
x <- rbgev(n = 200, mu = 1, sigma = 1, xi = 1, delta = 1)
fit <- bgev_mle(x)
fit$par

---

Profile log-likelihood for a BGEV parameter

Description

Computes the profile log-likelihood of one BGEV parameter over a grid, maximising over the remaining three at each grid point (Nelder-Mead). This is the diagnostic used in Otiniano et al. (2023) to check whether a fitted optimum is global and the parameter is well identified.

Usage

bgev_profile_likelihood(x, par, which, span = 0.5, n = 41, plot = TRUE)

Arguments

x	Numeric vector of observations.
par	Vector c(mu, sigma, xi, delta), e.g. bgev_mle(x)$par.
which	Index (1-4) or name of the parameter to profile.
span	Half-width of the grid, as a fraction of the parameter value.
n	Number of grid points.
plot	Logical; if TRUE, plot the profile curve.

Value

A data frame with the grid values of the profiled parameter and the profile log-likelihood.

Author(s)

Thiago do Rego Sousa

---

Compute the support of the BGEV distribution

Description

Returns the lower and upper limits of the support of the BGEV distribution

Usage

bgev_support(mu = 1, sigma = 1, xi = 0.3, delta = 2)

Arguments

mu	location parameter
sigma	scale parameter (sigma > 0)
xi	shape parameter in R
delta	shape parameter (delta > -1)

Details

It returns values with -Inf or Inf when the support is unbounded. When the shape parameter xi is different from zero, the support is truncated either at the left or at the right side of the real. Considering the support is particularly useful to estimating moments and to compute the likelihood function.

Value

A vector of length 2 with the lower and upper limits of the support

Author(s)

Thiago do Rego Sousa

---

Validate BGEV parameters

Description

Check if the provided parameters for the BGEV distribution are valid.

Usage

bgev_valid_params(mu, sigma, xi, delta)

Arguments

mu	location parameter, as in bgev
sigma	scale parameter (sigma > 0), as in bgev
xi	shape parameter in R, as in bgev
delta	shape parameter (delta > -1), as in bgev

Author(s)

Thiago do Rego Sousa

---