---
title: "Maximum Likelihood Estimation for the BGEV Distribution"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Maximum Likelihood Estimation for the BGEV Distribution}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.4)
library(bgev)
set.seed(1)
```

This vignette documents how `bgev_mle()` estimates the bimodal generalized
extreme value (BGEV) distribution: the likelihood and a singularity that forces
a restriction on `delta`, the starting-value strategy, a grouped likelihood for
discrete data, the diagnostics returned with every fit, and a Monte Carlo study
validating the estimator.

The parametrization follows the *revised* BGEV of Otiniano, Lisboa & Ribeiro
(2025), which adds a location parameter. With shape `xi`, location `mu`, scale
`sigma > 0` and bimodality parameter `delta > -1`, the CDF is
$F(x;\xi,\mu,\sigma,\delta) = F_{\mathrm{GEV}}(T(x);\xi,0,\sigma)$ with the
transform $T(x) = (x-\mu)\,|x-\mu|^{\delta}$ and derivative
$T'(x) = (\delta+1)\,|x-\mu|^{\delta}$.

## 1. The likelihood and its singularity at `x = mu`

The density is $f(x) = f_{\mathrm{GEV}}(T(x);\xi,0,\sigma)\,T'(x)$. The factor
$T'(x) = (\delta+1)|x-\mu|^{\delta}$ behaves very differently with the sign of
`delta`:

- for $\delta > 0$, $T'(x) \to 0$ as $x \to \mu$ (the bimodality *dip*);
- for $\delta < 0$, $T'(x) \to +\infty$ as $x \to \mu$ — **the density diverges
  at $x = \mu$.**

```{r singularity}
near_mu <- 10^-(1:4)
rbind(
  `delta<0` = sapply(near_mu, function(e) dbgev(2 + e, mu = 2, sigma = 1, xi = -0.3, delta = -0.3)),
  `delta>0` = sapply(near_mu, function(e) dbgev(2 + e, mu = 2, sigma = 1, xi = -0.3, delta =  0.5))
)
```

When `delta < 0` the likelihood is therefore **unbounded**: driving `mu` onto an
observation sends the log-likelihood to $+\infty$. This is the classic
unbounded-likelihood phenomenon — the same mechanism as the normal model with
$\sigma \to 0$ on a single point (Pawitan, 2001, §4.8). With continuous data it
has probability zero, but with tied or integer data the global maximizer parks
`mu` on a data atom, a spurious spike rather than a better fit.

**Consequence for estimation.** `bgev_mle()` restricts estimation to
`delta > 0`. This matches the revised reference, whose own estimation routine
uses `delta >= 0`, and it loses nothing of interest: bimodality occurs only for
`delta > 0` (`delta = 0` is the ordinary GEV). The *distribution* functions
(`dbgev()`, `pbgev()`, ...) still accept the full `delta > -1`; only the
estimator is restricted. Internally the search is reparametrized as
$(\mu, \log\sigma, \log\delta)$, so `sigma > 0` and `delta > 0` hold by
construction and no penalty cliff is needed for the parameter box.

## 2. Starting values: quantiles for the aim, a box for diversity

The optimizer is a **multistart local search** (Nelder–Mead, 1965). Two
complementary kinds of starting point feed it:

1. **One quantile-matching start** (`bgev_start_using_quantiles`), the *aimed*
   guess: it finds parameters whose model quantiles match the empirical
   quantiles at $p = (0.1, 0.3, 0.6, 0.9)$ by solving the four equations with a
   Newton-type solver.
2. **Several box starts** (`bgev_start_box`, `bgev_sample_starts`), the
   *diversity*: points drawn uniformly from a loose, data-driven box. They are
   deliberately scattered to probe for competing local maxima.

`bgev_mle()` runs a full local optimization from every start and keeps the best
*admissible* result (Section 4). The aimed start usually wins; the box starts
guard against multimodality.

This local-multistart design is **cheaper than a global optimizer** such as
differential evolution (`DEoptim`). The accepted pattern is to find the basin
with cheap starts and then polish locally (the global-to-local strategy;
Nocedal & Wright, 2006; Nash, 2014): the revised reference fits with plain
Nelder–Mead, and the Monte Carlo study below reaches the maximum in 93–100% of
replications, so a full global search is unnecessary for the regular regime. A
global fallback is reserved for cases where the starts disagree.

### Closed-form quantile estimators

The BGEV quantile function yields exact estimators at special probabilities.
Writing $Q(p)$ for the quantile function, at $p = e^{-1}$ we have
$-\log p = 1$, so the shape term $[(-\log p)^{-\xi}-1]/\xi$ vanishes and

$$Q(e^{-1}) = \mu \qquad \text{for all } \sigma, \xi, \delta.$$

```{r mu-closed-form}
x <- rbgev(2e5, mu = 5, sigma = 3, xi = -0.4, delta = 0.3)
c(closed_form = unname(quantile(x, exp(-1))), truth = 5)
```

So `quantile(x, exp(-1))` is an exact, parameter-free estimator of `mu`. For the
Gumbel-type case `xi = 0`, closed forms exist for the other two parameters as
well. Writing $\hat\mu = Q(e^{-1})$ and forming the **centered** quantiles
$q_k = Q(e^{-e^{k}}) - \hat\mu$ for $k = 1, 2$,

$$\delta = \frac{1}{\log_2(q_1/q_2)} - 1, \qquad \sigma = (-q_2)^{\delta+1}.$$

(Centering by $\hat\mu$ is essential: $q_k = -(\sigma k)^{1/(\delta+1)}$ only
after the location is removed.)

These are documented for completeness. They are **not** used as the estimator:
the general `xi` case has no such closed form, and experiments showed that
seeding `mu` with $Q(e^{-1})$ did not improve convergence or recovery over the
joint quantile solve — starts were not the bottleneck.

## 3. Discrete data: the grouped likelihood

BGEV is continuous, so a fit to rounded or integer data via the density is a
model mismatch (and, were `delta < 0` allowed, exactly where the atom spike
bites). The correct model records each value `x` to resolution `h` as the
interval probability

$$P\big(X \in [x - h/2,\ x + h/2]\big) = F(x + h/2) - F(x - h/2),$$

which is bounded by one and cannot blow up (Pawitan, 2001, §4.8). Select it with
the `likelihood` argument. To see that it does the right thing, we draw
continuous BGEV data, **record it to the nearest integer** (resolution
$h = 1$), and recover the generating parameters:

```{r grouped}
set.seed(1)
x_cont  <- rbgev(1000, mu = 0, sigma = 8, xi = 0.2, delta = 0.3)
x_obs   <- round(x_cont)                                  # rounded to nearest unit
fit_grp <- bgev_mle(x_obs, likelihood = "grouped_likelihood", h = 1)

rbind(truth    = c(mu = 0, sigma = 8, xi = 0.2, delta = 0.3),
      estimate = round(fit_grp$par, 2))
```

The grouped fit stays close to the truth despite the rounding. On continuous
data the grouped likelihood with a small `h` reproduces the continuous MLE; on
discrete data it is the appropriate choice and suppresses the tied-data warning
that the continuous density would raise.

## 4. Diagnostics returned with every fit

```{r fit-example}
x   <- rbgev(300, mu = 0, sigma = 1, xi = 0.2, delta = 1)
fit <- bgev_mle(x)
round(fit$par, 3)
round(fit$se, 3)
c(loglik = fit$loglik, convergence = fit$convergence,
  admissible = fit$admissible, agree = fit$agree)
```

- **`se`** — standard errors from the inverse observed-information Hessian
  (which is the Hessian of the negative log-likelihood at the estimate). They
  are returned only for an `admissible` optimum; near the parameter-dependent
  support boundary the regularity conditions fail, so `se` is `NA` there.
- **`convergence`** — the `optim` code (0 = success).
- **`agree`** — `TRUE` if several independent starts reached the same maximum
  (good evidence the optimum is global); `FALSE` if only one did (a weaker
  signal — the maximum may be local, worth inspecting a profile).
- **`admissible`** — whether the returned optimum is a *regular* interior
  maximum: a positive-definite Hessian with eigenvalue ratio above `1e-6`. The
  multistart keeps the best admissible optimum and rejects spurious / singular
  ones; `admissible = FALSE` warns that inference is not trustworthy there.
- **`boundary`** — proximity of the data to the parameter-dependent support
  endpoint, where the usual regularity conditions fail.
- **`optimum`** — gradient norm, Hessian positive-definiteness and eigenvalue
  ratio at the estimate.

`bgev_profile_likelihood()` complements these with a profile curve for any
parameter.

## 5. Monte Carlo validation

A study of `r "13,500"` fits (45 cells, 300 replications) simulated from known
parameters and refit, with `mu = 0`, `sigma = 1`,
`xi` in {−0.4, −0.2, 0, 0.2, 0.4}, `delta` in {0.25, 1, 3} and
`n` in {100, 250, 500}. The script is `benchmarks/mc_study.R`. Convergence was
93–100% in every cell (maximum failure rate 6%). The results below summarize the
two decision-relevant findings; they are static, taken from that run.

The estimator is **consistent and well-calibrated in the regular regime**
(`xi >= -0.2`): bias → 0, RMSE falls like `1/sqrt(n)`, and Wald coverage from the
observed-information Hessian is near the nominal 0.95.

```{r mc-coverage, echo = FALSE}
cov <- data.frame(
  xi  = c(-0.4, -0.2, 0.0, 0.2, 0.4),
  n100 = c(0.84, 0.98, 0.92, 0.91, 0.91),
  n250 = c(0.39, 0.97, 0.94, 0.94, 0.93),
  n500 = c(0.08, 0.97, 0.95, 0.93, 0.94)
)
knitr::kable(cov, col.names = c("xi", "n=100", "n=250", "n=500"),
             caption = "Wald 95% coverage of xi (delta = 1). Target 0.95.")
```

At the boundary (`xi = -0.4`) the model is **non-regular**: the
parameter-dependent support makes the MLE non-normal, so Wald coverage collapses
*and worsens with n* (0.84 → 0.39 → 0.08) and RMSE does not shrink. The package
flags this itself through the admissibility gate, which rejects most such fits:

```{r mc-admissible, echo = FALSE}
adm <- data.frame(
  xi = c(-0.4, -0.2, 0.0, 0.2, 0.4),
  d025 = c(0.46, 0.98, 1.00, 1.00, 1.00),
  d1   = c(0.12, 0.92, 1.00, 1.00, 1.00),
  d3   = c(0.10, 0.89, 1.00, 1.00, 1.00)
)
knitr::kable(adm, col.names = c("xi", "delta=0.25", "delta=1", "delta=3"),
             caption = "Admissible rate (positive-definite-Hessian gate), n = 500.")
```

Where Wald coverage holds the gate accepts ~100% of fits; where it fails it
rejects 88–90%. A user who checks `fit$admissible` is steered away from exactly
the fits whose Hessian-based standard errors cannot be trusted — and in that
regime `bgev_mle()` returns `se = NA` rather than an unreliable number.

## References

- Otiniano, C. E. G., Lisboa, M. N. S., & Ribeiro, T. K. A. (2025).
  *A Revised Bimodal Generalized Extreme Value Distribution: Theory and Climate
  Data Application.* Entropy, 27(7), 749.
  [doi:10.3390/e27070749](https://doi.org/10.3390/e27070749)
- Otiniano, C. E. G., et al. (2023). *A bimodal model for extremes data.*
  Environmental and Ecological Statistics.
  [doi:10.1007/s10651-023-00566-7](https://doi.org/10.1007/s10651-023-00566-7)
- Pawitan, Y. (2001). *In All Likelihood: Statistical Modelling and Inference
  Using Likelihood.* Oxford University Press. (§4.8, unbounded likelihood and the
  finite-precision / grouped likelihood.)
- Nelder, J. A., & Mead, R. (1965). *A simplex method for function minimization.*
  The Computer Journal, 7(4), 308–313.
- Nash, J. C. (2014). *Nonlinear Parameter Optimization Using R Tools.* Wiley.
- Nocedal, J., & Wright, S. J. (2006). *Numerical Optimization* (2nd ed.).
  Springer.