Creación de una página de ayuda personalizada de una función de R: «Ejemplo ecf (inversa percentil)»

  Por

En esta entrada se muestra el ejemplo de creación de una ayuda personalizada del uso de una función de R, en este caso de la función ecf(). Para ello, le he pedido a Calude Sonnet 4.5 el siguiente prompt:

Genera un documento quarto con R con el título 
Ejemplo ecf (inversa percentil)

Que sea un ejemplo sencillo del uso de ecf en R para calcular la inversa del percentil de una serie de caudles (responder a la pregunta ¿un caudal Q, a qué percentil corresponde de la serie X?
Explora las difernetes formas de cálculo de la función de distribución (especialmente entre las escalonadas y las interpoladas).
También algo de gráficas con ggplot

El resultado es un fichero de Quarto con el que se ha generado el siguiente html:

El Continuación se muestra el código de Quarto generado por Claude:

---
title: "Ejemplo ecf (inversa percentil)"
author: ""
format: 
  html:
    toc: true
    toc-depth: 3
    code-fold: show
    theme: cosmo
lang: es
---

## Introducción

En hidrología, es frecuente necesitar determinar **a qué percentil corresponde un caudal determinado** dentro de una serie histórica. Por ejemplo, para clasificar un evento como percentil 90, identificar caudales de estiaje (percentiles bajos) o evaluar la posición relativa de un caudal observado.

La función de distribución empírica acumulada (ecdf, *empirical cumulative distribution function*) es la herramienta natural para esta tarea. Sin embargo, existen diferentes métodos para calcular percentiles, especialmente al invertir la función de distribución.

Este documento explora:

1.  El uso básico de `ecdf()` en R
2.  Diferentes métodos de interpolación para el cálculo de percentiles
3.  La diferencia entre funciones escalonadas e interpoladas
4.  Aplicaciones prácticas con series de caudales

## Configuración inicial

```{r setup, message=FALSE, warning=FALSE}
library(tidyverse)
library(scales)

# Configuración de tema para gráficos
theme_set(theme_minimal(base_size = 11))
```

## Generación de datos de ejemplo

Creamos una serie sintética de caudales diarios para un río, con características hidrológicas realistas:

```{r datos}
set.seed(123)

# Serie de caudales (m³/s) con distribución log-normal
# Simulando un río con caudal medio de 15 m³/s
n_dias <- 1000

caudales <- tibble(
  fecha = seq.Date(from = as.Date("2020-01-01"), 
                   by = "day", 
                   length.out = n_dias),
  caudal = rlnorm(n_dias, meanlog = log(15), sdlog = 0.8)
) %>%
  arrange(caudal)

# Resumen estadístico
summary(caudales$caudal)
```

## Función de distribución empírica: ecdf()

La función `ecdf()` de R base crea una función escalonada que representa la probabilidad empírica acumulada:

```{r ecdf-basico}
# Crear la función ecdf
Fn <- ecdf(caudales$caudal)

# Ejemplos de uso: ¿Qué percentil corresponde a cada caudal?
caudales_prueba <- c(5, 10, 15, 20, 30, 50)
percentiles <- Fn(caudales_prueba) * 100

resultados_ecdf <- tibble(
  caudal = caudales_prueba,
  percentil = percentiles
)

print(resultados_ecdf)
```

La función `ecdf()` es **escalonada**: aumenta en saltos discretos en cada valor observado de la muestra.

## Métodos de cálculo de percentiles

R ofrece 9 métodos diferentes para calcular percentiles (cuantiles) mediante la función `quantile()`. Los más relevantes son:

-   **Tipo 1**: Función inversa de la ecdf (discontinua)
-   **Tipo 4**: Interpolación lineal de la ecdf empírica
-   **Tipo 5**: Interpolación lineal ponderada (hidrólogos)
-   **Tipo 7**: Interpolación lineal (por defecto en R)
-   **Tipo 8**: Interpolación lineal aproximadamente insesgada
-   **Tipo 9**: Interpolación con cuantiles insesgados

```{r metodos-quantile}
# Definir percentiles de interés
p_interes <- c(0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99)

# Calcular con diferentes métodos
metodos_comp <- map_dfr(1:9, function(tipo) {
  tibble(
    metodo = tipo,
    percentil = p_interes * 100,
    caudal = quantile(caudales$caudal, probs = p_interes, type = tipo)
  )
})

# Comparación de métodos para percentiles clave
metodos_comp %>%
  pivot_wider(names_from = metodo, 
              values_from = caudal,
              names_prefix = "Tipo_") %>%
  knitr::kable(digits = 2, 
               caption = "Caudales (m³/s) según método de cálculo")
```

## Inversión de la ecdf: de caudal a percentil

Para invertir la relación (dado un caudal Q, ¿qué percentil es?), necesitamos interpolar:

```{r inversion-ecdf}
# Función para calcular percentil de un caudal dado
calcular_percentil <- function(Q, datos, metodo = 7) {
  # Ordenar datos
  x_ord <- sort(datos)
  n <- length(x_ord)
  
  # Si Q está fuera del rango
  if (Q <= min(x_ord)) return(0)
  if (Q >= max(x_ord)) return(100)
  
  # Buscar posición mediante interpolación lineal de quantile()
  # Usamos approx() para interpolar la función inversa
  probs <- seq(0, 1, length.out = 100)
  cuantiles <- quantile(x_ord, probs = probs, type = metodo)
  
  percentil <- approx(x = cuantiles, y = probs * 100, xout = Q)$y
  
  return(percentil)
}

# Aplicar a varios caudales de prueba
caudales_test <- seq(5, 60, by = 5)

resultados_inversion <- tibble(
  caudal = caudales_test,
  percentil_tipo4 = map_dbl(caudales_test, 
                             ~calcular_percentil(., caudales$caudal, metodo = 4)),
  percentil_tipo7 = map_dbl(caudales_test, 
                             ~calcular_percentil(., caudales$caudal, metodo = 7)),
  percentil_tipo8 = map_dbl(caudales_test, 
                             ~calcular_percentil(., caudales$caudal, metodo = 8))
)

print(resultados_inversion)
```

## Visualización: Funciones escalonada vs interpolada

```{r grafico-comparacion, fig.width=10, fig.height=6}
# Preparar datos para gráfico
n_puntos <- length(caudales$caudal)
datos_plot <- tibble(
  caudal = sort(caudales$caudal),
  percentil_escalonado = (1:n_puntos) / n_puntos * 100
)

# Calcular versión interpolada (tipo 7)
percentiles_interp <- seq(0, 100, length.out = 500)
caudales_interp <- quantile(caudales$caudal, 
                             probs = percentiles_interp/100, 
                             type = 7)

datos_interpolados <- tibble(
  caudal = caudales_interp,
  percentil_interpolado = percentiles_interp
)

# Gráfico comparativo
ggplot() +
  # Función escalonada (ecdf)
  geom_step(data = datos_plot,
            aes(x = caudal, y = percentil_escalonado, 
                color = "Escalonada (ecdf)"),
            linewidth = 0.8) +
  # Función interpolada
  geom_line(data = datos_interpolados,
            aes(x = caudal, y = percentil_interpolado, 
                color = "Interpolada (tipo 7)"),
            linewidth = 0.8) +
  # Líneas de referencia para percentiles clave
  geom_hline(yintercept = c(10, 50, 90), 
             linetype = "dashed", 
             alpha = 0.3) +
  scale_color_manual(values = c("Escalonada (ecdf)" = "#2C3E50",
                                 "Interpolada (tipo 7)" = "#E74C3C")) +
  scale_x_continuous(labels = comma) +
  labs(
    title = "Comparación: Función de distribución escalonada vs interpolada",
    subtitle = "Serie de caudales diarios (n = 1000)",
    x = "Caudal (m³/s)",
    y = "Percentil (%)",
    color = "Método"
  ) +
  theme(legend.position = "bottom")
```

### Detalle en percentiles extremos

```{r grafico-extremos, fig.width=10, fig.height=5}
# Zoom en percentiles bajos (estiaje)
p1 <- ggplot() +
  geom_step(data = datos_plot %>% filter(percentil_escalonado <= 20),
            aes(x = caudal, y = percentil_escalonado),
            color = "#2C3E50", linewidth = 1) +
  geom_line(data = datos_interpolados %>% filter(percentil_interpolado <= 20),
            aes(x = caudal, y = percentil_interpolado),
            color = "#E74C3C", linewidth = 1, linetype = "dashed") +
  labs(title = "Percentiles bajos (estiaje)",
       x = "Caudal (m³/s)",
       y = "Percentil (%)") +
  theme_minimal()

# Zoom en percentiles altos (avenidas)
p2 <- ggplot() +
  geom_step(data = datos_plot %>% filter(percentil_escalonado >= 80),
            aes(x = caudal, y = percentil_escalonado),
            color = "#2C3E50", linewidth = 1) +
  geom_line(data = datos_interpolados %>% filter(percentil_interpolado >= 80),
            aes(x = caudal, y = percentil_interpolado),
            color = "#E74C3C", linewidth = 1, linetype = "dashed") +
  labs(title = "Percentiles altos (crecidas)",
       x = "Caudal (m³/s)",
       y = "Percentil (%)") +
  theme_minimal()

# Combinar gráficos
library(patchwork)
p1 + p2
```

## Ejemplo práctico: Clasificación de caudales

Clasificamos los caudales según su percentil:

```{r clasificacion}
# Función para clasificar caudal según percentil
clasificar_caudal <- function(Q, datos) {
  percentil <- calcular_percentil(Q, datos, metodo = 7)
  
  categoria <- case_when(
    percentil < 10 ~ "Estiaje extremo",
    percentil < 25 ~ "Estiaje",
    percentil < 75 ~ "Normal",
    percentil < 90 ~ "Crecida",
    percentil < 95 ~ "Crecida importante",
    TRUE ~ "Avenida excepcional"
  )
  
  return(list(percentil = percentil, categoria = categoria))
}

# Ejemplos de clasificación
caudales_ejemplo <- c(3, 8, 15, 25, 40, 60)

clasificaciones <- map_dfr(caudales_ejemplo, function(Q) {
  resultado <- clasificar_caudal(Q, caudales$caudal)
  tibble(
    caudal = Q,
    percentil = round(resultado$percentil, 1),
    categoria = resultado$categoria
  )
})

print(clasificaciones)
```

## Curva de duración de caudales

Una aplicación clásica de la ecdf en hidrología:

```{r curva-duracion, fig.width=10, fig.height=6}
# Curva de duración (ordenar de mayor a menor)
caudales_duracion <- caudales %>%
  arrange(desc(caudal)) %>%
  mutate(
    orden = row_number(),
    prob_excedencia = (orden) / (n() + 1) * 100
  )

ggplot(caudales_duracion, aes(x = prob_excedencia, y = caudal)) +
  geom_line(color = "#3498DB", linewidth = 1) +
  geom_hline(yintercept = quantile(caudales$caudal, c(0.10, 0.50, 0.90)),
             linetype = "dashed", alpha = 0.5, color = "#E74C3C") +
  annotate("text", x = 95, y = quantile(caudales$caudal, 0.90) * 1.1,
           label = "P90", color = "#E74C3C") +
  annotate("text", x = 95, y = quantile(caudales$caudal, 0.50) * 1.1,
           label = "Mediana", color = "#E74C3C") +
  annotate("text", x = 95, y = quantile(caudales$caudal, 0.10) * 1.1,
           label = "P10", color = "#E74C3C") +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_log10(labels = comma) +
  labs(
    title = "Curva de duración de caudales",
    subtitle = "Escala logarítmica en eje Y",
    x = "Probabilidad de excedencia (%)",
    y = "Caudal (m³/s)"
  )
```

## Conclusiones

1.  **La función `ecdf()`** proporciona la función de distribución empírica acumulada, que es escalonada por naturaleza.

2.  **Para invertir la relación** (de caudal a percentil), es recomendable usar métodos interpolados, especialmente el **tipo 7** (por defecto en R) o el **tipo 8** (aproximadamente insesgado).

3.  **La diferencia entre métodos** es más notable en los extremos de la distribución y en muestras pequeñas.

4.  **En hidrología práctica**, para series largas (n \> 100), las diferencias entre métodos son generalmente menores que las incertidumbres de medición.

5.  **Aplicaciones principales**:

    -   Clasificación de eventos (estiaje, avenida)
    -   Curvas de duración
    -   Análisis de frecuencias
    -   Evaluación de garantías de suministro

## Referencias

-   Hyndman, R.J. y Fan, Y. (1996). *Sample Quantiles in Statistical Packages*. The American Statistician, 50(4), 361-365.
-   Helsel, D.R. et al. (2020). *Statistical Methods in Water Resources*. U.S. Geological Survey.

Autor:

Antonio de Lucas Sepúlveda

Ingeniero de Caminos, Canales y Puertos (UPM) y Doctor por la Universidad de Alcalá en el programa Hidrología y Gestión de los Recursos Hídricos.
Ver todas las entradas de Antonio de Lucas Sepúlveda
Marcar el enlace permanente.

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *