*Australian Rainfall and Runoff* suggests the following formula when calculating confidence limits for estimates from flood frequency analysis.

Where:

CL = confidence limit

*Q _{y}* = discharge estimated from flood frequency analysis for the selected AEP of 1 in Y years

*F*= Frequency factor for the Normal distribution for the selected probability level

= parameter for determining the standard error of the Pearson III distribution

*S*= standard deviation of the logarithms of flows to base 10

*N*= number of floods used in analysis

The formula is based on work by Kite (1975; 1976; 2004).

To use the formula we need values for . These are tabulated in Australian Rainfall and Runoff (see below) but it would be much more convenient if they could be calculated.

The formula used to determine the values is given in Kite (2004), equation 9-54 and is based on the Pearson III frequency factor and its derivative with respect to skewness.

Where:

*K* = frequency factor for a given skew and AEP

*g* = skew

We can use the frequency factor function, developed earlier, and its numerical derivative to calculate these values. Alternatively, we can estimate the derivative based on the Wilson-Hilferty approximation.

Where, *z* is the standard normal deviate.

See the code below for functions that calculate delta based on these two approaches: 1) using a numerical derivative, and 2) using the Wilson-Hilferty approximation.

```
library(lmomco)
library(numDeriv)
Ky_Qp3 <- function(g, aep){
# Return the quantile of the Pearson III distribution given
# skewness of g, and annual exceedance probability, aep
param <- vec2par(c(0,1,g), type='pe3')
quape3(1 - aep, param)
}
Delta <- function(g, aep){
K <- Ky_Qp3(g, aep)
dK <- grad(Ky_Qp3, g, method="Richardson", aep = aep)
sqrt(1 + K*g + 0.5*(K^2)*(0.75*(g^2) + 1) + 3*K*dK*(g + 0.25*g^3) + 3*(dK^2)*(2 + 3*g^2 + (5/8)*g^4))
}
# Delta value using the Wilson-Hilferty approximation
Delta_wh <- function(g, aep){
K <- Ky_Qp3(g, aep)
z <- qnorm(1 - aep)
dK <- (z^2 - 1)/6 + g * 4*(z^3 - 6*z)/(6^3) - g^2 * 3*(z^2 - 1)/6^3 + g^3 * 4*z/(6^4) - g^4 * 10/(6^6)
sqrt(1 + K*g + 0.5*(K^2)*(0.75*(g^2) + 1) + 3*K*dK*(g + 0.25*g^3) + 3*(dK^2)*(2 + 3*g^2 + (5/8)*g^4))
}
```

Let’s test the calculated delta values against those tabulated in Australian Rainfall and Runoff. I’ve written a gist, DeltaTable, that will return the tabulated values and interpolate if necessary.

```
# Testing the function
library(lmomco)
library(numDeriv)
library(dplyr)
library(devtools)
library(e1071)
# Test Delta against tabulated values
# Source DeltaTable function
source_gist('https://gist.github.com/TonyLadson/08d6a1834d6f395c725e')
# We also need the functions: Ky_Qp3, Delta and Delta_wh from the code above
Delta(1, 0.01)
# [1] 5.137346
Delta_wh(1, 0.01)
# [1] 5.167474
DeltaTable(1, 0.01)
# [1] 5.1499 # So there are differences
# Test for a range of skews and aeps
my.skew <- seq(0,3.9,0.1)
my.aep <- c(0.998, 0.995, 0.990, 0.980, 0.950, 0.900, 0.800, 0.500, 0.200, 0.100, 0.050, 0.020, 0.010, 0.005, 0.002)
comparison.df <- tbl_df(expand.grid(my.skew = my.skew, my.aep = my.aep))
comparison.df %>%
mutate(delta.table = mapply(DeltaTable, my.skew, my.aep)) %>%
mutate(delta.calc = mapply(Delta, my.skew, my.aep)) %>%
mutate(pc.diff = round(abs(100 * (delta.table - delta.calc)/delta.table), 2)) %>%
arrange(desc(pc.diff)) # percentage difference
# my.skew my.aep delta.table delta.calc pc.diff
# (dbl) (dbl) (dbl) (dbl) (dbl)
# 1 3.7 0.002 30.5310 29.753060 2.55
# 2 3.0 0.002 23.7706 23.166982 2.54
# 3 3.7 0.005 20.5459 20.970925 2.07
# 4 0.4 0.998 3.1156 3.168725 1.71
# 5 3.8 0.002 31.2357 30.717047 1.66
# 6 3.9 0.005 22.4853 22.116728 1.64
# 7 2.9 0.002 21.9173 22.257135 1.55
# 8 3.5 0.002 28.2778 27.839306 1.55
# 9 2.5 0.002 18.4580 18.724809 1.45
# 10 2.8 0.002 21.6578 21.357200 1.39
# .. ... ... ... ... ...
# For the Wilson Hilferty version we'll do the comparison in the region where the
# approximation is considered reasonable and for values that are used in practice
# Test for a range of skews and aeps
my.skew <- seq(0,2,0.1)
my.aep <- c(0.200, 0.100, 0.050, 0.020, 0.010, 0.005)
comparison.df <- tbl_df(expand.grid(my.skew = my.skew, my.aep = my.aep))
comparison.df %>%
mutate(delta.table = mapply(DeltaTable, my.skew, my.aep)) %>%
mutate(delta.calc = mapply(Delta_wh, my.skew, my.aep)) %>%
mutate(pc.diff = round(abs(100 * (delta.table - delta.calc)/delta.table), 2)) %>%
arrange(desc(pc.diff)) # percentage difference
# Source: local data frame [126 x 5]
#
# my.skew my.aep delta.table delta.calc pc.diff
# (dbl) (dbl) (dbl) (dbl) (dbl)
# 1 2.0 0.05 3.8348 3.691666 3.73
# 2 1.9 0.05 3.7425 3.616571 3.36
# 3 2.0 0.02 6.3943 6.197861 3.07
# 4 1.8 0.05 3.6483 3.536904 3.05
# 5 1.9 0.02 6.1755 6.000034 2.84
# 6 1.7 0.05 3.5501 3.452598 2.75
# 7 1.6 0.05 3.4482 3.363735 2.45
# 8 1.8 0.02 5.9384 5.795547 2.41
# 9 1.5 0.05 3.3430 3.270531 2.17
# 10 1.7 0.02 5.7001 5.585454 2.01
# .. ... ... ... ... ...
```

The largest difference, using the numerical derivative approach, is about 2.55% which is probably because of the difference in the calculation method that I have applied compared with that used to generate the table in Australian Rainfall and Runoff. The Wilson-Hilferty version is ok for a smaller range of skews and AEPs.

The upshot is that the function `Delta`

is fit for purpose to calculate delta values for use in flood frequency analysis.

Edit 29 Feb 2016.

An additional source of delta values is Table 9-3 in Kite (2004) . These values are very similar to those described above but it is useful to have the actual values when checking your calculations with those published in Kite. Values are available via this gist.

require(devtools) source_gist('https://gist.github.com/TonyLadson/280fd3ae3e730b3bb370') DeltaTable_Kite(g = 0.1, aep = 0.02) # [1] 2.3425