Package {GLmom}

Type:	Package
Title:	Generalized L-Moments Estimation for Extreme Value Distributions
Version:	2.0.0
Description:	Provides generalized L-moments estimation methods for the generalized extreme value ('GEV') distribution. Implements both stationary 'GEV' and non-stationary 'GEV11' models where location and scale parameters vary with time. Includes various penalty functions (Martins-Stedinger, Park, Cannon, Coles-Dixon) for shape parameter regularization. Also provides model averaging estimation ('ma.gev') that combines MLE and L-moment methods with multiple weighting schemes for robust high quantile estimation. The 'GLME' methodology is described in Shin et al. (2025a) <doi:10.48550/arXiv.2512.20385>. The non-stationary L-moment method is based on Shin et al. (2025b) <doi:10.1007/s42952-025-00325-3>. The model averaging method is described in Shin et al. (2026) <doi:10.1007/s00477-025-03167-x>. See also Hosking (1990) <doi:10.1111/j.2517-6161.1990.tb01775.x> for L-moments theory and Martins and Stedinger (2000) <doi:10.1029/1999WR900330> for penalized likelihood methods.
License:	GPL (≥ 3)
URL:	https://github.com/sygstat/GLmom
BugReports:	https://github.com/sygstat/GLmom/issues
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.3
Imports:	lmomco, nleqslv, stats, methods, robustbase, ismev, Rsolnp, zoo, graphics
Depends:	R (≥ 3.5.0)
NeedsCompilation:	no
Packaged:	2026-07-05 13:01:28 UTC; yg
Author:	Yonggwan Shin [aut, cre], Seokkap Ko [aut, ctb], Jihong Park [ctb], Yire Shin [aut, dtc], Jeong-Soo Park [aut, ths]
Maintainer:	Yonggwan Shin <syg.stat@etri.re.kr>
Repository:	CRAN
Date/Publication:	2026-07-05 15:50:17 UTC

Beta function integrand for penalty calculation

Description

Beta function integrand for penalty calculation

Usage

Befun(w, al = al, bl = bl, p = p, q = q)

Arguments

Value

Integrand value.

Methods for fitted GLmom objects

Description

The estimation functions return classed objects: glme.gev returns an object of class "glme", glme.gev11 of class "glme11", lme.gev11 of class "lme11", and ma.gev of class "magev". The objects remain ordinary lists, so all documented fields (e.g., $para.glme) stay accessible as before.

print shows the parameter estimates and the penalty settings; summary additionally reports the auxiliary estimates stored in the object; plot draws a quantile-quantile diagnostic: against the fitted GEV distribution for "glme", on the standard Gumbel scale (after the transformation of Eq. 8 in the manuscript) for "glme11" and "lme11", and via magev.qqplot for "magev".

Usage

## S3 method for class 'glme'
print(x, ...)

## S3 method for class 'glme'
summary(object, ...)

## S3 method for class 'glme'
plot(x, ...)

## S3 method for class 'glme11'
print(x, ...)

## S3 method for class 'glme11'
summary(object, ...)

## S3 method for class 'glme11'
plot(x, ...)

## S3 method for class 'lme11'
print(x, ...)

## S3 method for class 'lme11'
summary(object, ...)

## S3 method for class 'lme11'
plot(x, ...)

## S3 method for class 'magev'
print(x, ...)

## S3 method for class 'magev'
summary(object, ...)

## S3 method for class 'magev'
plot(x, ...)

Arguments

Value

x (or object), invisibly.

Examples

data(haenam)
fit <- glme.gev(haenam$X1)
fit
summary(fit)
plot(fit)

GN16 estimation for the non-stationary GEV11 model

Description

Estimates the parameters of the non-stationary GEV11 model by the quantile-based GN16 method (Gilleland and Katz, 2016; as adapted in Shin et al., 2025b), with both the original and the modified specification returned.

Usage

GN16.gev11(xdat, rob = FALSE)

Arguments

Value

A list containing:

• para.gado.org - Original GN16 estimates (mu0, mu1, sigma0, sigma1, xi).

• para.gado.mdfy - Modified GN16 estimates.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Gilleland, E. & Katz, R. W. (2016). extRemes 2.0: An extreme value analysis package in R. Journal of Statistical Software, 72(8), 1-39.

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

glme.gev11, lme.gev11, strup.gev11.

Examples

data(PhliuAgromet)
result <- GN16.gev11(PhliuAgromet$prec)
print(result$para.gado.org)

Martins-Stedinger prior function

Description

Computes the Martins-Stedinger Beta(6,9) prior probability for the GEV shape parameter on the interval [-0.5, 0.5].

Usage

MS_pk(para = para, p = 6, q = 9)

Arguments

Value

Prior probability value (scalar).

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Martins, E. S. & Stedinger, J. R. (2000). Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data. Water Resources Research, 36(3), 737-744. doi: 10.1029/1999WR900330

See Also

pk.beta for the adaptive beta penalty, glme.gev which uses these penalty functions.

Examples

# Evaluate MS prior at xi = -0.2
MS_pk(para = c(100, 20, -0.2))

Phliu Agrometeorological Station Data

Description

Annual maximum daily precipitation and station metadata from the Phliu Agrometeorological Station, Chanthaburi province, Thailand (1984-2023). A Mann-Kendall test indicates a significant increasing trend in the annual maxima (tau = 0.235, p = 0.033), so this dataset is used for demonstrating non-stationary (GEV11) estimation; it is the heavy rainfall application of Shin et al. (2025).

Usage

PhliuAgromet

Format

A data frame with 40 rows and 9 columns:

Station.ID

Station identifier (character)

year

Year of observation (numeric, 1984-2023)

prec

Annual maximum daily precipitation in mm (numeric)

Name

Station name (character)

zone

Climate zone code (character)

latitude

Station latitude in degrees (numeric)

longitude

Station longitude in degrees (numeric)

Starting.year

Record start year (integer)

Ending.year

Record end year (numeric)

Source

Phliu Agrometeorological Station, Thailand.

References

Shin, Y., Shin, Y., Park, J., & Park, J. S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

Examples

data(PhliuAgromet)
head(PhliuAgromet)


# Fit non-stationary GEV11 model by GLME
result <- glme.gev11(PhliuAgromet$prec, ntry = 5)
print(result$para.glme)

Prescott-Walden expected information matrix for GEV with fixed xi

Description

Internal function that computes the Prescott-Walden expected information (Hessian) matrix for GEV parameters (mu, sigma) with fixed xi.

Usage

PrescottW(par2 = NULL, xifix = NULL, nsam = NULL)

Arguments

Value

A 2x2 expected information matrix.

Asymptotic variance of model-averaged quantile estimates

Description

Internal function that computes asymptotic standard errors for model-averaged quantile estimates under both fixed-weight and random-weight assumptions.

Usage

asymp.var(mywt, covint, qqq = NULL, order = 2)

Arguments

Value

A list containing:

fin.se.MA.qua

SE under fixed weights for MA

fin.se.bma.qua

SE under fixed weights for BMA

adj.se.MA.qua

SE under random weights for MA

adj.se.bma.qua

SE under random weights for BMA

MatC

Cross-covariance array (numk x numk x numq)

Bangkok Maximum Rainfall Data

Description

Annual maximum daily rainfall data from Bangkok, Thailand. This dataset is used for demonstrating model averaging methods for high quantile estimation in extreme value analysis.

Usage

bangkok

Format

A data frame with 58 rows and 5 columns:

X1

Annual maximum daily rainfall in mm (numeric)

X2

2nd largest annual daily rainfall in mm (numeric)

X3

3rd largest annual daily rainfall in mm (numeric)

X4

4th largest annual daily rainfall in mm (numeric)

X5

5th largest annual daily rainfall in mm (numeric)

Source

Thai Meteorological Department (TMD; https://www.tmd.go.th)

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

Examples

data(bangkok)
head(bangkok)

# Estimate high quantiles using model averaging
result <- ma.gev(bangkok$X1, quant = c(0.99, 0.995))
print(result$qua.ma)

Bootstrap covariance of sample L-moments

Description

Computes the covariance matrix of the first three sample L-moments by nonparametric bootstrap. Used as a fallback when the direct covariance estimate is singular.

Usage

boot.cov(xdat, BB)

Arguments

Value

A 3x3 covariance matrix.

Select candidate shape parameter values for model averaging

Description

Internal function that selects K candidate shape parameter (xi) values from the profile likelihood confidence interval. Falls back to bootstrap LME quantiles if the profile likelihood fails.

Usage

cand.xi(
  data,
  hosking = NULL,
  mle = NULL,
  pick0 = 0.95,
  nint = 256,
  start = "mle",
  numk = NULL,
  figure = TRUE,
  cov.lme = NULL,
  bma = FALSE,
  pen = "beta"
)

Arguments

Value

A list containing:

kpar

Numeric vector of K candidate xi values

start

Starting method actually used

get.ci

Profile CI result (if start="mle")

ymin

Minimum y value for plotting (if start="mle")

Coles-Dixon Penalty Function

Description

Computes the Coles-Dixon prior probability for the shape parameter.

Usage

cd.hos(sxi)

Arguments

Value

Prior probability value.

Validate and resolve penalty settings

Description

Normalizes the penalty name (e.g., "Beta" to "beta") and, when pen.choice is given, resolves the preset hyperparameters: for pen="beta", choices 1-6 give (p, c1, c2) = (6,3,1), (6,5,2), (6,7,3), (2,3,0.5), (2,5,1), (2,7,1.5); for pen="norm", choices 1-4 give (mu, std) = (-0.5,0.25), (-0.5,0.15), (-0.6,0.25), (-0.6,0.15). These presets follow Table 1 and Eq. (15) of Shin et al. (2025).

Usage

check.penalty(
  pen = NULL,
  pen.choice = NULL,
  p = NULL,
  c1 = NULL,
  c2 = NULL,
  mu = NULL,
  std = NULL
)

Arguments

Value

A list with elements pen, p, c1, c2, mu, std.

Compute generalized L-moment distance probabilities

Description

Internal function that computes generalized L-moment distance and median-based distance probabilities for each candidate submodel.

Usage

com.prdist(
  data = NULL,
  numk = NULL,
  kfix = NULL,
  Vinv = NULL,
  detV = NULL,
  numom = NULL,
  hosking = NULL,
  trim = NULL,
  cov.type = "lambda"
)

Arguments

Value

A list containing:

prob.mtx

Matrix (numk x 2) of probabilities (col 1 = gLd, col 2 = med)

gdd

Matrix (numk x 2) of generalized distances

Compute profile likelihood confidence intervals

Description

Internal function that extracts confidence intervals from a profile log-likelihood curve at specified confidence levels.

Usage

comp.prof.ci(d, v, conf = NULL)

Arguments

Value

A list containing:

vmax

Maximum log-likelihood value

nllh

Same as vmax

xmax

xi value at maximum likelihood

ci1

Lower CI bounds for each confidence level

ci2

Upper CI bounds for each confidence level

ci_length

CI lengths for each confidence level

Construct covariance matrix C for model averaging SE

Description

Internal function that constructs the cross-covariance matrix between submodel quantile estimators using the delta method and quantile correlation approximation.

Usage

cons.MatC(mywt, cov22, quant)

Arguments

Value

A list containing:

MatC

Array (numk x numk x numq) of cross-covariance values

fin.se

SE under fixed weights for MA

fin.se.bma

SE under fixed weights for BMA

numk

Number of submodels

numq

Number of quantiles

wtgd

MA weights

bmaw

BMA weights

Dirichlet covariance matrix for weights

Description

Internal function that computes the covariance matrix of weights assuming a Dirichlet-type distribution.

Usage

cov.dir(wtgd)

Arguments

Value

A (numk x numk) covariance matrix.

Interpolate missing covariance matrices across submodels

Description

Internal function that fills in missing 2x2 covariance matrices for submodels using natural spline interpolation across the shape parameter.

Usage

cov.interp(numk, para3, cov2 = NULL)

Arguments

Value

Array (2 x 2 x numk) of covariance matrices with missing values filled by interpolation.

Delta method variance and cross-covariance for GEV quantiles

Description

Internal function that computes delta method variances for GEV quantile estimates. Can compute either the variance from a full 3-parameter MLE (d3yes=TRUE) or the cross-covariance between two submodels with fixed xi (d3yes=FALSE).

Usage

delta.gev(
  gevf3 = NULL,
  mle3i = NULL,
  cov2i = NULL,
  mle3j = NULL,
  cov2j = NULL,
  quant = NULL,
  d3yes = FALSE
)

Arguments

Value

A list containing:

v3

Variance for each quantile (if d3yes=TRUE)

covij

Cross-covariance for each quantile (if d3yes=FALSE)

Fit submodels and compute distance-based probabilities

Description

Internal function that fits GEV submodels with fixed xi for each candidate and computes generalized L-moment distance or median-based probabilities for weight construction.

Usage

dist.noboot(
  data = NULL,
  numk = NULL,
  hosking = NULL,
  boot.lme = TRUE,
  kpar = NULL,
  numom = NULL,
  ntry = 5,
  varcom = NULL,
  cov.lme = NULL,
  trim = NULL,
  cov.type = "lambda"
)

Arguments

Value

A list containing:

aic

AIC values for each submodel

mle3

Matrix (numk x 3) of submodel parameter estimates

kfix

List of submodel fitting results

cov2

List of 2x2 covariance matrices (if varcom=TRUE)

prob.mtx

Matrix (numk x 3) of distance-based probabilities

gdd

Matrix (numk x 2) of generalized distances

Empirical prior for BMA based on MLE and LME

Description

Internal function that computes an empirical prior distribution for the shape parameter based on the range between MLE and LME estimates.

Usage

emp.prior(mle3, mle, lme)

Arguments

Value

Numeric vector of prior probability values for each candidate.

Find the shape parameter maximizing the beta penalty

Description

Finds the value of xi (over a grid within the adaptive beta penalty support) which maximizes the beta penalty function. Used as one of the initial points for optimization.

Usage

find_max_beta.pk(lme.center, p = p, c0 = c0, c1 = c1, c2 = c2)

Arguments

Value

The xi value maximizing the penalty.

L-moment distance function for GEV11 model

Description

Computes the vector of differences between the population L-moments of the standard Gumbel distribution and the sample L-moments of the data transformed by trans.gum01 under the current GEV11 parameters. Solving this system equal to zero gives the L-moment estimates (Shin et al., 2025b).

Usage

fun.lme.gev11(a, xdat = xdat, pretheta = pretheta)

Arguments

Value

Numeric vector of length 3 of L-moment equation residuals.

Comprehensive Non-stationary GEV Estimation (deprecated)

Description

Estimates parameters of a non-stationary GEV distribution using multiple methods: Weighted Least Squares (WLS), GN16 method, and the proposed L-moment method from Shin et al. (2025, J. Korean Stat. Soc.).

As of GLmom v2.0.0 this function is deprecated: it was renamed to lme.gev11 (for the proposed method), and the auxiliary methods are available as strup.gev11 and GN16.gev11. This wrapper reassembles the v1.x output format from those functions and will be removed in a future release.

Usage

gado.prop_11(xdat, ntry = 20, ftol = 1e-06)

Arguments

Value

A list containing:

• para.prop - L-moment based estimates (proposed method).

• para.gado - GN16 method estimates.

• para.wls - Weighted least squares estimates.

• strup.org - Original non-stationary WLSE by Strup method.

• lme.sta - Stationary L-moment estimates.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

lme.gev11 for the recommended interface, strup.gev11, GN16.gev11, glme.gev11 for the full GLME method, nsgev for the simple interface.

Examples

data(PhliuAgromet)
# deprecated wrapper; suppressWarnings() silences the deprecation notice
result <- suppressWarnings(gado.prop_11(PhliuAgromet$prec, ntry = 5))
print(result$para.prop)
print(result$lme.sta)

MLE for GEV distribution using constrained optimization

Description

Internal function that computes maximum likelihood estimates of GEV parameters using the Rsolnp constrained optimizer with multiple random starting points.

Usage

gev.max(xdat, ntry = 5)

Arguments

Value

A list containing:

nsample

Sample size

conv

Convergence status (0 = success)

nllh

Negative log-likelihood at the optimum

mle

MLE estimates (mu, sigma, xi) in Hosking style

Modified profile likelihood for GEV shape parameter

Description

Internal function that computes the profile log-likelihood for the GEV shape parameter with linear extrapolation beyond the observed range. Used for constructing confidence intervals on xi.

Usage

gev.profxi.mdfy(
  data = NULL,
  mle = NULL,
  xlow,
  xup,
  pick.v = NULL,
  nint = 256,
  figure = FALSE
)

Arguments

Value

A list containing:

fail

Logical. TRUE if profile likelihood is degenerate

start

Starting method recommendation

w1

Profile CI results from comp.prof.ci()

ymin

Minimum y value for plotting

ymax

Maximum y value for plotting

MLE with return level and delta method SE

Description

Internal function that computes GEV MLE using both constrained optimization (solnp) and ismev::gev.fit, selects the better fit, and returns the covariance matrix for delta method SE.

Usage

gev.rl.delta(data, ntry = 5, quant)

Arguments

Value

A list containing:

nllh

Negative log-likelihood

mle

MLE estimates (mu, sigma, xi) in Hosking style

qua.mle

Quantile estimates at quant probabilities

data

Input data

cov

3x3 covariance matrix from MLE

quant

Input quantile probabilities

MLE for GEV with fixed shape parameter (single candidate)

Description

Internal function that computes MLE of GEV location and scale parameters with a fixed shape parameter using L-BFGS-B optimization. Optionally computes Prescott-Walden Hessian for variance estimation.

Usage

gev.xifix.sing(xdat = NULL, xifix = -0.1, ntry = 5, varcom = NULL)

Arguments

Value

A list of class "gev.xifix" containing:

conv

Convergence status (0 = success)

nllh

Negative log-likelihood at the optimum

mle

MLE estimates (mu, sigma, xi)

cov

2x2 covariance matrix (if varcom=TRUE)

GEV negative log-likelihood with fixed xi

Description

Internal function that computes the negative log-likelihood for the GEV distribution with a fixed shape parameter (Hosking parameterization).

Usage

gev.xilik(a, xifix = NULL, xdat = NULL)

Arguments

Value

Negative log-likelihood value (scalar). Returns 10^6 if invalid.

GEV negative log-likelihood with fixed xi (wrapper)

Description

Internal wrapper function that calls gev.xilik().

Usage

gev.xilik2(a, xifix = xifix, xdat = xdat)

Arguments

Value

Negative log-likelihood value (scalar).

Initialize parameters for GEV with fixed xi

Description

Internal function that generates initial (mu, sigma) values for optimization with fixed shape parameter, using Gumbel L-moment estimates and random perturbations.

Usage

ginit.xifix(data, ntry)

Arguments

Value

A matrix with ntry rows and 2 columns (mu, sigma).

Generalized L-moments estimation for generalized extreme value distribution

Description

This function estimates the parameters of the Generalized Extreme Value distribution by the generalized method of L-moment estimation (GLME), which minimizes the generalized L-moment distance plus a penalty (preference) function on the shape parameter.

Usage

glme.gev(
  xdat,
  ntry = 5,
  pen = "beta",
  pen.choice = 1,
  mu = -0.5,
  std = 0.2,
  p = 6,
  c1 = 3,
  c2 = 1,
  c0 = 0.35,
  q = NULL,
  show = FALSE,
  method = "BFGS",
  maxit = 70,
  abstol = 1e-05
)

Arguments

Details

The equations for the L-moments for LME of the GEVD are

\underline{\bf \lambda} - \underline{\bf l} = \underline{\bf 0},

where \underline{\bf \lambda} =(\lambda_1,\; \lambda_2,\; \lambda_3)^t and \underline{\bf l} =(l_1,\; l_2,\; l_3)^t. Next, we define the generalized L-moments distance (GLD) as;

(\underline{\bf \lambda} -\underline{\bf l})^t V^{-1} (\underline{\bf \lambda} -\underline{\bf l}),

where V is the variance-covariance matrix of the sample L-moments up to the third order.

Value

The glme.gev function returns a list containing the following elements:

• para.glme - The estimated parameters of the Generalized Extreme Value distribution.

• para.lme - The L-moment estimates of the parameters.

• nllh.glme - The penalized negative log-likelihood of the GLME solution.

• convergence - 0 if the optimization converged, 5 otherwise.

• pen, pen_pen.choice - The penalization method (and choice) used.

• p_q - (for beta penalty) The p and q values used.

• c1_c2, c0_c1_c2 - (for beta penalty) The hyperparameters used.

• mu_std - (for norm penalty) The mu and std values used.

• covinv.lmom - The inverse of the covariance matrix of the L-moments.

• lcovdet - The log determinant of the covariance matrix.

• data - The input data (used by the plot method).

The object has class "glme"; see GLmom-methods.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

See Also

glme.gev11 for non-stationary GEV estimation, ma.gev for model averaging estimation, glme.like for the objective function, quagev.NS for quantile computation.

Examples

# Load example heavy-rainfall data
data(haenam)
x <- haenam$X1

# Estimate GEV parameters using beta penalty (default)
result <- glme.gev(x, ntry = 5)
print(result$para.glme)

# Using Martins-Stedinger penalty
result_ms <- glme.gev(x, ntry = 5, pen = "ms")
print(result_ms$para.glme)

Generalized L-moments estimation for non-stationary GEV11 model

Description

This function estimates parameters of the non-stationary GEV11 model where mu(t) = mu0 + mu1*t and sigma(t) = exp(sigma0 + sigma1*t), by the generalized method of L-moment estimation (GLME).

Starting from the WLS pre-estimate of strup.gev11, the penalized objective (generalized L-moment distance based on the asymptotic covariance of the sample L-moments, plus penalty) is minimized over (mu0, sigma0, xi) with mu1 and sigma1 held at their pre-estimated values. With pen="no" the result equals the pure L-moment method lme.gev11 of Shin et al. (2025b).

Usage

glme.gev11(
  xdat,
  ntry = 5,
  ftol = 1e-06,
  opt.choose = "nllh",
  pen = "beta",
  pen.choice = 1,
  mu = -0.55,
  std = 0.2,
  p = 6,
  c1 = 5,
  c2 = 2,
  c0 = 0.35,
  q = NULL,
  show = FALSE,
  init.rob = TRUE,
  glme.pre = NULL
)

Arguments

Value

A list containing:

• para.glme - Proposed GLME estimates (mu0, mu1, sigma0, sigma1, xi).

• nllh.glme - Penalized negative log-likelihood of the solution.

• convergence - 0 if converged, 5 otherwise.

• pen, pen_pen.choice - Penalty method (and choice) used.

• p_q, c1_c2, c0_c1_c2 - (for beta) hyperparameters used.

• mu_std - (for norm) hyperparameters used.

• strup.org - WLSE by the Strupczewski method (strup.para).

• para.wls - Modified WLS estimates (strup.mdfy).

• para.gado - GN16 method estimates.

• lme.sta - Stationary L-moment estimates.

• para.lme - Pure L-moment estimates; returned only when pen="no" (then identical to para.glme). For other penalties use lme.gev11 directly.

• precis - (only when pen="no") precision of the L-moment equation solution.

• data - The input data (used by the plot method).

The object has class "glme11"; see GLmom-methods.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

glme.gev for stationary GEV estimation, lme.gev11 for the pure L-moment method (no penalty), strup.gev11 for the WLS pre-estimation, GN16.gev11 for the GN16 method, quagev.NS for non-stationary quantile computation.

Examples

data(PhliuAgromet)
x <- PhliuAgromet$prec


# Estimate non-stationary GEV11 parameters
result <- glme.gev11(x, ntry = 5)
print(result$para.glme)  # Proposed GLME estimates

Calculate the likelihood for Generalized L-moments estimation of GEV distribution

Description

This function calculates the likelihood (or more precisely, a penalized negative log-likelihood) for the Generalized L-moments estimation of the Generalized Extreme Value (GEV) distribution.

Usage

glme.like(
  par,
  xdat = xdat,
  slmgev = slmgev,
  covinv = covinv,
  lcovdet = lcovdet,
  mu = mu,
  std = std,
  lme = lme,
  pen = pen,
  p = p,
  c1 = c1,
  c2 = c2,
  c0 = 0.35,
  q = NULL
)

Arguments

Details

The function performs the following steps: 1. Checks if the parameters are within valid ranges. 2. Calculates the expected L-moments based on the current parameters. 3. Computes the difference between expected and sample L-moments. 4. Calculates the generalized L-moments distance. 5. Applies a penalization term based on the specified method. 6. Returns the sum of the L-moments distance and the penalization term.

Value

A numeric value representing the penalized negative log-likelihood. A lower value indicates a better fit.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

See Also

glme.gev which calls this function for optimization.

Examples

data(haenam)
x <- haenam$X1
slm <- lmomco::lmoms(x, nmom = 3)
cov_mat <- lmomco::lmoms.cov(x, nmom = 3)
lme_par <- lmomco::pargev(slm)$para
glme.like(par = lme_par, xdat = x, slmgev = slm,
          covinv = solve(cov_mat), lcovdet = log(det(cov_mat)),
          mu = -0.5, std = 0.2, lme = lme_par, pen = "beta",
          p = 6, c1 = 3, c2 = 1)

Goodness-of-fit by number of exceedances

Description

Internal function computing a goodness-of-fit measure comparing the expected and observed numbers of exceedances of the T-year levels for a set of return periods.

Usage

gof.ene_all(xdat, vecT = c(5, 10, 20, 40, 80), para = NULL, model = NULL)

Arguments

Value

Sum of the relative absolute differences.

Haenam Maximum Rainfall Data

Description

Annual maximum daily rainfall data from Haenam, South Korea. This dataset is used for demonstrating model averaging methods for high quantile estimation in extreme value analysis.

Usage

haenam

Format

A data frame with 52 rows and 2 columns:

year

Year of observation (integer, 1971-2022)

X1

Annual maximum daily rainfall in mm (numeric)

Source

Korea Meteorological Administration (KMA; https://www.weather.go.kr)

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

Examples

data(haenam)
head(haenam)

# Estimate high quantiles using model averaging
result <- ma.gev(haenam$X1, quant = c(0.98, 0.99, 0.995))
print(result$qua.ma)

Initialize parameters for GEV MLE estimation

Description

This function initializes parameters for GEV maximum likelihood estimation.

Usage

init.gevmax(data = NULL, ntry = NULL)

Arguments

Details

The function generates 'ntry' sets of initial parameters for the GEV distribution. It uses L-moment estimates as a starting point and then generates additional sets of parameters using random perturbations.

Value

A matrix with 'ntry' rows and 3 columns, where each row represents a set of initial parameters (location, scale, shape) for the GEV distribution.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

Initialize starting values for GLME optimization

Description

Generates multiple starting parameter sets for multi-start optimization in GLME estimation, for both the stationary GEV (model="gev00") and the non-stationary GEV11 (model="gev11"). Uses L-moment estimates as a base and adds random perturbations.

For model="gev11" the second column is returned on the log scale (matching sigma0 of the GEV11 parameterization), and when pretheta is supplied the first candidate is a perturbation of the pre-estimate.

Usage

init.glme(
  data = NULL,
  ntry = NULL,
  model = "gev00",
  pretheta = NULL,
  xdat = NULL
)

Arguments

Value

A matrix with ntry rows and 3 columns, where each row is a candidate starting point for optimization.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

See Also

glme.gev, glme.gev11 which use this function internally.

Examples

data(haenam)
inits <- init.glme(haenam$X1, ntry = 5)
print(inits)

Check if a number is odd

Description

Check if a number is odd

Usage

is.odd.me(x)

Arguments

Value

Logical. TRUE if x is odd, FALSE otherwise.

L-moment estimation with bootstrap standard errors

Description

Internal function that computes L-moment estimates of GEV parameters and optionally performs bootstrap resampling to obtain covariance matrices and standard errors for quantile estimates.

Usage

lme.boots(data, B = NULL, quant, boot = TRUE, trim = NULL)

Arguments

Value

A list containing:

lme

L-moment estimates (mu, sigma, xi)

qua.lme

Quantile estimates from LME

quant

Input quantile probabilities

cov.par

3x3 covariance of bootstrap parameter estimates (if boot=TRUE)

cov.lambda

3x3 covariance of bootstrap L-moments (if boot=TRUE)

qua.lme.se

SE of quantile estimates from bootstrap (if boot=TRUE)

L-moment parameter estimation for the non-stationary GEV11 model

Description

Estimates the parameters of the non-stationary GEV11 model (mu(t) = mu0 + mu1*t, sigma(t) = exp(sigma0 + sigma1*t), constant xi) by the pure L-moment method of Shin et al. (2025b, JKSS). Starting from the WLS pre-estimate of strup.gev11, the three L-moment equations (matching the sample L-moments of the Gumbel-transformed data with the population L-moments of the standard Gumbel distribution) are solved by Broyden's method for (mu0, sigma0, xi), and the best solution is chosen by a goodness-of-fit measure.

This function replaces gado.prop_11 of GLmom v1.x (renamed in v2.0.0, with slightly changed input/output). It is also equivalent to glme.gev11(xdat, pen="no").

Usage

lme.gev11(xdat, ntry = 5, ftol = 1e-05, show = FALSE)

Arguments

Value

A list containing:

• lme.gev11 - The L-moment estimates (mu0, mu1, sigma0, sigma1, xi). If no solution is found, the WLS pre-estimate is returned.

• precis - Precision (mean absolute residual of the L-moment equations) of the selected solution.

• data - The input data (used by the plot method).

The object has class "lme11"; see GLmom-methods.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

glme.gev11 for the penalized (GLME) version, strup.gev11 for the WLS pre-estimation, GN16.gev11 for the GN16 method, nsgev and gado.prop_11 for the v1.x compatibility wrappers.

Examples

data(PhliuAgromet)
result <- lme.gev11(PhliuAgromet$prec, ntry = 5)
print(result$lme.gev11)

Modified L-moments computation with failure handling

Description

A slightly modified version of lmomco::lmoms (by J. Park) which returns an ifail flag instead of stopping when the L-moments cannot be computed.

Usage

lmoms.md.park(x, nmom = 5, mtrim = FALSE, no.stop = FALSE)

Arguments

Value

A list of L-moments with an additional ifail element.

Model Averaging for GEV High Quantile Estimation

Description

This function estimates high quantiles of the Generalized Extreme Value (GEV) distribution using model averaging with mixed criteria. It combines Maximum Likelihood Estimation (MLE) and L-moment Estimation (LME) to construct candidate submodels and assign weights effectively.

Usage

ma.gev(
  data = NULL,
  quant = c(0.98, 0.99, 0.995),
  weight = "like1",
  numk = 12,
  B = 200,
  varcom = TRUE,
  trim = 0,
  fig = FALSE,
  bma = FALSE,
  pen = "norm",
  CD = FALSE,
  remle = FALSE
)

Arguments

Details

The model averaging approach works as follows:

1. MLE and LME of GEV parameters are computed

2. K candidate shape parameters (xi) are selected from profile likelihood CI

3. For each candidate xi, MLE with fixed xi is computed

4. Weights are assigned based on the selected method

5. Final quantile estimates are weighted averages across submodels

The weighting schemes include:

• 'like': AIC-based weights using likelihood with fixed xi

• 'gLd': Weights based on generalized L-moment distance

• 'med': Weights based on median and L-moment distance

• 'cvt': Conventional AIC weights

When bma=TRUE, Bayesian model averaging is applied with prior specified by pen.

Value

A list containing:

• mle.hosking - MLE estimates in Hosking style (mu, sigma, xi)

• qua.mle - Quantile estimates from MLE

• mle.cov3 - Covariance matrix of MLE (3x3)

• qua.se.mle.delta - Standard errors of MLE quantiles (delta method)

• lme - L-moment estimates (mu, sigma, xi)

• lme.cov3 - Covariance matrix of LME (bootstrap)

• qua.lme - Quantile estimates from LME

• qua.se.lme.boots - Standard errors of LME quantiles (bootstrap)

• qua.ma - Model-averaged quantile estimates

• w.ma - Weights used for model averaging

• fixw.se.ma - Asymptotic SE under fixed weights

• ranw.se.ma - Asymptotic SE under random weights

• surr - Surrogate model parameters (mu, sigma, xi)

• pick_xi - Selected xi values for K submodels

• qua.bma - (if bma=TRUE) BMA quantile estimates

• w.bma - (if bma=TRUE) BMA weights

• mle.CD - (if CD=TRUE) Coles-Dixon penalized MLE

• qua.CD - (if CD=TRUE) Quantile estimates from CD-MLE

• remle1 - (if remle=TRUE) Restricted MLE (first constraint)

• qua.remle1 - (if remle=TRUE) Quantile estimates from remle1

• remle2 - (if remle=TRUE) Restricted MLE (second constraint)

• qua.remle2 - (if remle=TRUE) Quantile estimates from remle2

• quant - The quantile probabilities used

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

See Also

glme.gev for stationary GLME estimation, magev.ksensplot for K sensitivity analysis, magev.qqplot for Q-Q diagnostic plots, magev.rlplot for return level plots.

Examples

# Load example data
data(haenam)
x <- haenam$X1

# Basic usage with likelihood weights
result <- ma.gev(x, quant = c(0.95, 0.99), weight = 'like1', B = 100)
print(result$qua.ma)     # Model-averaged quantiles
print(result$qua.mle)    # MLE quantiles for comparison
print(result$qua.lme)    # LME quantiles for comparison

# Using generalized L-moment distance weights
result2 <- ma.gev(x, quant = c(0.95, 0.99), weight = 'gLd', B = 100)
print(result2$w.ma)      # Model weights

K Sensitivity Plot for MAGEV

Description

Plots return level estimates, standard errors, and first-order differences across different numbers of candidate submodels K. This helps identify stable regions where estimates converge and select an optimal K value.

Usage

magev.ksensplot(
  data = NULL,
  q.cut = 0.6,
  mink = 4,
  maxk = 20,
  quant = c(0.99, 0.995)
)

Arguments

Details

The function computes MAGEV estimates for K ranging from mink to maxk. For each K, it calculates:

• Return level estimates (black points)

• Normalized standard errors (blue line)

• First-order differences (red line with triangles)

The optimal K is selected as the smallest value where both the normalized standard error and first-order difference are below their respective q.cut quantile cutoffs. The selected K is indicated by a purple vertical line.

Value

The optimal K value (integer) selected by the algorithm.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

See Also

ma.gev for the main model averaging function.

Examples

data(haenam)
optimal_k <- magev.ksensplot(haenam$X1)
print(optimal_k)

Q-Q Diagnostic Plot for MAGEV

Description

Creates a 2x2 panel of Q-Q plots comparing observed vs. fitted quantiles for different estimation methods: MLE, LME, surrogate (MA), and REMLE.

Usage

magev.qqplot(data = NULL, zx = NULL)

Arguments

Details

The function creates four Q-Q plots:

• Upper left: MLE (Maximum Likelihood Estimation)

• Upper right: LME (L-moment Estimation)

• Lower left: Surrogate MA (Model Averaging surrogate)

• Lower right: REMLE (Restricted MLE, if available)

Points close to the 45-degree diagonal line indicate good model fit.

Value

NULL. The function produces a plot as a side effect.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

See Also

ma.gev for the main model averaging function, magev.rlplot for return level plots.

Examples

data(haenam)
qq <- c(seq(0.01, 0.99, by = 0.01), 0.995, 0.999)
zx <- ma.gev(haenam$X1, quant = qq, weight = 'like1',
             numk = 9, varcom = FALSE, remle = TRUE)
magev.qqplot(data = haenam$X1, zx = zx)

Return Level Plot for MAGEV

Description

Displays fitted return levels with 95 return period on a log scale.

Usage

magev.rlplot(par = NULL, se.vec = NULL, data = NULL)

Arguments

Details

The plot shows:

• Black line: Fitted return level curve

• Blue lines: 95

• Black points: Observed data at empirical return periods

The x-axis (return period) is on a log scale, ranging from 0.1 to 900 years.

Value

NULL. The function produces a plot as a side effect.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., & Park, J. S. (2026). Model averaging with mixed criteria for estimating high quantiles of extreme values: Application to heavy rainfall. Stochastic Environmental Research and Risk Assessment, 40(2), 47. doi: 10.1007/s00477-025-03167-x

See Also

ma.gev for the main model averaging function, magev.qqplot for Q-Q diagnostic plots.

Examples

data(haenam)
ff <- c(seq(0.01, 0.09, by = 0.01), 0.1, 0.2, 0.3, 0.4, 0.5,
        0.6, 0.7, 0.8, 0.9, 0.93, 0.95, 0.98, 0.99,
        0.993, 0.995, 0.998, 0.999)
zx <- ma.gev(haenam$X1, quant = ff, weight = 'like1',
             numk = 9, varcom = TRUE)
magev.rlplot(par = zx$surr$par, se.vec = zx$ranw.se.ma, data = haenam$X1)

Construct pseudo-maxima for the GN16 method

Description

Internal step of GN16.gev11: builds the adjusted residual series (qmax) and estimates the log-scale trend coefficients.

Usage

make.qmax.gev11(xdat = NULL, orig.para = NULL, rob = FALSE)

Arguments

Value

A list with sig0, sig1, and qmax.

Coles-Dixon Penalized MLE for GEV

Description

Computes maximum penalized likelihood estimates for GEV parameters using the Coles and Dixon (1999) prior penalty on the shape parameter.

Usage

mle.gev.CD(xdat, ntry = 10)

Arguments

Value

A list containing:

References

Coles, S., & Dixon, M. (1999). Likelihood-based inference for extreme value models. Extremes, 2(1), 5-23.

Moving average smoother for quantiles and weights

Description

Internal function that applies moving average smoothing to quantile estimates and weights across candidate submodels.

Usage

movave(order, numk, numq, zp, wt = NULL)

Arguments

Value

A list containing:

ez

Matrix (numq x numk) of smoothed quantile estimates

ew

Numeric vector of smoothed weights

Adaptively expand or prune candidate xi set

Description

Internal function that adaptively adjusts the candidate xi set by removing low-weight candidates and adding new candidates near boundary or dominant regions to improve weight coverage.

Usage

new.kpar2(
  wtgd = NULL,
  numk = NULL,
  kpar = NULL,
  remove = 0.004,
  dist = 0.02,
  tre1 = 0.01,
  tre8 = 0.5
)

Arguments

Value

A list containing:

kpar2

Updated candidate xi values

aw

Adjustment indicator (0 = no change needed)

numk

Updated number of candidates

GLME objective function for GEV11 model (mu0, sigma0, xi optimization)

Description

GLME objective function for GEV11 model (mu0, sigma0, xi optimization)

Usage

nllh.glme.gev11(
  a,
  xdat = xdat,
  newtheta = newtheta,
  covinv = covinv,
  lcovdet = lcovdet,
  pen = pen,
  mu = mu,
  std = std,
  p = p,
  c0 = c0,
  c1 = c1,
  c2 = c2,
  q = q
)

Arguments

Value

Penalized negative log-likelihood value.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

Non-stationary GEV Parameter Estimation

Description

Estimates parameters of a non-stationary Generalized Extreme Value (GEV) distribution using the L-moment-based algorithm from Shin et al. (2025, J. Korean Stat. Soc.).

This is a convenience wrapper around lme.gev11, providing compatibility with the original nsgev package interface (and with GLmom v1.x).

Usage

nsgev(xdat, ntry = 20, ftol = 1e-06)

Arguments

Value

A list containing:

• para.prop - The proposed L-moment based estimates (mu0, mu1, sigma0, sigma1, xi).

• precis - Precision of the optimization.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

lme.gev11 for the recommended interface, glme.gev11 for the full GLME method with penalty functions, gado.prop_11 for detailed estimation results.

Examples

data(PhliuAgromet)
result <- nsgev(PhliuAgromet$prec, ntry = 5)
print(result$para.prop)

Multi-start penalized optimization for the GEV11 model

Description

Internal driver of glme.gev11. For pen="no" it delegates to lme.gev11. Otherwise it minimizes nllh.glme.gev11 from several starting points, using the fixed asymptotic covariance matrix of the first three sample L-moments (computed for n=300 and rescaled by 300/n; see Eq. (26) of Shin et al., 2025).

Usage

opt.glme.gev11(
  xdat,
  ntry = 5,
  ftol = 1e-05,
  show = FALSE,
  pretheta = NULL,
  model = "gev11",
  pen = "beta",
  mu = NULL,
  std = NULL,
  p = 2,
  q = NULL,
  c0 = 0.35,
  c1 = 3,
  c2 = 2,
  opt.choose = "nllh"
)

Arguments

Value

A list with para.glme, nllh.glme, convergence (and precis when pen="no").

GEV parameter estimation with fixed shape parameter

Description

Estimates GEV location and scale parameters from L-moments while keeping the shape parameter fixed at a user-specified value. Modified from lmomco::pargev().

Usage

pargev.kfix(lmom, kfix = 0.1, checklmom = TRUE, ...)

Arguments

Value

A list with components:

type

Character "gev"

para

Numeric vector of GEV parameters (xi=location, alpha=scale, kappa=shape)

source

Character "pargev"

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Hosking, J. R. M. (1990). L-moments: Analysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society, Series B, 52(1), 105-124. doi: 10.1111/j.2517-6161.1990.tb01775.x

See Also

glme.gev for GLME estimation, pargev for the original L-moment GEV fitting.

Examples

data(haenam)
lmom <- lmomco::lmoms(haenam$X1, nmom = 3)
pargev.kfix(lmom, kfix = -0.1)

GEV parameter estimation with fixed shape parameter (MAGEV internal)

Description

Internal function that estimates GEV location and scale parameters given a fixed shape parameter, using L-moment equations.

Usage

pargev.xifix(lmom, xifix = 0.1, checklmom = TRUE, ...)

Arguments

Value

A list with components:

type

Character "gev"

para

Numeric vector of GEV parameters (mu, sigma, xi)

source

Character "pargev"

Penalty dispatcher (negative log prior)

Description

Returns the negative log of the preference (prior) function for the selected penalty type. Dispatches to pk.norm.stnary ("norm"), the adaptive pk.beta ("beta"), the fixed literature beta penalties ("ms": Beta(6,9); "park": Beta(2.5,2.5); "cannon": Beta(2,3.3)), the Coles-Dixon exponential penalty ("cd"), or zero ("no").

Usage

penalty.fun(
  par,
  mu = NULL,
  std = NULL,
  lme = NULL,
  pen = "beta",
  p = NULL,
  c0 = 0.35,
  c1 = NULL,
  c2 = NULL,
  q = NULL
)

Arguments

Value

Negative log prior value (scalar).

Beta preference function for the GEV shape parameter

Description

Computes the data-adaptive Beta preference (penalty) function for the GEV shape parameter, used by both the stationary (glme.gev) and the non-stationary (glme.gev11) GLME methods (v2.0.0 unified implementation).

When q is NULL (adaptive mode), the support is (\max(-1, \hat\xi - c_0),\; \min(0.3, \hat\xi + c_0)) where \hat\xi is lme.center[3], and q = p + \min(|\hat\xi| c_1, c_2) for \hat\xi \le -0.05 (otherwise q = p). When q is given (fixed mode), the support is (min.xi, max.xi); this mode is used for the fixed literature penalties ("ms", "park", "cannon").

pk.beta.stnary is kept as a backward-compatible alias. Note that since v2.0.0 both names use this unified implementation (the v1.x default hyperparameters c0=0.3, c2=5 changed to c0=0.35, c2=2, and out-of-support values now return 1e-100 instead of 1).

Usage

pk.beta(
  para = NULL,
  lme.center = NULL,
  p = NULL,
  q = NULL,
  min.xi = -1,
  max.xi = 0.5,
  c0 = 0.35,
  c1 = 10,
  c2 = 2
)

pk.beta.stnary(
  para = NULL,
  lme.center = NULL,
  p = NULL,
  q = NULL,
  min.xi = -1,
  max.xi = 0.5,
  c0 = 0.35,
  c1 = 10,
  c2 = 2
)

Arguments

Value

A list containing:

pk.one

Preference function value (scalar)

p

Shape parameter p used

q

Shape parameter q used

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

See Also

pk.norm.stnary for the normal penalty, MS_pk for the Martins-Stedinger penalty, glme.gev, glme.gev11.

Examples

# Adaptive beta preference for xi = -0.2 centered at LME xi = -0.15
pk.beta(para = c(100, 20, -0.2),
        lme.center = c(100, 20, -0.15), p = 6, c1 = 3, c2 = 1)

Normal preference function for shape parameter

Description

Computes a normal distribution-based preference (penalty) function value for the GEV shape parameter. The alias new_pf_norm is provided for backward compatibility.

Usage

pk.norm.stnary(para = NULL, mu = NULL, std = NULL)

new_pf_norm(para = NULL, mu = NULL, std = NULL)

Arguments

Value

Preference function value (scalar).

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

See Also

pk.beta for the beta penalty, MS_pk for the Martins-Stedinger penalty, glme.gev which uses these penalty functions.

Examples

# Normal preference with mean=-0.5, sd=0.2 at xi=-0.2
pk.norm.stnary(para = c(100, 20, -0.2), mu = -0.5, std = 0.2)

Beta prior for model averaging

Description

Internal function that computes a Beta distribution prior probability for a candidate shape parameter in model averaging.

Usage

prior.beta.stnary.ma(x = NULL, xi_lme = NULL, p = 2, q = 5, al, bl)

Arguments

Value

Prior probability value (scalar).

Quantile function for GEV11 model

Description

Internal wrapper that computes time-varying quantiles for the GEV11 model given a return period.

Usage

qns.gev11(Tp = NULL, para = NULL, year = NULL)

Arguments

Value

Numeric vector of quantile values (one per time point).

Quantile function for non-stationary GEV (return period input)

Description

Internal function that computes time-varying quantiles for non-stationary GEV models given a return period.

Usage

qns.gev_all(Tp = NULL, para = NULL, year = NULL, model = NULL)

Arguments

Value

Numeric vector of quantile values (one per time point).

Quantile function for non-stationary GEV models

Description

Calculates quantiles for non-stationary GEV models including GEV11, GEV10, GEV20, and stationary GEV00.

Usage

quagev.NS(f = NULL, para = NULL, nsample = NULL, model = NULL)

Arguments

Value

A matrix of quantiles (nsample x length(f)) or a vector if f is scalar.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y., Park, J. & Park, J.-S. (2025). Generalized method of L-moment estimation for stationary and nonstationary extreme value models. arXiv preprint arXiv:2512.20385. doi: 10.48550/arXiv.2512.20385

See Also

glme.gev11 for non-stationary GEV estimation, glme.gev for stationary GEV estimation, ran.gev_all for random number generation.

Examples

# GEV11 model: time-varying quantiles
para <- c(84.55, 1.03, 2.91, 0.009, -0.08)  # mu0, mu1, sigma0, sigma1, xi
q99 <- quagev.NS(f = 0.99, para = para, nsample = 53, model = "gev11")
print(q99)

Random number generation for non-stationary GEV models

Description

Generates a random sample from a stationary or non-stationary GEV model. For non-stationary models the location and/or scale vary with time t = 1, ..., nsample: mu(t) = mu0 + mu1*t + mu2*t^2 and sigma(t) = exp(sigma0 + sigma1*t).

Usage

ran.gev_all(nsample, para, model = NULL)

Arguments

Value

A numeric vector of length nsample.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

quagev.NS for the corresponding quantile function, glme.gev11, lme.gev11.

Examples

set.seed(1)
# Non-stationary GEV11 sample of size 50
x <- ran.gev_all(50, para = c(100, 0.5, 3, 0.005, -0.2), model = "gev11")
plot(x, type = "l")

Revised Coles-Dixon prior function

Description

Internal function that computes a revised Coles-Dixon prior probability for the shape parameter, used as a penalty in BMA.

Usage

rcd(xi)

Arguments

Value

Numeric vector of prior probability values.

Restricted MLE for GEV (Mixed Estimation)

Description

Computes restricted maximum likelihood estimates for GEV parameters with constraints on the mean or median matching the sample statistics.

Usage

remle.gev(
  xdat,
  ntry = 5,
  rest = "mean",
  quant = c(0.99, 0.995),
  trim = 0,
  CD.mle = NULL,
  mle = NULL,
  second = TRUE,
  w.mpse = FALSE
)

Arguments

Value

A list containing:

Select parameters using goodness-of-fit criterion

Description

Internal function that selects the best parameter set among candidate solutions using a goodness-of-fit measure based on the expected vs. observed number of exceedances (see Shin et al., 2025b, JKSS).

Usage

sel.para_all(xdat, para.sel = NULL, model = NULL, obj.fun = NULL)

Arguments

Value

A list containing:

para

Selected parameter vector

min.itry

Index of selected try

gof

Goodness-of-fit values

obj.fun

Objective function values

Set parameters based on non-stationary model type

Description

Internal function that maps a flat parameter vector to named location (mu), scale (sig), and shape (xi) components depending on the model specification (GEV00, GEV10, GEV11, GEV20, GEV01).

Usage

set.para.model(para, model = NULL)

Arguments

Value

A list containing:

mu

Numeric vector c(mu0, mu1, mu2) for location

sig

Numeric vector c(sigma0, sigma1) for log-scale

xi

Shape parameter (scalar)

Set BMA prior distribution for candidate shape parameters

Description

Internal function that sets up the prior distribution for Bayesian Model Averaging (BMA) weights. Supports both normal and beta priors, with hyperparameters adapted based on the weighting method and L-moment estimate of xi.

Usage

set.prior(pen = NULL, numk = NULL, xi_lme = NULL, kpar = NULL, weight = NULL)

Arguments

Value

A list containing:

prior

Numeric vector of prior probabilities (length numk)

prior_mu_std

(if pen="norm") Mean and std of the normal prior

p_q_beta

(if pen="beta") p and q parameters of the beta prior

Weighted least squares estimation for the non-stationary GEV11 model

Description

Estimates the parameters of the non-stationary GEV11 model by the weighted least squares (WLS) method of Strupczewski and Kaczmarek (2001), with the final-specification modification described in Shin et al. (2025b, JKSS). The result is used as the pre-estimate (starting value) by glme.gev11 and lme.gev11.

Usage

strup.gev11(xdat, init.rob = TRUE, wls.rob = FALSE, strup.fin.rob = TRUE)

Arguments

Value

A list containing:

• strup.para - Original WLS estimates (mu0, mu1, sigma0, sigma1, xi).

• strup.mdfy - Modified (re-specified) WLS estimates.

Author(s)

Yonggwan Shin, Seokkap Ko, Jihong Park, Yire Shin, Jeong-Soo Park

References

Strupczewski, W. G. & Kaczmarek, Z. (2001). Non-stationary approach to at-site flood frequency modelling II. Weighted least squares estimation. Journal of Hydrology, 248(1-4), 143-151.

Shin, Y., Shin, Y. & Park, J.-S. (2025). Building nonstationary extreme value model using L-moments. Journal of the Korean Statistical Society, 54, 947-970. doi: 10.1007/s42952-025-00325-3

See Also

glme.gev11, lme.gev11, GN16.gev11.

Examples

data(PhliuAgromet)
result <- strup.gev11(PhliuAgromet$prec)
print(result$strup.mdfy)

Find surrogate GEV parameters for model-averaged quantiles

Description

Internal function that fits a single GEV distribution to the model-averaged quantile curve, producing surrogate GEV parameters that approximate the model-averaged quantile function.

Usage

surrogate(zpf.surr = NULL, xqa, init.surr)

Arguments

Value

A list (from optim()) with additional components:

par

Surrogate GEV parameters (mu, sigma, xi)

zp.surrmodel

Quantiles from the surrogate model at xqa

zp.MA

The input model-averaged quantiles

Time-varying moment step of the GN16 method

Description

Internal step of GN16.gev11: fits a stationary GEV to the pseudo-maxima and converts the result to the GEV11 parameterization (both original and modified specification).

Usage

## S3 method for class 'm.gev11'
time(qmax = NULL, orig.para = NULL, rob = FALSE)

Arguments

Value

A list with para.gado.org and para.gado.mdfy.

Transform non-stationary GEV data to standard Gumbel

Description

Transforms observations from a non-stationary GEV11 model with parameters para = (mu0, mu1, sigma0, sigma1, xi) to the standard Gumbel scale.

Usage

trans.gum01(xdat, para = NULL, na.ok = FALSE)

Arguments

Value

Numeric vector on the standard Gumbel scale.

Compute model averaging weights

Description

Internal function that computes weights for model averaging across K candidate GEV submodels. Supports multiple weighting schemes including likelihood-based, generalized L-moment distance, and median-based methods. Optionally includes BMA prior integration.

Usage

weight.com(
  data = NULL,
  numk = NULL,
  hosking = NULL,
  kpar = NULL,
  numom = 3,
  xqa = NULL,
  varcom = TRUE,
  boot.lme = TRUE,
  cov.lme = NULL,
  surr = FALSE,
  type = "full",
  trim = NULL,
  cov.type = "ratio",
  bma = TRUE,
  pen = "norm"
)

Arguments

Value

A list containing:

weight

Weighting method name

numk

Number of submodels

wtgd

MA weight vector (length numk)

bmaw

BMA weight vector (length numk)

zp

Quantile matrix (numq x numk)

kpar

Candidate xi values

prob.call

Submodel fitting results

xzp

Surrogate quantile matrix (if surr=TRUE)

Likelihood-based weights with fixed xi

Description

Internal function that computes AIC-based model averaging weights using the profile likelihood with L-moment submodel estimates for each candidate xi.

Usage

wlik.xifix(
  data = NULL,
  numk = NULL,
  kpar = NULL,
  weight = NULL,
  pertr = 0.85,
  varcom = FALSE,
  type = NULL,
  prior = NULL,
  trim = NULL
)

Arguments

Value

A list containing:

bmaw

BMA weight vector

wt

MA weight vector

aic

AIC values for each submodel

kfix

Matrix (3 x numk) of submodel parameters

cov2

List of 2x2 covariance matrices (if varcom=TRUE)

Log-scale regression step of the Strupczewski WLS method

Description

Internal step of strup.gev11: regresses the log absolute residuals on time to estimate the log-scale coefficients (sigma0, sigma1).

Usage

wls.gev11(xdat, res = NULL, rob = NULL)

Arguments

Value

A list with element sig (the two coefficients).