Blended GEV: a tutorial using R-INLA
Version 1.2
Daniela Castro-Camilo • daniela.castrocamilo@glasgow.ac.uk

A reparametrisation of the GEVThe blended Generalised extreme value (bGEV) model is an alternative to the usual GEV distribution when
the tail parameter ξ is positive, and it is designed to tackle the artificial boundary restrictions inherited by
the GEV. By artificial boundary restrictions, we mean the following: in practical applications, we assume that
the GEV is a reasonable approximation for the distribution of maxima over blocks, and we fit it accordingly.
This implies that GEV properties, such as finite lower endpoint in the case ξ > 0, are inherited by theoriginal maxima distribution, which might not be bounded-supported. This is particularly problematic in thepresence of covariates.
Before defining the bGEV distribution, we need to introduce an alternative parametrisation for the GEV.
The GEV distribution parametrised in terms of the location qα ∈ R, the spread sβ ∈ (0, ∞) and the tailparameter ξ ∈ R has the form( 
F (y | qα , sβ , ξ) = exp −y − qα+ `a,ξ
sβ (`1−β/2,ξ − `β/2,ξ )−1

−1/ξ ),+where a+ = max(a, 0) and for any a > 0, `a,ξ = {− log a}−ξ . Note that the case ξ = 0 simplifies to

 

y − qαF (y | qα , sβ ) = exp − exp −−`,αsβ (`1−β/2 − `β/2 )−1with `a = log{− log a}. There is a one-to-one mapping between (qα , sβ , ξ) and the usual location-scale-shape
GEV parameters, (µ, σ, ξ) (see, e.g., Coles, 2001, Chapter 3). For the case ξ 6= 0, the mapping is given byµ = qα −sβ (`α,ξ − 1),ξ(`1−β/2,ξ − `β/2,ξ )

σ=

sβ
,
(`1−β/2,ξ − `β/2,ξ )ξ = ξ.

The case ξ = 0 is interpreted as the limit when ξ → 0, i.e.,µ = qα +sβ `α,`β/2 − `1−/β/2

σ=

sβ
.
`β/2 − `1−/β/2This reparametrisation is proposed to provide a more meaningful interpretation of the parameters. In
statistics, the location-scale parametrisation is quite popular as it relates to the mean and the standard
deviation of the distribution. In skewed distributions such as the GEV, the mean is no longer a reasonableproxy for the location of the distribution. Moreover, the mean and variance of the GEV are only definedwhen ξ < 1 and ξ < 0.5, respectively. This effect of the tail parameter over the mean and variance does
not allow us to interpret the location-scale (and tail) GEV parametrisation as we do in other models. This
problem is particularly troublesome for the case where the parameters vary according to a set of covariates.
Assigning sensible priors to the GEV distribution with the usual parametrisation is also tricky when themean and variance are not defined.
1

The blended GEV model
The bGEV distribution is defined as
H(x | θ) = F(x | qαF , sβF , ξ)p(x|sp1 ,sp2 ,a,b) G(x | qαG , sβG )1−p(x|sp1 ,sp2 ,a,b) ,(1)

where θ = (qαF , sβF , ξ, qαG , sβG , sp1 , sp2 , a, b), F is the Frechét (or type II GEV) distribution with location
parameter qαF , spread parameter sβF , and tail parameter ξ, and G is the Gumbel (or type I GEV) distribution
with location parameter qαG and spread parameter sβG . The function p is a weight function defined as the
cumulative distribution function of a Beta distribution with shape parameters sp1 > 1 and sp2 > 1, evaluatedin the point (x − a)/(b − a), i.e.,
x−ap(x | sp1 , sp2 , a, b) = Pr Y ≤| sp1 , sp2 ,
b−a

(2)

where Y follows a Beta distribution with shape parameters sp1 > 1 and sp2 > 1. The Beta weight controls
the way the distributions F (Frechét) and G (Gumbel) influence the model H. The lower and upper boundsof the weight function (a and b, respectively) define the mixing area, i.e., where F and G are merged (seeFigure 1). Here, we choose them as quantiles of the Frechét distribution, i.e., a = F −1 (pa ) and b = F −1 (pb ),with 0 < pa , pb < 1. Below, pa and pb will be refered as the mixing area quantiles. The current INLA
implementation assumes that sp1 = sp2 = 5.

Figure 1: bGEV distribution (H, black) constructed from distributions F (Frechét, red), G (Gumbel, green)and Beta weight function p (purple). The shaded area is the mixing area, where F and G are merged.
In the following sections, we will learn how to fit the bGEV using R-INLA using three simulated examples
with increasing level of complexity.

Simulated example 1
To get familiar with the bGEV R-INLA implementation, we consider a simple model where the linear predictor
is linked to the α−quantile, qα . The model we want to fit isqα = η(x) = 1 + 0.4x

Data simulation
We start by generating n = 1000 samples from (3)2(3)

n = 1000
x = rnorm(n, sd=0.5) # we generate values for x from a N(0,0.5ˆ2) dist.
eta.x = 1 + 0.4*x
The spread and tail parameters are assumed to be covariate-free and unknown and are treated as hyperparameters within the INLA framework. We assume that the true spread and tail parameters are 0.3 and 0.1,
respectively.
spread = 0.3
tail = 0.1
To generate the GEV samples, we need to define the probabilities α and β that define the location (qα ) andspread (sβ ) parameters, respectively. In our case, they will be fixed to α = 0.5 (the median) and β = 0.25:
p.alpha = 0.5
p.beta = 0.25
Now we are ready to generate the samples. We use the function giveme.gev.par (See Section 5) to obtainthe usual GEV parameters (µ, σ, ξ) and plug them into the evd::rgev function to generate the samples
par = giveme.gev.par(q = eta.x, sbeta = spread, alpha = p.alpha, beta = p.beta,
xi = tail)y = numeric(n)for(i in 1:n)y[i] = rgev(1, loc = par$mu[i], scale = par$sigma, shape = par$xi)

Prior specification
The default prior distribution for the spread sβ is a Gamma with shape and rate parameters equal to 3 (notethat a log-scale is used below). For the tail parameter we consider a PC prior approach with parameters
λ = 7, low = 0 and high = 0.5. For more details on the PC prior for the tail GEV parameter see Section 4.2.
For the sake of illustration, we here define the priors for all the parameters involved. The priors for the
spread is
hyper.spread = list(initial = 1,
fixed=FALSE,
prior = "loggamma",
param = c(3, 3))
The prior for the tail parameter requires a bit more explanation. For computational reasons, R-INLA uses aninternal parametrisation of the tail parameter, which is unbounded. The hyperparameter specification isdefined for this internal parameter instead of the usual one, so in order to know how to specify a prior for the
tail parameter, we need to understand how both parametrisations are connected.
The map φ : [0, ∞) → R specifies the link between the usual and the internal parametrisations and it isdefined asexp(ξint ):= ξ,φ(ξint ) = low + (high − low)
1 + exp(ξint )where ξ ∈ [0, ∞) is the usual tail parameter and ξint refers to the unbounded internal tail parameter. Theinterval (low, high) constraints the possible values for ξ. The map and its inverse are defined in the function
map.tail (See Section 5).
In the internal parametrisation, the default initial value is -4 with (low, high) = (0, 0.5), which correspond toξ ≈ 0.04. If we want to provide an initial value of ξ = 0.1, then we can do
tail.interval = c(0, 0.5)
tail.intern = map.tail(tail, tail.interval, inverse=TRUE)3

We kept (low, high) = (0, 0.5) to ensure the existence of second moments1 .
The prior for the tail parameter is defined as
hyper.tail = list(initial = tail.intern,
prior = "pc.gevtail",param = c(7, tail.interval),
fixed= FALSE)
Note: if we have reasons to believe that ξ = 0 (Gumbel case), then a good initial value is ∞ (in the internalparametrisation). We can generalise the prior specification for the tail to allow ξ > 0 or ξ = 0 as follows
hyper.tail = list(initial = if (tail == 0.0) -Inf else tail.intern,
prior = "pc.gevtail",param = c(7, tail.interval),
fixed= if (tail == 0.0) TRUE else FALSE)
Therefore, the (default) hyperparameter specification for the bGEV model is
hyper.bgev = list(spread = hyper.spread,
tail = hyper.tail)Control variables
As part of the control.family argument in R-INLA, the argument control.bgev allows us to include
additional bGEV parameters. Specifically, the probabilities α and β, the mixing area quantiles pa and pb ,and the Beta weight function parameters, which we know are equal and fixed to 5 (for now, this cannot bechanged).
control.bgev = list(q.location = p.alpha,q.spread = p.beta,# quantile levels for the mixing partq.mix= c(0.05, 0.20),# the Beta(s1, s2) mixing distribution parameters.
# Hard-coded for the moment: s1=s2=5
beta.ab = 5)

INLA fit
The INLA formula for the bGEV model uses the function inla.mdata to allow the inclusion of simple linear
models in the spread and the tail parameters (for details on the inla.mdata function see Section 4.3). A nullmatrix can be used to indicate that these parameters are covariate-free (as it is the case in this example).
null.matrix = matrix(nrow = n, ncol= 0)
spread.x = null.matrix
tail.x = null.matrix
Note: matrices for the spread and tail covariates (here spread.x and tail.x) should always be defined andpassed to the R-INLA formula, even if they are empty.
Then, the INLA data and formula can be defined asdata.bgev = data.frame(y = y, intercept = 1, x = x, spread.x = spread.x, tail.x = tail.x)
formula = inla.mdata(y, spread.x, tail.x) ~ -1 + intercept + x
Note: data.bgev only contains three columns, as null matrices cannot be passed to data frames in R. As
spread.x and tail.x are not defined in data.bgev, INLA will search for these variables in the global R
1 Recall that the variance for the GEV is defined for ξ < 0.5.

4

environment. Therefore, spread.x and tail.x should be defined as null matrices in the global environment(as we do here). We write data.bgev this way for consistency with the following examples.
Finally, we fit the model
r1 = inla(formula,family = "bgev",
data = data.bgev,
control.family = list(hyper = hyper.bgev,
control.bgev = control.bgev),
control.predictor = list(compute = TRUE),
control.fixed = list(prec=1000),
control.compute = list(cpo = TRUE),
control.inla = list(int.strategy = "eb"),
verbose=FALSE, safe=TRUE)
A summary of the fitted fixed effects and hyperparameters can be obtained as follows
round(r1$summary.fixed,4)##mean
sd 0.025quant 0.5quant 0.975quantmode kld## intercept 0.9867 0.0031
0.98070.9867
0.9927 0.9867
0## x0.3798 0.0058
0.36840.37980.3911 0.3798
0round(r1$summary.hyperpar,4)##
mean
sd 0.025quant 0.5quant 0.975quant## spread for BGEV observations 0.2838 0.0078
0.26870.28370.2994## tail for BGEV observations
0.1071 0.0243
0.06390.1056
0.1590##mode## spread for BGEV observations 0.2835## tail for BGEV observations
0.1032

Simulated example 2
Although R-INLA does not allow more than one linear predictor, the bGEV implementation does allow for
simpler regression models on the spread and tail parameters. We then extend model (3) as follows

qα = η(x1 ) = 1 + 0.4x1
sβ = exp(0.1 + 0.3x2 )ξ = 0.1 + 0.2x3

Data simulation
As before, we start by generating n = 1000 samples from (4)n = 1000
x1 = rnorm(n)eta.x = 1 + 0.4*x1
The spread and tail parameter can be simulated as follows
x2 = rnorm(n, sd= 0.2)
s.x = exp(0.1 + 0.3*x2)
x3 = runif(n,-0.25,1)

5(4)

t.x = 0.1 + 0.2*x3
tail.intern = map.tail(t.x, tail.interval, inverse=TRUE) # internal xi
Note that since −0.5 < x3 < 2, we have that 0 < ξ < 0.5. As before, we use the function giveme.gev.par to
obtain the usual GEV parameters and plug them into the evd::rgev function to generate the samples (notethat α and β are the same as before)
par = giveme.gev.par(q = eta.x, sbeta = s.x, alpha = p.alpha, beta = p.beta,
xi = t.x)y = numeric(n)for(i in 1:n)y[i] = rgev(1, loc = par$mu[i], scale = par$sigma[i], shape = par$xi[i])Prior specification
Additional to what was discussed in Section 3.2, we can also adjust the priors for the regression coefficient of
the covariates for the spread and tail parameters, which we will call β1 and β2 . Within INLA, β1 and β2 aretreated as hyperparameters, with default prior given by a zero-mean Gaussian distribution with precisionequal to 300, as specified below.
hyper.beta1 = hyper.beta2 = list(prior = "normal",param = c(0, 300),
initial = 0)
In this case, the (default) hyperparameter specification for the bGEV model is
hyper.bgev = list(spread = hyper.spread,
tail = hyper.tail,beta1 = hyper.beta1,
beta2 = hyper.beta2)

INLA fit
As mentioned in Section 3.4, we can use the inla.mdata function to define linear models for the spread andthe tail parameters.
spread.x = x2
tail.x = x3
formula = inla.mdata(y, spread.x, tail.x) ~ -1 + intercept + x
Then, the INLA data can be defined asdata.bgev = data.frame(y = y, intercept = 1, x = x, spread.x = spread.x, tail.x = tail.x)
We fit the model using the same priors and mixing area quantiles (pa , pb ) as before.
r2 = inla(formula,family = "bgev",
data = data.bgev,
control.family = list(hyper = hyper.bgev,
control.bgev = control.bgev),
control.predictor = list(compute = TRUE),
control.fixed = list(prec=1000),
control.compute = list(cpo = TRUE),
control.inla = list(int.strategy = "eb"),
verbose=FALSE, safe=TRUE)
A summary of the fitted fixed effects and hyperparameters can be obtained as follows
6round(r2$summary.fixed,4)##mean
sd 0.025quant 0.5quant 0.975quantmode kld## intercept 0.7288 0.0174
0.69470.72880.7629 0.7288
0## x
0.0194 0.0254-0.0304
0.01940.0693 0.0194
0round(r2$summary.hyperpar,4)##
mean
sd 0.025quant 0.5quant## spread for BGEV observations
1.4583 0.0364
1.3884
1.4576## tail for BGEV observations
0.0554 0.0190
0.0236
0.0537## beta1 (spread) for BGEV observations 0.0453 0.0524-0.05800.0452## beta2 (tail) for BGEV observations
0.0073 0.0578-0.1067
0.0074##
0.975quantmode## spread for BGEV observations
1.5317 1.4560## tail for BGEV observations
0.0972 0.0502## beta1 (spread) for BGEV observations
0.1487 0.0451## beta2 (tail) for BGEV observations
0.1208 0.0078

Simulated example 3
The linear predictor can include more complicate structures defined as functions that depend on a set ofcovariates. By varying the form of these functions, we can accommodate a wide range of models, fromstandard and hierarchical regression to spatial and spatio-temporal models (Rue et al., 2009). Here we extend
the model in (4) by assuming the following structureqα = η(x1 , z1 , z2 ) = 1 + 0.4x1 + f1 (z1 ) + f2 (z2 )sβ = exp(0.1 + 0.3x2 + x4 )
ξ = 0.1 + 0.2x3(5)where f1 is a random walk of order 1 and f2 is an autoregressive process of order 2. Note that we also
extended the linear model for the spread parameter.

Data simulation
There are many ways we can simulate dara to fit (5). One alternative isn = 1000x = rnorm(n)z1 = seq(0, 6, length.out = n)z2 = 1:np = 2 # AR orderpacf = runif(p)phi = inla.ar.pacf2phi(pacf)
eta.x = 1 + 0.4*x + sin(z1) + c(scale(arima.sim(n, model = list(ar = phi))))The spread and tail parameter are simulated as before:x2 = rnorm(n, sd = 0.2)
x4 = rnorm(n, sd = 0.2)
s.x = exp(0.1 + 0.3*x2 + x4)
x3 = runif(n,-0.2, 0.2)
t.x = 0.1 + 0.2*x3
tail.intern = map.tail(t.x, tail.interval, inverse=TRUE) # internal xi

7

The data are simulated as
par = giveme.gev.par(q = eta.x, sbeta = s.x, alpha = p.alpha, beta = p.beta,
xi = t.x)y = numeric(n)for(i in 1:n)y[i] = rgev(1, loc = par$mu[i], scale = par$sigma[i], shape = par$xi[i])INLA fit
We have two covariates for the spread parameter and one for the tail, and we can have priors for each of theircoefficients. Below, hyper.beta1 and hyper.beta2 are the priors for the coefficients of the covariates in the
spread parameters, while hyper.beta3 is the prior for the coefficient of the covariate in the tail parameter.
hyper.beta1 = hyper.beta2 = hyper.beta3 = list(prior = "normal",param = c(0, 300),
initial = 0)
hyper.bgev = list(spread = hyper.spread,
tail = hyper.tail,beta1 = hyper.beta1,
beta2 = hyper.beta2,
beta3 = hyper.beta3)
Using the same mixing area quantiles (qa , qb ) as before, we only need to specify the INLA formula, the newdata, and run the model.
spread.x = x2
spread.xx = x4
tail.x = x3
# With this change of variable it is easier to keep track of the effect# of the covariates in each parameter, but it is not needed.
formula = inla.mdata(y, cbind(spread.x, spread.xx), tail.x) ~ -1 + intercept + x +f(z1, model = "rw1", scale.model=TRUE, constr=TRUE,
hyper = list(prec = list(prior = "pc.prec",param = c(0.1, 0.01)))) +f(z2, model = 'ar', order = 1,
hyper=list(prec=list(prior="pc.prec", constr=TRUE,
param=c(0.1,0.01)),
pacf1=list(param=c(0.5,0.8),
pacf2=list(param=c(0.5,0.8)))))
data.bgev = data.frame(y = y, intercept = 1, x = x, z1 = z1, z2 = z2,
spread.x = spread.x, spread.xx = spread.xx, tail.x = tail.x)
r3 = inla(formula,family = "bgev",
data = data.bgev,
control.family = list(hyper = hyper.bgev,
control.bgev = control.bgev),
control.predictor = list(compute = TRUE),
control.fixed = list(prec=1000,prec.intercept=1000),
control.compute = list(cpo = TRUE),
control.inla = list(int.strategy = "eb",cmin=0,b.strategy="keep"),
verbose=FALSE, safe=TRUE)
A summary of the fitted effects and hyperparameters can be obtained as follows
8round(r3$summary.fixed,4)##mean
sd 0.025quant 0.5quant 0.975quantmode kld## intercept 0.0154 0.0314-0.0460
0.01540.0769 0.0154
0## x0.2687 0.0163
0.23680.26870.3006 0.2687
0round(r3$summary.hyperpar,4)##
mean
sd 0.025quant 0.5quant## spread for BGEV observations
1.3613 0.0494
1.26831.3597## tail for BGEV observations
0.0048 0.0034
0.00060.0040## beta1 (spread) for BGEV observations
0.0312 0.0543-0.07430.0307## beta2 (spread) for BGEV observations
0.0774 0.0546-0.02990.0772## beta3 (tail) for BGEV observations
0.0031 0.0573-0.11170.0039## Precision for z1
264.7710 83.2254
132.1716 254.8788## Precision for z2
0.9897 0.1537
0.73730.9721## PACF1 for z2
0.9671 0.0063
0.9526
0.9679##0.975quantmode## spread for BGEV observations
1.46281.3555## tail for BGEV observations
0.01310.0017## beta1 (spread) for BGEV observations
0.13970.0285## beta2 (spread) for BGEV observations
0.18530.0767## beta3 (tail) for BGEV observations
0.11410.0067## Precision for z1
455.8932 235.7569## Precision for z2
1.33920.9311## PACF1 for z2
0.97740.9695

Notes
Linear predictor for the bGEV model
The linear predictor can take any combinations of the latent structures currently implemented. To see a listof the available latent models, we can dolibrary(INLA)
inla.list.models('latent')

PC prior for the tail parameter
Although non-informative priors are a common choice when little expert knowledge is available, the PC
prior approach allows us to select moderately informative prior distributions in a reasonable way. This
procedure penalises excessively complex models at a constant rate by putting an exponential prior on a
distance (specifically, the Kullback-Leibler distance, or KLD) to a simpler baseline model.
The PC prior for the tail GEV parameter is defined in terms of the KLD ξ 2 /(1 − ξ) for 0 ≤ ξ < 1. In R-INLA,
the prior is specified as
hyper = list(<theta> = list(prior="pc.gevtail", param=c(<lambda>, <low>, <high>)))where• <lambda>: the constant penalisation rate
• <low>,<high>: interval to restrict possible values for the tail parameter. The default is [0,0.5].
For more details on the PC prior for the tail GEV parameter see inla.doc("pc.gevtail").

9About inla.mdatainla.mdata(y, x1, x2,...) is a matrix where each row are replicates, and responses that are NAs are
ignored. The function inla.mdata accept covariates that are one or many vectors, matrices or data frames.
If we pass m covariates x1 , . . . xm , then each row (xi1 , xi2 , . . . , xim ) defines the covariates used for the i-th
row of y.

Tools
The function giveme.gev.par computes the usual GEV parameters (µ, σ, ξ) given the bGEV parameters(q, sβ , ξ). Note that in both parametrisations, the tail parameter is the same.
library(evd)giveme.gev.par = function(q, sbeta, alpha, beta, xi){.mu = function(q, sbeta, alpha, beta, xi) {a = -log(1-beta/2)
b = -log(beta/2)
c = -log(alpha)if (all(xi > 0.0)) {tmp0 = (cˆ(-xi) - 1)/xi
tmp1 = aˆ(-xi)
tmp2 = bˆ(-xi)
dbeta = (tmp1 - tmp2)/xi
return(q - (sbeta/dbeta) * tmp0)} else if (all(xi == 0.0)) {dbeta = log(b) - log(a)tmp0 = log(c)return(q + (sbeta/dbeta) * tmp0)} else {stop("mixed case not implemented")}}.sigma = function(q, sbeta, alpha, beta, xi) {a = -log(1-beta/2)
b = -log(beta/2)
if (all(xi > 0.0)) {tmp1 = aˆ(-xi)
tmp2 = bˆ(-xi)
dbeta = (tmp1 - tmp2)/xi
return(sbeta/dbeta)} else if (all(xi == 0.0)) {dbeta = log(b) - log(a)return(sbeta/dbeta)} else {stop("mixed case not implemented")}}
return(list(mu = .mu(q, sbeta, alpha, beta, xi),
sigma = .sigma(q, sbeta, alpha, beta, xi),
xi = xi))}

10The function map.tail specifies the link between the internal and usual parametrisations. In the code below,
interval constraints the possible values for the tail parameter, while inverse is a boolean variable indicating
whether φ or φ−1 should be computed.
map.tail = function(x, interval, inverse = FALSE) {if (!inverse) {return (interval[1] + (interval[2] - interval[1]) * exp(x)/(1.0 + exp(x)))} else {return (log((x-interval[1])/(interval[2]-x)))}}

References
1. Coles, S. (2001). An introduction to statistical modeling of extreme values (Vol. 208, p. 208). London:
Springer.
2. Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models
by using integrated nested Laplace approximations. Journal of the royal statistical society: Series B(statistical methodology), 71(2), 319-392.

11