library(dplyr)
library(lubridate)
library(fitdistrplus)
library(PearsonDS)

En este script se ajustarán las funciones de distribución Pearson tipo III y LogPearson tipo III a una serie anual de precipitaciones máximas diarias. Los métodos utilizados para obtener los parámetros del ajuste son: Método de los Momentos (MME) y Máxima Verosimilitud (MLE) mediante los paquetes fitdistrplus y PearsonDS.

Cargar datos

La estación pluviométrica con la que se trabaja corresponde a Coltauco.

Coltauco <- group_by(Coltauco, Año = year(Fecha)) %>%
  summarise(PMax = max(valor, na.rm = TRUE))

plot(Coltauco, type = "h", xlab = "Año", ylab = "Precipitación máxima diaria [mm]")

Ajuste Pearson III por método de máxima verosimilitud (MLE)

Funciones d, p y q

Estas funciones están incluidas en el paquete PearsonDS.

dPIII <- function(x, shape, location, scale) dpearsonIII(x, shape, location, scale, log = FALSE)
pPIII <- function(q, shape, location, scale) ppearsonIII(q, shape, location, scale, lower.tail = TRUE, log.p = FALSE)
qPIII <- function(p, shape, location, scale) qpearsonIII(p, shape, location, scale, lower.tail = TRUE, log.p = FALSE)

Estadígrafos

# Note that the above forgot to mulitply the scale by the sign of skewness .
# Refer to Page 24 of the Bulletin 17c

m <- mean(Coltauco$PMax)
v <- var(Coltauco$PMax)
s <- sd(Coltauco$PMax)
g <- e1071::skewness(Coltauco$PMax, type = 2)

Valores iniciales para método MLE

# This can be corrected, but HEC Bulletin 17b does not do these corrections
# Correct the sample skew for bias using the recommendation of
# Bobee, B. and R. Robitaille (1977). "The use of the Pearson Type 3 and Log Pearson Type 3 distributions revisited."
# Water Resources Reseach 13(2): 427-443, as used by Kite

n <- length(Coltauco$PMax)
g <- g*(sqrt(n*(n-1))/(n-2))*(1+8.5/n)

# We will use method of moment estimates as starting Values for the MLE search

my_shape <- (2 / g) ^ 2
my_scale <- sqrt(v / my_shape) * sign(g)
my_location <- m - my_scale * my_shape

start <- list(shape = my_shape, location = my_location, scale = my_scale)

Se ajusta la función de distribución

PIII <- fitdist(Coltauco$PMax, "PIII", start = start, control = list(maxit = 1000))

## $start.arg
## $start.arg$shape
## [1] 18.83867
## 
## $start.arg$location
## [1] -75.47069
## 
## $start.arg$scale
## [1] 8.577303
## 
## 
## $fix.arg
## NULL

Para verificar que está todo bien, se procede a hacer los gráficos de densidad, densidad acumulada, P-P y Q-Q.

plot(PIII)

Ajuste LogPearson III por método de máxima verosimilitud (MLE)

Este caso es casi idéntico al anterior. Se debe transformar la variable y sus estadígrafos a Log10 o Ln y hacer ciertos cambios en las funciones del paquete PearsonDS.

Funciones d, p y q

Nótese el cambio logarítmico en el argumento de p y p, además de dividir d por la variable x y trabaajr con el exponencial de la función q.

dlPIII <- function(x, shape, location, scale) dpearsonIII(log(x), shape, location, scale, log=FALSE)/x
plPIII <- function(q, shape, location, scale) ppearsonIII(log(q), shape, location, scale, lower.tail = TRUE, log.p = FALSE)
qlPIII <- function(p, shape, location, scale) exp(qpearsonIII(p, shape, location, scale, lower.tail = TRUE, log.p = FALSE))

Estadígrafos

m <- mean(log(Coltauco$PMax))
v <- var(log(Coltauco$PMax))
s <- sd(log(Coltauco$PMax))
g <- e1071::skewness(log(Coltauco$PMax), type = 2)

Valores iniciales para método MLE

#n <- length(Coltauco$PMax)
#g <- g*(sqrt(n*(n-1))/(n-2))*(1+8.5/n)

my_shape <- (2 / g) ^ 2
my_scale <- sqrt(v) / sqrt(my_shape) * sign(g)
my_location <- m - my_scale * my_shape

start <- list(shape = my_shape, location = my_location, scale = my_scale)

Se ajusta la función de distribución

LPIII <- fitdist(Coltauco$PMax, "lPIII", start = start, control = list(maxit = 1000))

## $start.arg
## $start.arg$shape
## [1] 41.86267
## 
## $start.arg$location
## [1] 7.368915
## 
## $start.arg$scale
## [1] -0.07196876
## 
## 
## $fix.arg
## NULL

plot(LPIII)

Ajuste Pearson III por método de los momentos (MME)

Estadígrafos

m <- mean(Coltauco$PMax)
v <- var(Coltauco$PMax)
s <- sd(Coltauco$PMax)
g <- e1071::skewness(Coltauco$PMax, type = 2)

Valores iniciales

my_shape <- (2 / g) ^ 2
my_scale <- sqrt(v / my_shape) * sign(g)
my_location <- m - my_scale * my_shape

start <- list(shape = my_shape, location = my_location, scale = my_scale)

Funciones generadoras de momentos

En estos casos, el paquete fitdistrplus requiere que se ingresen, como argumento, las funciones de momentos empíricos y teóricos de la distribución que se busca ajustar. Estas funciones (mPIII y ePIII) las obtuve del código del paquete fasstr.

mPIII <- function(order, shape, location, scale) # compute the empirical first 3 moments of the PIII distribution
  {
  if(order == 1) return(location + shape * scale)
  if(order == 2) return(scale * scale * shape)
  if(order == 3) return(2 / sqrt(shape) * sign(scale))
  }

ePIII <- function(x, order) # compute (centered) empirical centered moments of the data
  {
  if(order == 1) return(mean(x))
  if(order == 2) return(var(x))
  if(order == 3) return(e1071::skewness(x, type = 2))
  }

Ajuste de función

PIIIMME <- fitdist(Coltauco$PMax, "PIII", start = start, method = "mme", order = 1:3, memp = ePIII, control = list(maxit = 1000))

## $start.arg
## $start.arg$shape
## [1] 29.31214
## 
## $start.arg$location
## [1] -115.4435
## 
## $start.arg$scale
## [1] 6.876256
## 
## 
## $fix.arg
## NULL

plot(PIIIMME)

Ajuste LogPearson III por método de los momentos (MME)

Estadígrafos

m <- mean(log(Coltauco$PMax))
v <- var(log(Coltauco$PMax))
s <- sd(log(Coltauco$PMax))
g <- e1071::skewness(log(Coltauco$PMax), type = 2)

Valores iniciales

my_shape <- (2 / g) ^ 2
my_scale <- sqrt(v / my_shape) * sign(g)
my_location <- m - my_scale * my_shape

start <- list(shape = my_shape, location = my_location, scale = my_scale)

Funciones generadoras de momentos (logarítmicas)

mlPIII <- function(order, shape, location, scale)
  {
  # compute the empirical first 3 moments of the PIII distribution
  if(order == 1) return(location + shape * scale)
  if(order == 2) return(scale * scale * shape)
  if(order == 3) return(2 / sqrt(shape) * sign(scale))
  }

elPIII <- function(x, order)
  {
  # compute (centered) empirical centered moments of the data
  if(order == 1) return(mean(log(x)))
  if(order == 2) return(var(log(x)))
  if(order == 3) return(e1071::skewness(log(x), type = 2))
  }

Ajuste de distribución

LPIIIMME <- fitdist(Coltauco$PMax, "lPIII", start = start, method = "mme", order = 1:3, memp = elPIII, control = list(maxit = 1000))

## $start.arg
## $start.arg$shape
## [1] 41.86267
## 
## $start.arg$location
## [1] 7.368915
## 
## $start.arg$scale
## [1] -0.07196876
## 
## 
## $fix.arg
## NULL

plot(LPIIIMME)

Comparación entre ajustes

Se procede a comparar los ajustes, de forma visual y numérica, para percatarnos de las diferencias.

Gráficos de densidad

par(mfrow=c(2,2))

denscomp(PIII, addlegend = FALSE, main = "Pearson 3 MLE")
denscomp(PIIIMME, addlegend = FALSE, main = "Pearson 3 MME")
denscomp(LPIII, addlegend = FALSE, main = "LogPearson 3 MLE")
denscomp(LPIIIMME, addlegend = FALSE, main = "LogPearson 3 MME")

Gráficos de densidad acumulada

par(mfrow=c(2,2))

cdfcomp(PIII, addlegend = FALSE, main = "Pearson 3 MLE")
cdfcomp(PIIIMME, addlegend = FALSE, main = "Pearson 3 MME")
cdfcomp(LPIII, addlegend = FALSE, main = "LogPearson 3 MLE")
cdfcomp(LPIIIMME, addlegend = FALSE, main = "LogPearson 3 MME")

Gráficos Q-Q

par(mfrow=c(2,2))

qqcomp(PIII, addlegend = FALSE, main = "Pearson 3 MLE")
qqcomp(PIIIMME, addlegend = FALSE, main = "Pearson 3 MME")
qqcomp(LPIII, addlegend = FALSE, main = "LogPearson 3 MLE")
qqcomp(LPIIIMME, addlegend = FALSE, main = "LogPearson 3 MME")

Gráficos P-P

par(mfrow=c(2,2))

ppcomp(PIII, addlegend = FALSE, main = "Pearson 3 MLE")
ppcomp(PIIIMME, addlegend = FALSE, main = "Pearson 3 MME")
ppcomp(LPIII, addlegend = FALSE, main = "LogPearson 3 MLE")
ppcomp(LPIIIMME, addlegend = FALSE, main = "LogPearson 3 MME")

Gráficos combinados para Pearson tipo III

denscomp(list(PIII, PIIIMME))
cdfcomp(list(PIII, PIIIMME))
qqcomp(list(PIII, PIIIMME))
ppcomp(list(PIII, PIIIMME))

Gráficos combinados para LogPearson tipo III

denscomp(list(LPIII, LPIIIMME))
cdfcomp(list(LPIII, LPIIIMME))
qqcomp(list(LPIII, LPIIIMME))
ppcomp(list(LPIII, LPIIIMME))

En el caso de las distribuciones Pearson III, no se observan grandes variaciones en los gráficos P-P, Q-Q y CDF. Por el contrario, en el caso de las distribuciones LogPearson III, hay mayores diferencias en todos los gráficos.

Magnitudes de retorno

Se compararán las magnitudes de retorno entregadas por cada distribución, para periodos de retorno de 2 a 500 años.

T <- 2:500
Pexc <- 1/T

Magnitudes <- data.frame(T,
                         Pexc,
                         qPIII(1 - Pexc, PIII$estimate[1], PIII$estimate[2], PIII$estimate[3]),
                         qPIII(1 - Pexc, PIIIMME$estimate[1], PIIIMME$estimate[2], PIIIMME$estimate[3]),
                         qlPIII(1 - Pexc, LPIII$estimate[1], LPIII$estimate[2], LPIII$estimate[3]),
                         qlPIII(1 - Pexc, LPIIIMME$estimate[1], LPIIIMME$estimate[2], LPIIIMME$estimate[3]))

colnames(Magnitudes) <- c("T", "Pexc", "P3MLE", "P3MME", "LP3MLE", "LP3MME")

matplot(Magnitudes[,c(3:6)], type = c("l"), col = 1:4, ylab = "Precipitación Máxima Diaria [mm]", xlab = "T [años]")

legend("bottomright", 
       legend = c("Pearson III MLE", "Pearson III MME", "LogPearson III MLE", "LogPearson III MME"), 
       col=1:4, pch=15, cex = 0.8)

Según Rivano, Fulvio (2004), el método de los momentos entrega estimadores más robustos para las distribuciones Pearson III y LogPearson III, por lo que se selecciona este método para hacer los ajustes a las series hidrológicas de máximos anuales.

Ajuste de distribuciones Pearson III y Log Pearson III

Fernando Carreño

17-02-2022

Cargar datos

Ajuste Pearson III por método de máxima verosimilitud (MLE)

Funciones d, p y q

Estadígrafos

Valores iniciales para método MLE

Ajuste LogPearson III por método de máxima verosimilitud (MLE)

Funciones d, p y q

Estadígrafos

Valores iniciales para método MLE

Ajuste Pearson III por método de los momentos (MME)

Estadígrafos

Valores iniciales

Funciones generadoras de momentos

Ajuste LogPearson III por método de los momentos (MME)

Estadígrafos

Valores iniciales

Funciones generadoras de momentos (logarítmicas)

Comparación entre ajustes

Gráficos de densidad

Gráficos de densidad acumulada

Gráficos Q-Q

Gráficos P-P

Gráficos combinados para Pearson tipo III

Gráficos combinados para LogPearson tipo III

Magnitudes de retorno