# Install pacakge
install.packages("GLmom")

# Install from GitHub
remotes::install_github("sygstat/GLmom")

library(GLmom)

# Stationary GEV
data(streamflow)
result <- glme.gev(streamflow$r1)
result$para.glme  # GLME estimates: (mu, sigma, xi)

# Non-stationary GEV11
data(Trehafod)
result <- glme.gev11(Trehafod$r1)
result$para.glme  # (mu0, mu1, sigma0, sigma1, xi)

library(GLmom)
data(streamflow)
x <- streamflow$r1

# Default: beta penalty (adaptive)
result <- glme.gev(x)
result$para.glme  # GLME estimates
#> [1]  55.5346   8.7206  -0.4508
result$para.lme   # Traditional L-moment estimates
#>       xi    alpha    kappa
#>  55.4021   9.5229  -0.4075

# Compare different penalty functions
glme.gev(x, pen = "beta")$para.glme[3]    # xi = -0.4508
glme.gev(x, pen = "ms")$para.glme[3]      # xi = -0.2855 (Martins-Stedinger)
glme.gev(x, pen = "park")$para.glme[3]    # xi = -0.3474
glme.gev(x, pen = "no")$para.glme[3]      # xi = -0.4075 (no penalty = L-moment)

library(GLmom)
data(Trehafod)
x <- Trehafod$r1  # 53 years of river flow data

# Estimate with GLME (default: beta penalty)
result <- glme.gev11(x, ntry = 10)

# GLME estimates (with penalty)
result$para.glme
#>        mu0        mu1     sigma0     sigma1         xi
#>  79.312877   1.243031   2.788088   0.013611  -0.087837

# L-moment estimates (no penalty) - Shin et al. (2025b)
result$para.lme
#>        mu0        mu1     sigma0     sigma1         xi
#>  79.330671   1.243031   2.792161   0.013611  -0.076654

# Interpretation:
# - mu0 = 79.31: baseline location at t=0
# - mu1 = 1.24: location increases ~1.2 units per year
# - sigma0, sigma1: log-scale parameters
# - xi = -0.09: bounded upper tail (Weibull-type)

library(GLmom)
data(Trehafod)

# Simple interface - returns proposed L-moment estimates
result1 <- nsgev(Trehafod$r1, ntry = 10)
result1$para.prop
#>        mu0        mu1     sigma0     sigma1         xi
#>  79.330671   1.243031   2.792161   0.013611  -0.076654

# Comprehensive output - multiple estimation methods
result2 <- gado.prop_11(Trehafod$r1, ntry = 10)

# Compare methods:
result2$para.prop    # Proposed L-moment method
result2$para.gado    # GN16 method
result2$para.wls     # Weighted Least Squares
result2$lme.sta      # Stationary L-moments (mu, sigma, xi)

library(GLmom)
data(Trehafod)

# All penalty options for non-stationary model
penalties <- c("beta", "norm", "ms", "park", "cannon", "cd", "no")

results <- sapply(penalties, function(p) {
  r <- glme.gev11(Trehafod$r1, ntry = 10, pen = p)
  r$para.glme[5]  # shape parameter xi
})

print(round(results, 4))
#>  beta.xi norm.xi   ms.xi park.xi cannon.xi  cd.xi   no.xi
#>  -0.0878 -0.0904 -0.0895 -0.0692   -0.0901 -0.0667 -0.0767

library(GLmom)
data(streamflow)

# Beta penalty with custom hyperparameters
glme.gev(streamflow$r1, pen = "beta", p = 6, c1 = 20, c2 = 7)

# Or use preset choices (1-6 for beta, 1-4 for norm)
glme.gev(streamflow$r1, pen = "beta", pen.choice = 2)
glme.gev(streamflow$r1, pen = "norm", pen.choice = 1)

# Normal penalty with custom mean and std
glme.gev(streamflow$r1, pen = "norm", mu = -0.5, std = 0.2)

library(GLmom)
data(streamflow)
x <- streamflow$r1

# Model averaging with likelihood weights (default)
result <- ma.gev(x, quant = c(0.95, 0.99, 0.995), weight = 'like1', B = 200)

# Compare estimates
result$qua.mle    # MLE quantiles
#> [1]  107.28  147.02  167.35
result$qua.lme    # L-moment quantiles
#> [1]  110.43  184.35  234.27
result$qua.ma     # Model-averaged quantiles (recommended)
#> [1]  109.55  187.10  241.77

# Standard errors
result$fixw.se.ma # SE under fixed weights
result$ranw.se.ma # SE under random weights

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

# Using Bayesian Model Averaging
result3 <- ma.gev(x, quant = c(0.99), bma = TRUE, pen = "norm")
print(result3$qua.bma)

library(GLmom)
data(streamflow)
x <- streamflow$r1

# CD-MLE: Coles-Dixon penalized MLE
result_cd <- ma.gev(x, quant = c(0.99, 0.995), weight = "like1",
                    B = 50, CD = TRUE)
result_cd$mle.CD   # Coles-Dixon penalized MLE parameters (mu, sigma, xi)
result_cd$qua.CD   # Quantiles from CD-MLE

# REMLE: Restricted MLE with mean constraint
result_re <- ma.gev(x, quant = c(0.99, 0.995), weight = "like1",
                    B = 50, remle = TRUE)
result_re$remle1      # REMLE stage 1 parameters (mean constraint)
result_re$qua.remle1  # Quantiles from REMLE stage 1
result_re$remle2      # REMLE stage 2 parameters (mean + L-scale constraints)
result_re$qua.remle2  # Quantiles from REMLE stage 2

# Combined: all methods at once
result_all <- ma.gev(x, quant = c(0.98, 0.99, 0.995), weight = "like1",
                     B = 50, CD = TRUE, remle = TRUE)
result_all$qua.mle     # Standard MLE quantiles
result_all$qua.CD      # CD-penalized MLE quantiles
result_all$qua.remle1  # REMLE (mean) quantiles
result_all$qua.remle2  # REMLE (mean + L-scale) quantiles
result_all$qua.ma      # Model-averaged quantiles

# Load and explore datasets
data(Trehafod)
head(Trehafod)
#>   Year    r1
#> 1 1968 72.61
#> 2 1969 52.82
#> 3 1970 94.80