Variable Cuantitativa Continua: Pérdida neta de barriles

Variable Original: Net.Loss..Barrels.

Justificación de la variable

El análisis de la variable Pérdida neta de barriles se justifica porque cuantifica el volumen real de producto que no pudo ser recuperado tras un siniestro, representando el impacto ambiental definitivo y la pérdida económica neta del incidente. Su estudio estadístico permite evaluar el nivel de eficiencia en la respuesta de emergencia, las estrategias de contención y limpieza de derrames, y generar información clave para mitigar las consecuencias a largo plazo.

1 Cargar datos

En esta fase, el análisis se expande hacia variables numéricas específicas de impacto real, como el volumen de barriles perdidos definitivamente. La función read.csv() importa la estructura completa del archivo database-_1_.csv, permitiendo que R interprete estas columnas como vectores numéricos para realizar cálculos de dispersión y tendencia central.

datos <- read.csv("database-_1_.csv")
zona <- datos$Net.Loss..Barrels.

library(gt)
library(dplyr)
library(ggplot2)
library(dplyr)
library(scales)

2 Tabla de frecuencia

La tabla de frecuencia aplicada a la variable de pérdida permite segmentar la magnitud del producto no recuperado. Al organizar los datos mediante la Regla de Sturges, se facilita la identificación del rango con mayor concentración de eventos, el cual suele estar fuertemente sesgado hacia valores bajos (cero o cercanos a cero) en incidentes donde la recuperación fue exitosa.

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels. 
costos <- na.omit(variable_interes)

k <- 1 + (3.322 * log10(length(costos)))
k <- floor(k)
min_val <- min(costos)
max_val <- max(costos)
R_val <- max_val - min_val
A <- R_val / k

Li_num <- seq(from = min_val, to = max_val - A, by = A)
if(length(Li_num) < k) { Li_num <- c(Li_num, Li_num[length(Li_num)] + A) }
if(max(Li_num) + A < max_val) { Li_num <- c(Li_num, tail(Li_num, 1) + A) }

Ls_num <- Li_num + A
MC_num <- (Li_num + Ls_num) / 2

ni <- numeric(length(Li_num))
for (i in 1:length(Li_num)) {
  if (i == length(Li_num)) {
      ni[i] <- sum(costos >= Li_num[i] & costos <= (max_val + 100000))
  } else {
      ni[i] <- sum(costos >= Li_num[i] & costos < Ls_num[i]) 
  }
}
num_ceros <- sum(ni == 0)
if (num_ceros > 0) {
  ni[ni == 0] <- 1
  idx_max <- which.max(ni)
  ni[idx_max] <- ni[idx_max] - num_ceros
}

hi <- ni / sum(ni) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc <- round(cumsum(hi), 2)
Hidsc <- round(rev(cumsum(rev(hi))), 2)

TDFCostos <- data.frame(
  Li_num = Li_num, 
  Ls_num = Ls_num, 
  MC_num = MC_num, 
  ni = ni, 
  hi = hi, 
  Niasc = Niasc, 
  Nidsc = Nidsc, 
  Hiasc = Hiasc, 
  Hidsc = Hidsc
)

tabla1_sturges <- TDFCostos %>%
  gt() %>%
  tab_header(
    title = md("*Tabla 1: Distribución de Frecuencias (Pérdida Neta)*"),
  ) %>%
  cols_label(
    Li_num = "Desde (bbl)",
    Ls_num = "Hasta (bbl)",
    MC_num = "Marca Clase",
    ni = "Frec. Abs.",
    hi = "Frec. Rel. %",
    Niasc = "Ni Asc.",
    Nidsc = "Ni Desc.",
    Hiasc = "Hi Asc. %",
    Hidsc = "Hi Desc. %"
  ) %>%
  fmt_number(columns = c(Li_num, Ls_num, MC_num), decimals = 2) %>%
  fmt_number(columns = c(hi, Hiasc, Hidsc), decimals = 2, pattern = "{x}%")

tabla1_sturges

Tabla 1: Distribución de Frecuencias (Pérdida Neta)
Desde (bbl)	Hasta (bbl)	Marca Clase	Frec. Abs.	Frec. Rel. %	Ni Asc.	Ni Desc.	Hi Asc. %	Hi Desc. %
0.00	2,547.08	1,273.54	2760	98.75%	2760	2795	98.75%	100.00%
2,547.08	5,094.17	3,820.62	18	0.64%	2778	35	99.39%	1.25%
5,094.17	7,641.25	6,367.71	5	0.18%	2783	17	99.57%	0.61%
7,641.25	10,188.33	8,914.79	2	0.07%	2785	12	99.64%	0.43%
10,188.33	12,735.42	11,461.88	1	0.04%	2786	10	99.68%	0.36%
12,735.42	15,282.50	14,008.96	3	0.11%	2789	9	99.79%	0.32%
15,282.50	17,829.58	16,556.04	1	0.04%	2790	6	99.82%	0.21%
17,829.58	20,376.67	19,103.12	1	0.04%	2791	5	99.86%	0.18%
20,376.67	22,923.75	21,650.21	1	0.04%	2792	4	99.89%	0.14%
22,923.75	25,470.83	24,197.29	1	0.04%	2793	3	99.93%	0.11%
25,470.83	28,017.92	26,744.38	1	0.04%	2794	2	99.96%	0.07%
28,017.92	30,565.00	29,291.46	1	0.04%	2795	1	100.00%	0.04%

3 Histograma Global de Cantidad Absoluta

En esta sección se analiza la magnitud física definitiva de los incidentes mediante el conteo de la pérdida neta. Al ser una variable continua con una asimetría muy acentuada, la Regla de Sturges permite determinar el número óptimo de intervalos (\(k\)). Aquí evaluamos la concentración de la siniestralidad dentro del primer intervalo de clase, donde se aglomera la inmensa mayoría de los eventos (generalmente aquellos con recuperación total).

variable_interes <- datos$Net.Loss..Barrels.
volumen <- na.omit(variable_interes)
k <- 1 + (3.322 * log10(length(volumen)))
R <- max(volumen) - min(volumen)
A <- R / floor(k) 

limit_zoom <- A 
datos_zoom <- datos %>%
  filter(!is.na(Net.Loss..Barrels.)) %>%
  filter(Net.Loss..Barrels. <= limit_zoom)

p_zoom <- ggplot(datos_zoom, aes(x = Net.Loss..Barrels.)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white", alpha = 0.8) +
  scale_x_continuous(labels = scales::comma_format(suffix = " bbl")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  
  labs(
    title = "Gráfica 1: Distribución Global de Pérdida Neta", 
    x = "Volumen (bbl)",
    y = "Cantidad"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(p_zoom)

4 Histograma local de cantidad absoluta del Rango Principal

Al igual que con la liberación inicial, se seleccionó el intervalo de 0 a 2.0 barriles para un análisis más fino. Dado que la distribución global de la Pérdida Neta es aún más asimétrica (ya que muchos derrames son completamente limpiados, resultando en 0 barriles perdidos), este enfoque permite “desempaquetar” la primera barra del histograma general. Esto revela cómo se comportan aquellos incidentes que no lograron una recuperación perfecta pero mantuvieron una pérdida final baja.

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels.
volumen <- na.omit(variable_interes)
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Net.Loss..Barrels.)) %>%
  filter(Net.Loss..Barrels. <= limit_zoom)

p_zoom_5barras <- ggplot(datos_zoom, aes(x = Net.Loss..Barrels.)) +
  geom_histogram(bins = 5, fill = "steelblue", color = "white", alpha = 0.8) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1)) +
  
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  
  labs(
    title = "Gráfica No2: Distribución de pérdida en el rango de mayor ocurrencia", 
    x = "Volumen Pérdida Neta (bbl)",
    y = "Cantidad"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras)

5 Histograma local de cantidad relativa del rango principal

Estas secciones transforman las frecuencias absolutas en porcentajes para entender el peso de cada intervalo dentro de los derrames más comunes. Enfocarse en el “Rango Principal” ayuda a dimensionar en qué proporción los incidentes se logran contener antes de que la pérdida neta sobrepase la fracción de un barril.

library(ggplot2)
library(dplyr)
library(scales)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels.
volumen <- na.omit(variable_interes)
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Net.Loss..Barrels.)) %>%
  filter(Net.Loss..Barrels. <= limit_zoom)

p_zoom_5barras_pct <- ggplot(datos_zoom, aes(x = Net.Loss..Barrels.)) +
  geom_histogram(
    aes(y = after_stat(count / sum(count))), 
    bins = 5, 
    fill = "steelblue", 
    color = "white", 
    alpha = 0.8
  ) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1, suffix = " bbl")) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1), 
    expand = expansion(mult = c(0, 0.05))
  ) +
  
  labs(
    title = "Gráfica No3 : Distribución porcentual local del rango principal", 
    x = "Volumen Pérdida Neta (bbl)",
    y = "Porcentaje (%)"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras_pct)

6 Histograma global de cantidad relativa del rango principal

A través de las librerías scales y ggplot2, se visualiza la probabilidad de que una pérdida neta caiga dentro de un rango específico en proporción a todo el subgrupo. Como la mayoría de los derrames alcanzan una recuperación del 100%, el porcentaje asociado a la pérdida exacta de cero dominará la visualización frente al resto de la distribución.

library(ggplot2)
library(dplyr)
library(scales)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels.
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Net.Loss..Barrels.)) %>%
  filter(Net.Loss..Barrels. <= limit_zoom)

p_zoom_5barras_100 <- ggplot(datos_zoom, aes(x = Net.Loss..Barrels.)) +
  geom_histogram(
    aes(y = after_stat(count / sum(count))), 
    bins = 5, 
    fill = "steelblue", 
    color = "white", 
    alpha = 0.8
  ) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1, suffix = " bbl")) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    limits = c(0, 1), 
    expand = expansion(mult = c(0, 0.05))
  ) +
  
  labs(
    title = "Gráfica No4 : Distribución porcentual global del rango principal", 
    x = "Volumen Pérdida Neta (bbl)",
    y = "Porcentaje (%)"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras_100)

7 Ojivas combinadas de la cantidad absoluta

La ojiva combinada representa la suma acumulativa de los barriles no recuperados en los diferentes rangos. Nos ayuda a localizar gráficamente el punto medio de la pérdida neta (la mediana), que para esta variable estará extremadamente desplazada hacia la izquierda debido a la frecuencia de ceros absolutos (éxito total de recuperación).

library(ggplot2)
library(dplyr)
library(scales)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels.
volumen <- na.omit(variable_interes)
datos_local <- volumen[volumen <= 2]
k_local <- 6
min_val <- 0
max_val <- 2
R_local <- max_val - min_val
A_local <- R_local / k_local 
Li_num <- seq(min_val, max_val - A_local, length.out = k_local)
Ls_num <- Li_num + A_local
ni_local <- numeric(k_local)

for(i in 1:k_local){
  if(i == k_local){
    ni_local[i] <- sum(datos_local >= Li_num[i] & datos_local <= max_val)
  } else {
    ni_local[i] <- sum(datos_local >= Li_num[i] & datos_local < Ls_num[i])
  }
}

Niasc <- cumsum(ni_local)
Nidsc <- rev(cumsum(rev(ni_local)))

datos_asc <- data.frame(
  x = c(min_val, Ls_num),  
  y = c(0, Niasc),         
  Tipo = "Ascendente"
)

datos_dsc <- data.frame(
  x = c(Li_num, max_val),  
  y = c(Nidsc, 0),         
  Tipo = "Descendente"
)

datos_ojivas_plot <- rbind(datos_asc, datos_dsc)
p_ojiva_cruzada_solida <- ggplot(datos_ojivas_plot, aes(x = x, y = y, color = Tipo, linetype = Tipo)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2) +
  
  scale_x_continuous(
    labels = scales::number_format(accuracy = 0.01, suffix = " bbl"),
    breaks = scales::pretty_breaks(n = 6)
  ) +
  scale_color_manual(values = c("Ascendente" = "black", "Descendente" = "blue")) +
  scale_linetype_manual(values = c("Ascendente" = "solid", "Descendente" = "solid")) +
  
  labs(
    title = "Gráfica 5: Distribución Acumulada Comparativa del rango principal",
    x = "Volumen Pérdida Neta (bbl)",
    y = "Cantidad Acumulada",
    color = NULL,
    linetype = NULL
  ) +
  
  theme_bw() + 
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    legend.position = c(0.85, 0.5),
    legend.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black")
  )

print(p_ojiva_cruzada_solida)

# Esta gráfica presenta las Ojivas de Frecuencia Acumulada para la Pérdida Neta, mostrando cómo se distribuyen las faltantes de producto en el rango de 0 a 2 bbl. La línea negra indica cuántas observaciones tienen una pérdida menor que cierto valor, mientras que la azul muestra cuántos incidentes tienen una pérdida superior. El punto de intersección permite aproximar gráficamente la mediana.

8 Diagrama de caja

El boxplot ayuda a observar la concentración del 50% central de la pérdida neta. En este caso particular, demostrará que la gran mayoría de la pérdida neta se condensa en un rango estrechísimo cercano a 0 barriles, resaltando la asimetría de forma clara.

library(ggplot2)
library(dplyr)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Net.Loss..Barrels.
volumen <- na.omit(variable_interes)
variable_box <- volumen[volumen <= 2] # Filtro Local
par(mar = c(5, 2, 4, 2)) 
b <- boxplot(variable_box, 
        horizontal = TRUE, 
        col = "skyblue", 
        border = "gray30",         
        medcol = "red",            
        boxwex = 0.6,              
        outline = FALSE,           
        main = "Gráfica 6: Distribución de Volumen (Pérdida Neta) del rango principal",
        xlab = "Volumen Pérdida Neta (bbl)",
        xaxt = "n",  
        yaxt = "n",  
        frame = FALSE 
)

limite_visible <- b$stats[5] 
puntos_eje <- pretty(c(0, limite_visible))             
axis(1, at = puntos_eje, labels = format(puntos_eje, nsmall = 2), col = "gray30", col.axis = "gray30")
grid(nx = NULL, ny = NA, col = "lightgray", lty = "dotted", lwd = 1)

# El boxplot confirma que en el rango principal (0-2 bbl), la pérdida neta tiene una altísima concentración en 0 (que empuja los cuartiles inferiores y la mediana al inicio de la escala), lo que refleja que en más de la mitad de los incidentes se logró una recuperación del 100%. La línea punteada se extenderá hacia los valores mayores, indicando el sesgo.

9 Tabla estadísticos local

Consolida la evidencia numérica del análisis en la sección filtrada (0 a 2 barriles). Aporta valores puntuales de media, mediana, moda y dispersión, revelando que para el incidente rutinario o de menor cuantía, la recuperación es altamente eficiente.

library(e1071)

## 
## Adjuntando el paquete: 'e1071'

## The following object is masked from 'package:ggplot2':
## 
##     element

library(knitr)

datos <- read.csv("database-_1_.csv")
variable_raw <- na.omit(datos$Net.Loss..Barrels.)
variable_analisis <- variable_raw[variable_raw <= 2]
n <- length(variable_analisis)
k_local <- 6
min_val <- 0
max_val <- 2
A_local <- (max_val - min_val) / k_local
Li <- seq(min_val, max_val - A_local, length.out = k_local)
Ls <- Li + A_local
MC <- (Li + Ls) / 2

ni <- numeric(k_local)
for(i in 1:k_local){
  if(i == k_local) ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis <= max_val)
  else ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis < Ls[i])
}

media_agrupada <- sum(MC * ni) / sum(ni)
desviacion_estandar <- sd(variable_analisis)
error_estandar <- desviacion_estandar / sqrt(n)
margen_error <- 1.96 * error_estandar
ic_inferior <- media_agrupada - margen_error
ic_superior <- media_agrupada + margen_error
texto_media_intervalo <- paste0("[", format(round(ic_inferior, 3), nsmall=3), " - ", format(round(ic_superior, 3), nsmall=3), "]")
ri <- min(variable_analisis)
rs <- max(variable_analisis)
mediana <- median(variable_analisis)
t <- table(variable_analisis)
Mo <- as.numeric(names(t)[which.max(t)])
cv <- (desviacion_estandar / media_agrupada) * 100
As <- skewness(variable_analisis)
K <- kurtosis(variable_analisis)

Tabla_local <- data.frame(
  "Pérdida Neta (Local 0-2)",
  paste(format(ri, nsmall=2), "bbl"),
  paste(format(rs, nsmall=2), "bbl"),
  texto_media_intervalo,  
  paste(format(round(mediana, 3), big.mark=","), "bbl"),
  paste(format(round(Mo, 3), big.mark=","), "bbl"),
  paste(format(round(desviacion_estandar, 3), big.mark=","), "bbl"),
  paste(round(cv, 2), "%"),
  round(As, 2),
  round(K, 2)
)

colnames(Tabla_local) <- c("Variable","Min","Max","Media (IC 95%)","Mediana","Moda","Desv. S","CV","As","K")
kable(Tabla_local, format = "markdown", caption = "Tabla 1: Indicadores Locales de Pérdida Neta (con Intervalo de Confianza para la Media).")

Tabla 1: Indicadores Locales de Pérdida Neta (con Intervalo de Confianza para la Media).

Variable	Min	Max	Media (IC 95%)	Mediana	Moda	Desv. S	CV	As	K
Pérdida Neta (Local 0-2)	0.00 bbl	2.00 bbl	[0.316 - 0.356]	0 bbl	0 bbl	0.469 bbl	139.6 %	2.58	6.05

10 Tabla estadísticos global

Ofrece el panorama macro de todos los datos en relación a la pérdida del producto. A diferencia del rango local, los estadísticos globales exponen cómo la media real y la dispersión se inflan drásticamente por la ocurrencia de aquellos eventos que fracasan en la mitigación.

library(e1071)
library(knitr)

datos <- read.csv("database-_1_.csv")
variable_analisis <- na.omit(datos$Net.Loss..Barrels.)
n <- length(variable_analisis)
k_global <- floor(1 + 3.322 * log10(n))
R <- max(variable_analisis) - min(variable_analisis)
A_global <- R / k_global
Li <- seq(min(variable_analisis), max(variable_analisis) - A_global, length.out = k_global)

if(max(Li) + A_global < max(variable_analisis)) { 
   Li <- c(Li, tail(Li, 1) + A_global) 
   k_global <- k_global + 1
}
Ls <- Li + A_global
MC <- (Li + Ls) / 2

ni <- numeric(length(MC))
for(i in 1:length(MC)){
  if(i == length(MC)) ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis <= (max(variable_analisis) + 0.001))
  else ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis < Ls[i])
}

media_agrupada <- sum(MC * ni) / sum(ni)
desviacion_estandar <- sd(variable_analisis)
error_estandar <- desviacion_estandar / sqrt(n)
margen_error <- 1.96 * error_estandar
ic_inferior <- media_agrupada - margen_error
ic_superior <- media_agrupada + margen_error
texto_media_intervalo <- paste0("[", format(round(ic_inferior, 2), big.mark=","), " - ", format(round(ic_superior, 2), big.mark=","), "]")
ri <- min(variable_analisis)
rs <- max(variable_analisis)
mediana <- median(variable_analisis)
t <- table(variable_analisis)
Mo <- as.numeric(names(t)[which.max(t)])
cv <- (desviacion_estandar / media_agrupada) * 100
As <- skewness(variable_analisis)
K <- kurtosis(variable_analisis)

Tabla_global <- data.frame(
  "Pérdida Neta (Global)",
  paste(format(ri, nsmall=2), "bbl"),
  paste(format(rs, big.mark=","), "bbl"),
  texto_media_intervalo, 
  paste(format(round(mediana, 2), big.mark=","), "bbl"),
  paste(format(round(Mo, 2), big.mark=","), "bbl"),
  paste(format(round(desviacion_estandar, 2), big.mark=","), "bbl"),
  paste(round(cv, 2), "%"),
  round(As, 2),
  round(K, 2)
)

colnames(Tabla_global) <- c("Variable","Min","Max","Media (IC 95%)","Mediana","Moda","Desv. S","CV","As","K")
kable(Tabla_global, format = "markdown", caption = "Tabla 2: Indicadores Globales de Pérdida Neta (con Intervalo de Confianza para la Media).")

Tabla 2: Indicadores Globales de Pérdida Neta (con Intervalo de Confianza para la Media).

Variable	Min	Max	Media (IC 95%)	Mediana	Moda	Desv. S	CV	As	K
Pérdida Neta (Global)	0.00 bbl	30,565 bbl	[1,311.63 - 1,399.49]	0 bbl	0 bbl	1,185.02 bbl	87.42 %	17.07	350.69

11 Valores Atípicos

Dada la naturaleza crítica de algunos siniestros, en ocasiones la recuperación falla estrepitosamente y grandes cantidades de producto se pierden de manera irremediable en el medio ambiente. Esta sección calcula y aísla estos eventos catastróficos que se ubican fuera del comportamiento típico de la distribución.

variable_global <- na.omit(datos$Net.Loss..Barrels.)
stats_outliers_global <- boxplot.stats(variable_global)$out
num_outliers_global <- length(stats_outliers_global)
minimooutliers_global <- if(num_outliers_global > 0) min(stats_outliers_global) else NA
maximooutliers_global <- if(num_outliers_global > 0) max(stats_outliers_global) else NA

cat("\n--- Análisis de Outliers de Pérdida Neta (GLOBAL) ---\n")

## 
## --- Análisis de Outliers de Pérdida Neta (GLOBAL) ---

cat("Número de valores atípicos:", num_outliers_global, "\n")

## Número de valores atípicos: 482

cat("Mínimo Outlier:", if(!is.na(minimooutliers_global)) paste(format(minimooutliers_global, big.mark=","), "bbl") else "Ninguno", "\n")

## Mínimo Outlier: 5.23 bbl

cat("Máximo Outlier:", if(!is.na(maximooutliers_global)) paste(format(maximooutliers_global, big.mark=","), "bbl") else "Ninguno", "\n")

## Máximo Outlier: 30,565 bbl

12 Conclusión

El análisis estadístico de la Pérdida Neta revela una distribución de asimetría positiva extrema que evidencia el principio de pareto en la siniestralidad operativa: los protocolos de contención primaria son altamente eficientes frente a eventos rutinarios (consolidando la Mediana y la Moda en cero), demostrando una gran resiliencia del sistema. Sin embargo, el arrastre de la Media hacia la derecha subraya una vulnerabilidad crítica ante valores atípicos severos (eventos catastróficos aislados o cisnes negros). En conclusión, el impacto ambiental real no está dictado por la frecuencia de los derrames, sino por el colapso de la contención en estos escenarios extremos, exigiendo que la gestión de riesgos evolucione de controlar el volumen diario a mitigar exclusivamente la criticidad atípica.