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

Kappa-fn

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%.

2 thoughts on “Estimating parameters from quantiles, GEV distribution

  1. Pingback: Interpolating design rainfall intensities | tonyladson

  2. Pingback: Estimating parameters from quantiles, Pearson III distribution | tonyladson

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s