--- 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.