Fit a GEV CDN model via nonlinear optimization of the generalized maximum likelihood cost function.
1 2 3 4 5 6  gevcdn.fit(x, y, iter.max = 1000, n.hidden = 2,
Th = gevcdn.logistic, fixed = NULL,
init.range = c(0.25, 0.25),
scale.min = .Machine$double.eps, beta.p = 3.3,
beta.q = 2, sd.norm = Inf, n.trials = 5,
method = c("BFGS", "NelderMead"), max.fails = 100, ...)

x 
covariate matrix with number of rows equal to the number of samples and number of columns equal to the number of variables. 
y 
column matrix of target values with number of rows equal to the number of samples. 
iter.max 
maximum number of iterations of optimization algorithm. 
n.hidden 
number of hidden nodes in the GEV CDN model. 
Th 
hidden layer transfer function; defaults to 
fixed 
vector indicating GEV parameters to be held constant; elements chosen from 
init.range 
range for random weights on [ 
scale.min 
minimum allowable value for the GEV scale parameter. 
beta.p 

beta.q 

sd.norm 

n.trials 
number of repeated trials used to avoid shallow local minima during optimization 
method 
optimization method for 
max.fails 
maximum number of repeated exceptions allowed during optimization 
... 
additional arguments passed to the 
Fit a nonstationary GEV CDN model (Cannon, 2010) by minimizing a cost function
based on the generalized maximum likelihood of Martins and Stedinger (2000).
The hidden layer transfer function Th
should be set to
gevcdn.logistic
for a nonlinear model and to
gevcdn.identity
for a linear model. In the nonlinear case, the
number of hidden nodes n.hidden
controls the overall complexity of the
model. GEV parameters can be held constant (i.e., stationary) via the
fixed
argument. The form of the shifted beta distribution prior for the
GEV shape parameter is controlled by the beta.p
and beta.q
arguments. By default, these are set to values used in Cannon (2010).
Other alternatives include values recommended by Martins and Stedinger (2000)
(beta.p = 9
and beta.q = 6
) or values following from
the normal distribution reported by Papalexiou and Koutsoyiannis (2013)
(beta.p = 71.1
and beta.q = 44.7
). Weight penalty
regularization for the magnitude of the inputhidden layer weights can
be applied by setting sd.norm
to a value less than Inf
.
Values of the Akaike information criterion (AIC), Akaike information criterion
with small sample size correction (AICc), and Bayesian information criterion
(BIC) are calculated to assist in model selection. It is possible for such
criteria to fail in the face of overfitting, for example with a nonlinear model
and n.hidden
set too high, as the GEV distribution may converge on one or
more samples. This can usually be diagnosed by inspecting the scale parameter
for near zero values. In this case, one can apply a weight penalty
(via sd.norm
), although this rules out the use of AIC/AICc/BIC for model
selection as the effective number of model parameters will no longer equal the
number of weights in the GEV CDN model. Alternatively, a lower threshold (as a
proportion of the range of y
) for the scale parameter can be set
via scale.min
. Finally, bootstrap aggregation is available via
gevcdn.bag
as a second method for fitting GEV CDN models.
Note: values of x
and y
need not be standardized or rescaled
by the user. All variables are automatically scaled to range between 0 and 1
prior to fitting and parameters are automatically rescaled by
gevcdn.evaluate
.
a list consisting of
W1 
inputhidden layer weights 
W2 
hiddenoutput layer weights 
Attributes indicating the minimum/maximum values of x
and y
; the
values of Th
, fixed
, scale.min
; the negative of the
logarithm of the generalized maximum likelihood GML
, the negative of
the logarithm of the likelihood NLL
, the value of the penalty term
penalty
, the Bayesian information criterion BIC
, the Akaike
information criterion with (AICc
) and without (AIC
) small sample
size correction; and the number of model parameters k
are also
returned.
Cannon, A.J., 2010. A flexible nonlinear modelling framework for nonstationary generalized extreme value analysis in hydroclimatology. Hydrological Processes, 24: 673685. DOI: 10.1002/hyp.7506
Martins, E.S. and J.R. Stedinger, 2000. Generalized maximumlikelihood generalized extremevalue quantile estimators for hydrologic data. Water Resources Research, 36:737744. DOI: 10.1029/1999WR900330
Papalexiou, S.M. and Koutsoyiannis, D., 2013. Battle of extreme value distributions: A global survey on extreme daily rainfall. Water Resources Research, 49(1), 187201. DOI: 10.1029/2012WR012557
gevcdn.cost
, gevcdn.evaluate
,
gevcdn.bag
, gevcdn.bootstrap
,
dgev
, optim
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57  ## Generate synthetic data, quantiles
x < as.matrix(seq(0.1, 1, length = 50))
loc < x^2
scl < x/2
shp < seq(0.1, 0.3, length = length(x))
set.seed(100)
y < as.matrix(rgev(length(x), location = loc, scale = scl,
shape = shp))
q < sapply(c(0.1, 0.5, 0.9), qgev, location = loc, scale = scl,
shape = shp)
## Define a hierarchy of models of increasing complexity
models < list()
# Stationary model
models[[1]] < list(Th = gevcdn.identity,
fixed = c("location", "scale", "shape"))
# Linear model
models[[2]] < list(Th = gevcdn.identity)
# Nonlinear, 1 hidden node
models[[3]] < list(n.hidden = 1, Th = gevcdn.logistic)
# Nonlinear, 2 hidden nodes
models[[4]] < list(n.hidden = 2, Th = gevcdn.logistic)
## Fit models
weights.models < list()
for(i in seq_along(models)){
weights.models[[i]] < gevcdn.fit(x = x, y = y, n.trials = 1,
n.hidden = models[[i]]$n.hidden,
Th = models[[i]]$Th,
fixed = models[[i]]$fixed)
}
## Select model with minimum AICc
models.AICc < sapply(weights.models, attr, which = "AICc")
weights.best < weights.models[[which.min(models.AICc)]]
parms.best < gevcdn.evaluate(x, weights.best)
## 10th, 50th, and 90th percentiles
q.best < sapply(c(0.1, 0.5, 0.9), qgev,
location = parms.best[,"location"],
scale = parms.best[,"scale"],
shape = parms.best[,"shape"])
## Plot data and quantiles
matplot(x, cbind(y, q, q.best), type = c("b", rep("l", 6)),
lty = c(1, rep(c(1, 2, 1), 2)),
lwd = c(1, rep(c(3, 2, 3), 2)),
col = c("red", rep("orange", 3), rep("blue", 3)),
pch = 19, xlab = "x", ylab = "y", main = "gevcdn.fit")

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.