Interpolating design rainfall intensities

Background:  The design rainfall data from the Bureau of Meteorology provides rainfall depths for a standard set of seven Annual Exceedance Probabilities.   These are AEPs of 1%, 2%, 5%, 10%, 20%, 50% and 1EY (1EY means one exceedance per year).  For any particular location in Australia, data are supplied as a table with durations as rows and the the seven standard AEPs as columns (Table 1).

These design rainfall depths have been determined by the Bureau by fitting a GEV distribution to rainfall data (see ARR Book 2, Chapter 3).  A separate GEV model has been developed for each of the standard durations (the rows of Table 1).  To calculate rainfall depths other than at the standard durations and AEPs we need to interpolate the values in the table.  In this post I am looking at how to find rainfall depths for AEPs other than the 7 standard values.  That is, I’m restricting the discussion here to the standard durations i.e. interpolating between the columns, not between the rows.  I’ll look at interpolating between durations in a later post.

Example: Let’s look at how we could estimate the 24 hour rainfall depth for the 30% Annual Exceedance Probability.  As shown in Table 1, no rainfall depths are provided for an AEP of 30% but there are 7 quantile estimates for each of the 7 annual exceedance probabilities. For example, the 24 hour, 1% rainfall depth is 169.4 mm, the 24 hour 2% rainfall is 150.6 mm etc.

Thomson-RainfallIFD copy 2

Table 1: design rainfall depths

To calculate rainfall depths at other exceedance probabilities we need to work out the parameters of the GEV distribution that the Bureau has used and then calculate the required quantiles.  In a previous post I developed a general method to calculate GEV parameters given 3 quantiles.  We could use that method, but here I want to take a different tack.  For the current problem, we have two pieces of specific information:

  • We know the quantile for an exceedance probability of 1EY (63.21%).  This provides a direct estimate of the location parameter of GEV
  • There is information available for seven quantiles (rather than the case of three quantiles as discussed in the previous post).  Using a numerical search, we can find the GEV parameters that provide a best fit to these seven quantiles.

The quantile function for the GEV distribution is:

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

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

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.  Where the shape parameter is zero, the GEV simplifies to the Gumbel distribution (equation 1b above).

As shown previously, it is possible to estimate the parameters of the GEV distribution using three quantiles, but the Bureau has kindly provided 7 quantiles.  If we can get good initial estimates of the parameter values, a numerical optimisation approach can be used to polish the estimates to provide a best fit to all the available data.

A probability plot of these quantiles provides some insight.  For the Gumbel distribution a lot of Q against -log[-log(1 – p)] (referred to as the Gumbel reduced variate) will plot as a straight line with a slope of \alpha .  For the GEV this line will be curved because of the influence of \kappa .

The probability plot for the 24 hour rainfall depths from Table 1 is shown below.

Gumbel_rv

Figure 1: Probability plot for the 24-hour rainfall depths from Table 1

The plot can be used to suggest some initial values for the GEV parameters.

First, the plotted line is reasonably straight which suggests that \kappa is near zero.

Second, consider the quantile with an Annual Exceedance Probability of “1EY” i.e. 1 exceedance per year.  According to Table 1.2.1 in ARR (Book 1, Section 2), 1EY corresponds to an AEP of 63.21 %; this comes from 1 – 1/e.

Substitute p = 1 - 1/e into equation 1 and all terms other than \xi are zero i.e.

Q(1EY) = \xi

This is the point where the Gumbel reduced variate is zero on Figure 1.

Therefore, the column headed “1EY” in Table 1 provides an estimate of the location parameter for all of the standard durations.  For the 24 hour rainfall plotted in Figure 1, the location parameter is about 64.5 mm.

Third, the \alpha parameter is the slope of Figure 1 near where the Gumbel reduced variate is zero.  Let’s calculate the slope between an AEP of 1EY and an AEP of 50%.  There are the two leftmost points in Figure 1.

\alpha \approx \frac{71.4 - 64.5}{-log[-log(1 - 0.5)]} = 18.83

So we have initial estimates of the three GEV parameters for 24-hour rainfall totals:

  • location (\xi) = 64.5
  • scale (\alpha ) = 18.83
  • shape (\kappa ) = 0

We can use these initial estimates to find optimal parameter values.  Here I’m defining optimal in the sense of minimising the sum of the squared residuals between calculated and tabulated quantiles.

First a quantile function is required that provides sensible estimates near zero values of shape (\kappa ).  See this discussion for details.

Qgev <- function(p, par){
 # p is an exceedance probability
 
 location <- par[1]
 scale <- par[2]
 k <- par[3]
 
 .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))
 }
}

Next a function is required that will require the sum of squares of the difference between tabulated and calculated quantiles.

SSQuantDiff <- function(par, tab.quantile, aep){
# par = parameters of the GEV distribution, location, scale shape
# tab.quantile = vector of the tabulated quantiles
# aep = vector of annual exceedance probabilities corresponding to the quantiles
 
 calc.quantile <- sapply(aep, Qgev, par = par )
 sum((tab.quantile - calc.quantile)^2)
 
}

 

Now use optimisation to find the best parameter values (i.e. the minimise the sum of squares).  Here I use the optimx function in R.

library(optimx)

par.init <- c(64.5, 18.83, 0) # initial parameter values

Param.est <- optimx(par.init, SSQuantDiff, 
  method = "Nelder-Mead", 
  tab.quantile = gev.data$tab.quantile, # gev.data stores the tabulated quantiles and aeps
  aep = gev.data$aep)

Param.est

# p1 p2 p3 
# 64.49929 18.47522 0.0885286 
# 

So are best estimates are:

  • location (\xi) = 64.499
  • scale (\alpha ) = 18.475
  • shape (\kappa ) = 0.0885

The residuals, when using these parameters, are shown in Figure 2.  All the residuals are much less than 1 mm so using the parameters for interpolation is fine and a moderate amount of extrapolation is probably ok.  For very rare rainfalls, a different approach should be used which is discussed here.

GEV-resid

Figure 2: Residuals (tabulated – calculated quantiles)

 

These parameters can be used to estimate rainfall with a 30% AEP which was the focus of the example at the start of this post.

Param.est <- unlist(Param.est[ ,1:3])
Qgev(0.3, Param.est)
# 84.44

The 30% AEP 24 hour rainfall is 84.4 mm.

All the code is available as a gist.

One thought on “Interpolating design rainfall intensities

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