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:

(1)

Where is the quantile, is an exceedance probability, is the location parameter, is the scale parameter and 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 at exceedance probabilities and that

(2)

(3)

(4)

To eliminate subtract (3) from (2).

(5)

Similarly, remove by division to achieve an equation in alone.

(5)

Equation (5) can be solved numerically to obtain then substitute back into (4) to obtain and then into (1) to obtain

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

(6)

Note that this has a singularity at . We’ll address this by interpolating equation (6) through zero based on positive and negative 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 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 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 () = 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 = 0. This can cause numerical issues as 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%.

Pingback: Interpolating design rainfall intensities | tonyladson

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