Clase de Hidrología Estadística: Distribuciones en R

Bienvenido. Vamos a replicar en R lo que hicieron en Python con scipy.stats. La buena noticia es que R tiene un sistema mucho más elegante y unificado para trabajar con distribuciones: el sistema de prefijos.

1. El sistema de prefijos en R

En R, cada distribución tiene 4 funciones asociadas, identificadas por una letra inicial. Es un patrón que memorizas una vez y aplica a todas las distribuciones (norm, gamma, pois, binom, lnorm, exp, weibull, etc.).

Prefijo Nombre Equivalente en scipy.stats Qué hace
d density / distribución pdf (continuas) o pmf (discretas) Devuelve la función de densidad (o masa) en un punto x
p probabilidad acumulada cdf Devuelve \(P(X \le x)\)
q quantil ppf Dado \(p \in (0,1)\), devuelve \(x\) tal que \(P(X\le x)=p\)
r random rvs Genera valores aleatorios de la distribución

Ejemplos con norm, gamma y pois

# ---------- NORMAL ----------
dnorm(x = 0, mean = 0, sd = 1)     # densidad en x = 0  → 0.3989
pnorm(q = 1.96, mean = 0, sd = 1)  # P(X ≤ 1.96)        → 0.975
qnorm(p = 0.95, mean = 0, sd = 1)  # cuantil 0.95       → 1.6449
rnorm(n = 10, mean = 0, sd = 1)    # 10 valores aleatorios

# ---------- GAMMA ----------
dgamma(x = 100, shape = 4, rate = 0.05)   # pdf en x = 100
pgamma(q = 100, shape = 4, rate = 0.05)   # P(X ≤ 100)
qgamma(p = 0.90, shape = 4, rate = 0.05)  # cuantil 0.90
rgamma(n = 1000, shape = 4, rate = 0.05)  # 1000 simulaciones

# ---------- POISSON ----------
dpois(x = 3, lambda = 2.5)      # P(X = 3)  → pmf (¡no pdf!)
ppois(q = 3, lambda = 2.5)      # P(X ≤ 3)  → cdf
qpois(p = 0.90, lambda = 2.5)   # cuantil 0.90
rpois(n = 100, lambda = 2.5)    # 100 conteos simulados

Ojo hidrólogo: en gamma puedes parametrizar por rate (tasa \(\beta\)) o por scale (= \(1/\beta\)). Es la fuente #1 de errores. Verifica siempre con la media: \(\mu = \text{shape}/\text{rate}\).

2. Código completo: histograma simulado + curva teórica Gamma con ggplot2::stat_function

Supongamos que estamos modelando precipitación mensual (mm) con una Gamma de forma \(k=4\) y tasa \(\beta=0.05\) (media = 80 mm).

# ============================================================
# Simulación Gamma + curva teórica con ggplot2
# ============================================================
library(ggplot2)

# --- 1. Parámetros de la distribución Gamma ---
shape_k <- 4
rate_b  <- 0.05          # media = 4/0.05 = 80 mm

# --- 2. Simulación (equivalente a .rvs en scipy) ---
set.seed(2024)                          # reproducibilidad
n_sim   <- 5000
sims    <- rgamma(n = n_sim, shape = shape_k, rate = rate_b)

# --- 3. Histograma + curva teórica ---
ggplot(data.frame(x = sims), aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 40, fill = "steelblue",
                 color = "white", alpha = 0.6) +
  stat_function(
    fun  = dgamma,                       # función a graficar
    args = list(shape = shape_k, rate = rate_b),
    color = "firebrick", linewidth = 1.2
  ) +
  labs(
    title = "Distribución Gamma simulada vs. curva teórica",
    subtitle = sprintf("shape = %.1f, rate = %.2f  (media = %.1f mm)",
                       shape_k, rate_b, shape_k / rate_b),
    x = "Precipitación mensual (mm)",
    y = "Densidad"
  ) +
  theme_minimal(base_size = 13)

¿Por qué funciona stat_function?

  • stat_function no necesita datos: evalúa la función fun sobre el rango del eje x y la dibuja.
  • El truco clave es aes(y = after_stat(density)) en el histograma: sin esto, el eje y sería conteo y la curva teórica (densidad) no coincidiría en escala.

Alternativa 100% base R (equivalente):

hist(sims, breaks = 40, freq = FALSE,
     col = "steelblue", border = "white",
     main = "Gamma: simulación vs teoría",
     xlab = "Precipitación (mm)")
curve(dgamma(x, shape = shape_k, rate = rate_b),
      add = TRUE, col = "firebrick", lwd = 2)

3. Probabilidad acumulada y cuantil

Con la misma Gamma (\(k=4\), \(\beta=0.05\)):

# --- P(X ≤ 100)  →  equivale a scipy.stats.gamma.cdf(100, ...) ---
p_leq_100 <- pgamma(q = 100, shape = shape_k, rate = rate_b)
p_leq_100
#> [1] 0.5665299

# Interpretación hidrológica:
# "Hay ~56.65% de probabilidad de que la precipitación mensual
#  no supere los 100 mm."

# --- Cuantil 0.90  →  equivale a scipy.stats.gamma.ppf(0.90, ...) ---
x_90 <- qgamma(p = 0.90, shape = shape_k, rate = rate_b)
x_90
#> [1] 120.6787

# Interpretación:
# "El evento de precipitación mensual asociado a un período de
#  retorno de 10 años (no excedencia 0.90) es ~120.7 mm."
#  Nota: T = 1 / (1 - F) = 1 / 0.10 = 10 años.

Verificación de consistencia (¡buena práctica!)

# Si pgamma es correcto, qgamma debe invertirlo:
qgamma(pgamma(100, shape_k, rate_b), shape_k, rate_b)
#> [1] 100

# Y el cuantil 0.90 evaluado en la cdf debe dar 0.90:
pgamma(x_90, shape_k, rate_b)
#> [1] 0.9

Esta propiedad de inversa mutua entre p* y q* es la que te permite, en hidrología, pasar de probabilidad → magnitud del evento (diseño) y de magnitud → probabilidad (riesgo).

4. Tabla resumen: Python ↔︎ R

Acción scipy.stats (Python) R
Densidad / masa gamma.pdf(x, a, scale) / poisson.pmf(k, mu) dgamma(x, shape, rate) / dpois(k, lambda)
Acumulada gamma.cdf(x, a, scale) pgamma(x, shape, rate)
Cuantil gamma.ppf(p, a, scale) qgamma(p, shape, rate)
Simulación gamma.rvs(a, scale, size=n) rgamma(n, shape, rate)

5. Ejercicio propuesto para casa

  1. Simula \(n=10000\) valores Poisson con \(\lambda = 3.5\) (número de tormentas por mes) y grafica con ggplot2 el diagrama de barras de frecuencias relativas frente a dpois.
  2. Calcula \(P(X \le 5)\) y el cuantil 0.95 con ppois y qpois.
  3. Repite para la Normal con \(\mu=80\), \(\sigma=20\) y compara con pnorm/qnorm.

Pista: para Poisson usa geom_col sobre una tabla de frecuencias, no geom_histogram (es discreta).