Find the GEV parameters used to generate design rainfalls

FindGEVparameters.R

#
# Interpolating design rainfall intensities
#

library(dplyr)
library(stringr)
library(optimx)


# Assemble data

tab.quantile <- c(64.5, 71.4, 94.1, 110.5, 127.3, 150.6, 169.4) # tabulated quantiles
aep <- c(0.6321, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01) # tabulated annual exceedance probabilities


gev.data <- data_frame(tab.quantile = tab.quantile, aep = aep) # keep together in a data frame


# Calculate Gumble reduced variate
Gumbel_rv <- function(aep) {
  
  -log(-log(1 - aep))
}

# Add to data frame

gev.data <- gev.data %>% 
  mutate(Gumbel.rv  = Gumbel_rv(aep))

# Create labels for plotting and add to data frame

my.labels <- str_c(formatC(100*aep, digits = 0, format = 'f' ), '%')
my.labels[1] <- '63.21%'

gev.data <- gev.data %>% 
   mutate(my.labels = my.labels)


# Probability plot

par(oma = c(7,2,0,0))
plot(rain.quantile ~ GEV.rv, data = gev.data,
     xlab = "Gumbel reduced variate",
     ylab = "Rainfall (mm)",
     type = 'b',
     las = 1, 
     pch = 21,
     bg = 'grey',
     ylim = c(60, 180))

# add a second axis

with(gev.data, 
  axis(side = 1, at = gev.data$GEV.rv, label = my.labels, outer = TRUE, cex.axis = 0.8)
)

mtext(text = 'AEP(%)', outer = TRUE, side = 1, line = 4)


## Find the parameters that minimimise the sum of squares between the
# calculated and tabulated quantiles



# GEV quantile function that works close to zero
# See http://stackoverflow.com/questions/36110293/floating-point-comparison-with-zero

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))
  }
}


# Sum of squares of differences between tabulated and calculated quantiles

SSQuantDiff <- function(par, tab.quantile, aep){
  calc.quantile <- sapply(aep, Qgev, par = par )
  sum((tab.quantile - calc.quantile)^2)
  
}

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

# Search for optimal parameter estimats

Param.est <- optimx(par.init, SSQuantDiff, method = "Nelder-Mead", tab.quantile = gev.data$tab.quantile, aep = gev.data$aep)

#Param.est
# p1       p2        p3       value fevals gevals niter convcode  kkt1 kkt2 xtimes
# Nelder-Mead 64.49929 18.47522 0.0885286 0.002863346    126     NA    NA        0 FALSE TRUE  0.077
# 

Param.est <- unlist(Param.est[ ,1:3]) # extract the 3 optimal parameter values: location, scale shape

# add estimates and residuals to data frame
gev.data <- gev.data %>% 
  mutate(calc.quantile = sapply(aep, Qgev, Param.est)) %>%
  mutate(quantile.resid = tab.quantile - calc.quantile)



# plot residuals
# 

par(oma = c(7,2,0,0))
plot(quantile.resid ~ Gumbel.rv, 
     data = gev.data, 
     ylim = c(-1, 1),
     ylab = 'Residual (mm)',
     pch = 21,
     bg = 'grey',
     las = 1)

abline(0, 0, lty = 2)

with(gev.data, 
     axis(side = 1, at = gev.data$Gumbel.rv, label = gev.data$my.labels, outer = TRUE, cex.axis = 0.8)
)

mtext(text = 'AEP(%)', outer = TRUE, side = 1, line = 4)


Qgev(0.3, Param.est)
# 84.44