# Estimating parameters from quantiles, GEV distribution

If we know the quantiles of the GEV distribution how can we determine the parameters?  John Cook explores this problem in relation to eliciting prior probabilities from subject matter experts.  If subject matter experts can provide quantile estimates then these can be used to construct prior probability distributions for use in Bayesian analysis.

Here I have a hydrological application in mind.  When using design rainfalls we generally have a published set of  rainfall depths (or intensities) for a corresponding set of annual exceedance probabilities, that is, we have a set of quantiles.  Often the GEV distribution is used to generate these quantiles but the parameters are not provided.  If the parameters can be estimated then it would be possible to calculate the design rainfall depths for any exceedance probability not just those that are tabulated.

The quantile function for the GEV distribution is:

$Q(p) = \xi + \frac{\alpha}{\kappa} \left[ (-\log(1-p))^{-\kappa} - 1\right] \,\,\, \kappa \ne 0$

$Q(p) = \xi - \alpha (\log (-\log (1 - p)))\,\,\, \kappa = 0$      (1)

Where $Q$ is the quantile, $p$ is an exceedance probability, $\xi$ is the location parameter, $\alpha$ is the scale parameter and $\kappa$ is the shape parameter.  Actually there are several slightly different ways the GEV distribution is parameterised but I’ll adopt the equations above.

We have three parameters so we need three quantiles to undertake the estimation.

Assume three quantiles $Q_1, Q_2, Q_3$ at exceedance probabilities $p_1, p_2, p_3$ and that $k \ne 0$

$Q_1(p_1) = \xi + \frac{\alpha}{\kappa} \left[ (-\log(1-p_1))^{-\kappa} - 1\right]$       (2)

$Q_2(p_2) = \xi + \frac{\alpha}{\kappa} \left[ (-\log(1-p_2))^{-\kappa} - 1\right]$       (3)

$Q_3(p_3) = \xi + \frac{\alpha}{\kappa} \left[ (-\log(1-p_3))^{-\kappa} - 1\right]$       (4)

To eliminate $\xi$ subtract (3) from (2).

$Q_1(p) - Q_2(p_2) = \frac{\alpha}{\kappa} \left[ (-\log(1-p_1))^{-\kappa} - (-\log(1-p_2))^{-\kappa} \right]$       (5)

Similarly, remove $\alpha$ by division to achieve an equation in $\kappa$ alone.

$\frac{Q(p_1) - Q(p_2)}{Q(p_1) - Q(p_3)} = \frac {\left[ (-\log(1-p_1))^{-\kappa} - (-\log(1-p_2))^{-\kappa} \right] }{\left[ (-\log(1-p_1))^{-\kappa} - (-\log(1-p_3))^{-\kappa} \right] }$       (5)

Equation (5) can be solved numerically to obtain $\kappa$ then substitute back into (4) to obtain $\alpha$ and then into (1) to obtain $\xi$

I’ll use a root finding approach to determine $\kappa$ from equation (5), which can be rearranged to give.

$f(\kappa)= \frac {\left[ (-\log(1-p_1))^{-\kappa} - (-\log(1-p_2))^{-\kappa} \right] }{\left[ (-\log(1-p_1))^{-\kappa} - (-\log(1-p_3))^{-\kappa} \right] } - \frac{Q(p_1) - Q(p_2)}{Q(p_1) - Q(p_3)}$       (6)

Note that this has a singularity at $k = 0$.  We’ll address this by interpolating equation (6) through zero based on positive and negative $\kappa$ values close to zero.

An example:

Consider a subset of three quantiles for  24 hour rainfall depths.

AEP Quantile (mm)
0.5 71.4
0.1 110.5
0.01 169.4

Plotting equation (6) for these values suggests that a root finding approach will work ok

FindKappa <- function(k, Q1, Q2, Q3, p1, p2, p3){

F <- function(k, Q1, Q2, Q3, p1, p2, p3){
((-log(1-p1))^-k - (-log(1-p2))^-k)/((-log(1-p1))^-k - (-log(1-p3))^-k) - (Q1 - Q2)/(Q1 - Q3)
}

# interpolate through zero

k1 <- -1e-7
k2 <- 1e-7
y1 <- F(k1, Q1, Q2, Q3, p1, p2, p3)
y2 <- F(k2, Q1, Q2, Q3, p1, p2, p3)

F_nearZero <- approxfun(c(k1, k2), c(y1, y2))

if(k > k1 & k < k2) {
F_nearZero(k)
} else {
F(k, Q1, Q2, Q3, p1, p2, p3)
}
}

# plot graph

Q1 <- 71.4
Q2 <- 110.5
Q3 <- 169.4

p1 <- 0.5
p2 <- 0.1
p3 <- 0.01

k.seq <- seq(-1, 1, length.out = 99)
par(oma = c(0,2,1,1))
Y <- sapply(k.seq, FindKappa, Q1, Q2, Q3, p1, p2, p3)
plot(k.seq, Y,
type = 'l',
xlab = expression(kappa),
ylab = expression(paste('f(', kappa, ')')))
abline(h=0, lty = 2)



The uniroot function will find the required $k$ value

uniroot(FindKappa, c(-1, 1), Q1, Q2, Q3, p1, p2, p3 )

$root [1] 0.08884866$f.root
[1] -2.131346e-06



The required value of $\kappa$ is 0.0888.

Functions to return the scale and location are as follows.

CalcScale <- function(Q1, Q2, p1, p2, k){
k*(Q1-Q2)/((-log(1-p1))^-k - (-log(1-p2))^-k)
}

CalcLocation <- function(Q1, p1, scale, k){
Q1 - (scale/k) * ((-log(1 - p1))^-k - 1)
}



Try them with the example data.

Q1 <- 71.4
Q2 <- 110.5
Q3 <- 169.4

p1 <- 0.5
p2 <- 0.1
p3 <- 0.01

my.k <- uniroot(FindKappa, c(-1, 1), Q1, Q2, Q3, p1, p2, p3 )\$root
my.k
#[1] 0.08884866

my.scale <- CalcScale(Q1, Q2, p1, p2, k = my.k)
my.scale
#[1] 18.45587

my.location <- CalcLocation(Q1, p1, scale = my.scale, k = my.k)
my.location
#[1] 64.52434

We now have the three parameters:

• shape ($\kappa$) = 0.0885
• scale = 18.456
• location = 64.524

Lets use these to calculate the quantile for the 2% AEP.

First, we need a quantile function based on equation 1 above then input the required parameters.  The quantile function needs to take account of the change in formula when $\kappa$  = 0.  This can cause numerical issues as $\kappa$ approaches 0 which can be overcome by interpolating between points on either side of zero.  A similar approach was used in the FindKappa function above.  For further details, see this Stackoverflow discussion.

Qgev <- function(p, location, scale, k){
# p is an exceedance probability

.F <- function(p, location, scale, k) {
(location + (scale/k)*((-log(1-p))^-k - 1))
}

k1 <- -1e-7
k2 <- 1e-7
y1 <- .F(p, location, scale, k1)
y2 <- .F(p, location, scale, k2)

F_nearZero <- approxfun(c(k1, k2), c(y1, y2))

if(k > k1 & k < k2) {
return(F_nearZero(k))
} else {
return(.F(p, location, scale, k))
}
}

Qgev(0.02, 64.52434, 18.45587, 0.08884866)
# [1] 150.5977



The predicted 24 hour, 2% rainfall depth is 150.6 mm.

Probability of a given amount of rainfall

So far we’ve been working to predict the rainfall depth for a given Annual Exceedance Probability.  It is also possible to calculate the Annual Exceedance Probability given a rainfall depth.

Example:

On a wet day there was 140 mm in 24 hours.  What was the AEP of this event?

Here we need the distribution function for the GEV which is available in the evd package.  Using the parameters calculated before:

library(evd)
evd::pgev(q = 140, loc = 64.50045, scale = 18.47377, shape = 0.08856304, lower.tail = FALSE )
# [1] 0.0300979



Therefore, 140 mm of rain in 24 hours has an AEP of 3%.