Find parameters of the Pearson III distribution given three quantile values

FindP3Params.R

# Find parameters of the Pearson III distribution given three quantile values


# Function for the Pearson III frequency factor
# details at https://tonyladson.wordpress.com/2015/03/10/frequency-factors-2/

Ky_gamma <- function(g, aep){
  
  if(abs(g) < 1e-8) { # Use the Wilson-Hilferty approximation to avoid numerical
    # issues near zero
    z <- qnorm(1-aep)
    b <- g/6
    return(z + (z^2 - 1)*b + (1/3)*(z^3 - 6*z)*b^2 - 
             (z^2 - 1)*b^3 + z*b^4 - (1/3)*(b)^5 )
  }
  
  if(g < 0) aep  <- 1 - aep
  a  <- 4/(g^2)
  b  <- 1/sqrt(a)
  sign(g)*(qgamma(1 - aep, shape = a, scale = b) - a*b)
}


# function to find g.  When this function is zero, we have a reasonable estimate of g (skew)
# Quantiles are Q1, Q2, Q3, at Annual exceedance probabilities aep1, aep2, aep3

Find_g <- function(g, Q1, Q2, Q3, aep1, aep2, aep3){
  
  # Function for the Pearson III frequency factor
  # details at https://tonyladson.wordpress.com/2015/03/10/frequency-factors-2/
  
  my.offset <- (log(Q1) - log(Q2))/(log(Q2) - log(Q3))
  (Ky_gamma(g, aep1) - Ky_gamma(g, aep2))/(Ky_gamma(g, aep2) - Ky_gamma(g, aep3)) - my.offset
  
  
}

# Calculate the PIII standard deviation parameter given g and information on two quantiles

Calc_s <- function(g, Q1, Q2, aep1, aep2){
  (log(Q1) - log(Q2))/(Ky_gamma(g, aep1) - Ky_gamma(g, aep2))
  
}

# Calculate the PIII mean parameter give g, s and information on 1 quantile

Calc_m <- function(g, s, Q, aep){
  log(Q) - s * Ky_gamma(g, aep)
}


# Example data

Q1 = 346.59
Q2 = 131.95
Q3 = 38.21
aep1 = 0.01
aep2 = 0.1
aep3 = 0.5


# Plot the graph from the blog

g.seq <- seq(-1, 1, length.out = 99)
par(oma = c(0,2,1,1))
Y <- sapply(g.seq, Find_g, Q1, Q2, Q3, aep1, aep2, aep3)

plot(g.seq, Y, 
     type = 'l',
     xlab = 'g, (skew)',
     ylab = 'f(g)',
     las = 1)
abline(h=0, lty = 2)


# find parameter values
uniroot(Find_g, c(-1,1), Q1, Q2, Q3, aep1, aep2, aep3)
# $root
# [1] -0.1131508
# 
# $f.root
# [1] -4.017024e-07
# 
# $iter
# [1] 4
# 
# $init.it
# [1] NA
# 
# $estim.prec
# [1] 6.103516e-05

Calc_s <- function(g, Q1, Q2, aep1, aep2){
  (log(Q1) - log(Q2))/(Ky_gamma(g, aep1) - Ky_gamma(g, aep2))
  
}

Calc_s(-0.1131508, Q1, Q2, aep1, aep2)
# [1] 0.9914827


Calc_m(-0.1131508, 0.9914827, Q1, aep1)
# [1] 3.624402


# Function to estimate three parameters given three quantiles and corresponding exceedance probabilities
# g.lower and g.upper specify interval for root finding



FindPIIIParams <- function(Q1, Q2, Q3, aep1, aep2, aep3, g.lower = -2, g.upper = +2){
  
  # Helper functions
  
  Ky_gamma <- function(g, aep){
    
    if(abs(g) < 1e-8) { # Use the Wilson-Hilferty approximation to avoid numerical
      # issues near zero
      z <- qnorm(1-aep)
      b <- g/6
      return(z + (z^2 - 1)*b + (1/3)*(z^3 - 6*z)*b^2 - 
               (z^2 - 1)*b^3 + z*b^4 - (1/3)*(b)^5 )
    }
    
    if(g < 0) aep  <- 1 - aep
    a  <- 4/(g^2)
    b  <- 1/sqrt(a)
    sign(g)*(qgamma(1 - aep, shape = a, scale = b) - a*b)
  }
  
  
  Find_g <- function(g, Q1, Q2, Q3, aep1, aep2, aep3){
    
    # Function for the Pearson III frequency factor
    # details at https://tonyladson.wordpress.com/2015/03/10/frequency-factors-2/
    
    my.offset <- (log(Q1) - log(Q2))/(log(Q2) - log(Q3))
    (Ky_gamma(g, aep1) - Ky_gamma(g, aep2))/(Ky_gamma(g, aep2) - Ky_gamma(g, aep3)) - my.offset
    
    
  }
  
  # Calculate the PIII standard deviation parameter given g and information on two quantiles
  
  Calc_s <- function(g, Q1, Q2, aep1, aep2){
    (log(Q1) - log(Q2))/(Ky_gamma(g, aep1) - Ky_gamma(g, aep2))
    
  }
  
  # Calculate the PIII mean parameter give g, s and information on 1 quantile
  
  Calc_m <- function(g, s, Q, aep){
    log(Q) - s * Ky_gamma(g, aep)
  }
  
  
  # Find parameters
  
  g <- uniroot(Find_g, c(g.lower,g.upper), Q1, Q2, Q3, aep1, aep2, aep3, extendInt = 'yes')$root
  s <- Calc_s(g, Q1, Q2, aep1, aep2)
  m <- Calc_m(g, s, Q1, aep1)
  
  return(data.frame(m = m, s = s, g = g))
  
  
}


FindPIIIParams(Q1 = 346.59, Q2 = 131.95, Q3 = 38.21, aep1 = 0.01, aep2 = 0.1, aep3 = 0.5)




# Check by estimating a quantile not used in derivation

# Q5% = 186.74


exp(3.6244 + 0.9914827 * Ky_gamma(-0.1131508, 0.05))
# 185.441 Ok