0.1 Carga de datos

setwd("/cloud/project")
datos <- read.csv("derrames_globales_.csv", header = TRUE, sep = ";" , dec = ".")
str(datos)
## 'data.frame':    3550 obs. of  23 variables:
##  $ Id                                 : int  6786 6250 8220 6241 6216 6620 6262 6229 6201 6221 ...
##  $ Dia                                : int  19 3 21 16 19 7 10 12 18 29 ...
##  $ Mes                                : int  1 6 4 3 12 10 2 5 3 1 ...
##  $ ANo                                : chr  "A" "1979" "2010" "1978" ...
##  $ Nombre                             : chr  "Arabian Gulf Spills; Persian Gulf, Kuwait" "IXTOC I; Bahia de Campeche, Mexico" "Deepwater Horizon; Gulf of Mexico" "Amoco Cadiz; Brittany, France" ...
##  $ Ubicacion                          : chr  "Persian Gulf, Kuwait" "Bahia de Campeche, Mexico" "Gulf of Mexico" "Brittany, France" ...
##  $ Latitud                            : chr  "29,5" "19,4083" "28,7367" "48,5833" ...
##  $ Longuitud                          : chr  "48" "-92,325" "-88,3872" "-4,71667" ...
##  $ Amenaza                            : chr  "Oil" "Oil" "Oil" "Oil" ...
##  $ Etiquetas                          : chr  "" "Collision" "" "Grounding" ...
##  $ Tipo_de_crudo                      : chr  "Kuwait crude oil" "IXTOC I crude oil" "Diesel, crude oil" "Arabian light crude, Iranian light crude, Bunker C" ...
##  $ Recuperacion_en_superficie_aplicada: int  NA NA 1 NA NA NA NA NA NA NA ...
##  $ Recuperacion_en_costas_aplicada    : int  NA NA 1 NA NA NA NA NA NA NA ...
##  $ Tratamiento_biologico_aplicado     : int  1 NA 1 1 NA NA NA NA NA NA ...
##  $ Dispersion_quimica_aplicada        : int  NA 1 1 1 NA NA NA 1 1 1 ...
##  $ Quema_aplicada                     : int  NA 1 1 NA NA NA NA 1 1 1 ...
##  $ Maximo_liberacion_galones          : int  336000009 -1 205000000 68000017 -1 -1 -1 9240000 36100000 -1 ...
##  $ Barreras_de_contencion_flotantes   : int  35 12 182 17 3 3 7 8 5 6 ...
##  $ Causa_principal                    : chr  "DaNo del tanque  " "Incendio y explosion " "Incendio y explosion " "DaNo del tanque " ...
##  $ Volumen_derramados_galones         : chr  "336.000.000" "365.000.000" "600.000.000" "68.000.000" ...
##  $ Respuesta_actual_galones           : chr  "336000000" "252000000" "168000000" "68700000" ...
##  $ Fuente_respuesta                   : chr  "description and posts" "posts" "description" "posts" ...
##  $ etiqueta_actualizacion             : chr  "RA updated" "RA newly acquired" "RA updated" "RA updated" ...

0.2 Cargar librerías

if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
if (!require("fitdistrplus")) install.packages("fitdistrplus")
## Loading required package: fitdistrplus
## Loading required package: MASS
## Loading required package: survival
if (!require("MASS")) install.packages("MASS")
if (!require("knitr")) install.packages("knitr")
## Loading required package: knitr
library(ggplot2)
library(fitdistrplus)
library(MASS)
library(knitr)

0.3 Extraer variables

datos <- read.csv("derrames_globales_.csv",
                  header = TRUE, sep = ";", dec = ".", stringsAsFactors = FALSE)

datos$Latitud <- as.numeric(gsub(",", ".", datos$Latitud))
datos$Longuitud <- as.numeric(gsub(",", ".", datos$Longuitud))
datos$Maximo_liberacion_galones <- as.numeric(gsub(",", ".", datos$Maximo_liberacion_galones))
datos$Respuesta_actual_galones <- as.numeric(gsub(",", ".", datos$Respuesta_actual_galones))

0.4 Función auxiliar

generar_tabla_inferencia <- function(nombre_var, correlacion, chi2, umbral) {
  data.frame(
    Variable = nombre_var,
    `Test Pearson (%)` = round(correlacion, 2),
    `Chi Cuadrado` = round(chi2, 2),
    `Umbral de aceptación` = round(umbral, 2),
    Decisión = ifelse(chi2 < umbral, "Se acepta H0", "Se rechaza H0")
  )
}

1 Latitud

1.1 Modelo normal

latitud <- na.omit(datos$Latitud)
n_lat <- length(latitud)
mu_lat <- mean(latitud)
sd_lat <- sd(latitud)

hist_lat <- hist(latitud, breaks = nclass.Sturges(latitud), plot = FALSE)

plot(hist_lat, freq = FALSE, col = "lightblue",
     main = "Distribución Normal – Latitud",
     xlab = "Latitud", ylab = "Densidad")
curve(dnorm(x, mu_lat, sd_lat), add = TRUE, col = "red", lwd = 3)

Fo <- hist_lat$counts
breaks <- hist_lat$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  pnorm(breaks[i+1], mu_lat, sd_lat) -
    pnorm(breaks[i], mu_lat, sd_lat))

Fe <- P * n_lat
Fe[Fe < 1] <- 1

cor_lat <- cor(Fo, Fe) * 100
x2_lat <- sum((Fe - Fo)^2 / Fe)
gl_lat <- h - 1 - 2
umbral_lat <- qchisq(0.95, df = max(1, gl_lat))

print(kable(generar_tabla_inferencia("Latitud", cor_lat, x2_lat, umbral_lat),
            caption = "Test de bondad de ajuste – Latitud"))
## 
## 
## Table: Test de bondad de ajuste – Latitud
## 
## |Variable | Test.Pearson....| Chi.Cuadrado| Umbral.de.aceptación|Decisión      |
## |:--------|----------------:|------------:|--------------------:|:-------------|
## |Latitud  |            91.58|       670.87|                22.36|Se rechaza H0 |

1.2 FO vs FE

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Latitud",
     xlab = "Frecuencia Observada",
     ylab = "Frecuencia Esperada")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

1.3 Probabilidad en intervalo

li_lat <- mu_lat - sd_lat
ls_lat <- mu_lat + sd_lat
prob_lat <- pnorm(ls_lat, mu_lat, sd_lat) -
  pnorm(li_lat, mu_lat, sd_lat)

1.4 Gráfica probabilidad

x <- seq(min(latitud), max(latitud), length.out = 1000)
plot(x, dnorm(x, mu_lat, sd_lat), type = "l", col = "darkgreen", lwd = 2,
     main = "Probabilidad en intervalo – Latitud",
     xlab = "Latitud", ylab = "Densidad")

x_fill <- seq(li_lat, ls_lat, length.out = 500)
polygon(c(x_fill, rev(x_fill)),
        c(dnorm(x_fill, mu_lat, sd_lat), rep(0, length(x_fill))),
        col = rgb(1,0,0,0.4))

1.5 Intervalo de confianza

e_lat <- sd_lat / sqrt(n_lat)
IC_lat <- c(mu_lat - 1.96 * e_lat, mu_lat + 1.96 * e_lat)

cat("\n--- DATOS LATITUD ---\n")
## 
## --- DATOS LATITUD ---
cat("Modelo: Normal\n")
## Modelo: Normal
cat("Media:", round(mu_lat, 2), "\n")
## Media: 36.32
cat("Desviación estándar:", round(sd_lat, 2), "\n")
## Desviación estándar: 12.28
cat("Intervalo de confianza 95%: [",
    round(IC_lat[1], 2), ",",
    round(IC_lat[2], 2), "]\n")
## Intervalo de confianza 95%: [ 35.92 , 36.72 ]
cat("Probabilidad entre (μ ± σ):",
    round(prob_lat * 100, 2), "%\n")
## Probabilidad entre (μ ± σ): 68.27 %

1.6 Conclusión

La variable Latitud se explica adecuadamente mediante un modelo de probabilidad de distribución normal, con una media de 36.32 y una desviación estándar de 12.28. El intervalo de confianza del 95 % para la media poblacional se encuentra entre 35.92 y 36.72, lo que permite afirmar con un 95 % de confianza que la latitud promedio de los derrames petroleros se ubica dentro de dicho rango. Asimismo, la probabilidad de que la variable tome valores dentro del intervalo (μ ± σ) es del 68.27 %, lo cual es consistente con las propiedades teóricas de la distribución normal.

2 Longitud

2.1 Modelo Log-Normal

longitud <- na.omit(datos$Longuitud)
longitud_pos <- longitud - min(longitud) + 1

fit_long <- fitdist(longitud_pos, "lnorm")
ml <- fit_long$estimate["meanlog"]
sl <- fit_long$estimate["sdlog"]

hist_long <- hist(longitud_pos, breaks = nclass.Sturges(longitud_pos), plot = FALSE)

plot(hist_long, freq = FALSE, col = "thistle",
     main = "Distribución Log-normal – Longitud",
     xlab = "Longitud ajustada")
curve(dlnorm(x, ml, sl), add = TRUE, col = "darkgreen", lwd = 3)

Fo <- hist_long$counts
breaks <- hist_long$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  plnorm(breaks[i+1], ml, sl) -
    plnorm(breaks[i], ml, sl))

Fe <- P * length(longitud_pos)

2.2 FO vs FE

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Longitud",
     xlab = "FO", ylab = "FE")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

2.3 Probabilidad en intervalo

li_long <- quantile(longitud_pos, 0.25)
ls_long <- quantile(longitud_pos, 0.75)

prob_long <- plnorm(ls_long, ml, sl) -
  plnorm(li_long, ml, sl)

2.4 Gráfica probabilidad

x <- seq(min(longitud_pos), max(longitud_pos), length.out = 1000)
plot(x, dlnorm(x, ml, sl), type = "l", col = "darkgreen", lwd = 2,
     main = "Probabilidad en intervalo – Longitud",
     xlab = "Longitud ajustada", ylab = "Densidad")

x_fill <- seq(li_long, ls_long, length.out = 500)
polygon(c(x_fill, rev(x_fill)),
        c(dlnorm(x_fill, ml, sl), rep(0, length(x_fill))),
        col = rgb(1,0,0,0.4))

2.5 Media y desviación en escala original

media_long <- exp(ml + 0.5 * sl^2)
sd_long <- sqrt((exp(sl^2) - 1) * exp(2 * ml + sl^2))

cat("\n--- DATOS LONGITUD ---\n")
## 
## --- DATOS LONGITUD ---
cat("Modelo: Log-normal\n")
## Modelo: Log-normal
cat("Media:", round(media_long, 2), "\n")
## Media: 103.59
cat("Desviación estándar:", round(sd_long, 2), "\n")
## Desviación estándar: 44.77
cat("Intervalo de probabilidad (Q1–Q3):",
    round(prob_long * 100, 2), "%\n")
## Intervalo de probabilidad (Q1–Q3): 42.56 %

2.6 Conclusión

La variable Longitud presenta un comportamiento asimétrico positivo, por lo que se modela mediante una distribución log-normal. En la escala original, se obtiene una media de 103.59 y una desviación estándar de 44.77. La probabilidad de que la longitud se encuentre entre el primer y tercer cuartil (Q1–Q3) es del 42.56 %, lo que indica una concentración moderada de observaciones en el rango central de la distribución. Este modelo resulta adecuado para representar la variabilidad espacial longitudinal de los derrames analizados.

3 Maxima Liberación

max_lib <- na.omit(datos$Maximo_liberacion_galones)
max_lib <- max_lib[max_lib > 0]

3.1 Rango 1: Gamma (≤ 40.000)

r1 <- max_lib[max_lib <= 40000]
fit_g <- fitdist(r1, "gamma", method = "mme")
shape_g <- fit_g$estimate["shape"]
rate_g <- fit_g$estimate["rate"]

hist_r1 <- hist(r1, breaks = nclass.Sturges(r1), plot = FALSE)

plot(hist_r1, freq = FALSE, col = "orange",
     main = "Max Liberación – Rango 1 (Gamma)",
     xlab = "Galones")
curve(dgamma(x, shape_g, rate_g), add = TRUE, col = "blue", lwd = 3)

Fo <- hist_r1$counts
breaks <- hist_r1$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  pgamma(breaks[i+1], shape_g, rate_g) -
    pgamma(breaks[i], shape_g, rate_g))

Fe <- P * length(r1)

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Max Liberación (Rango 1)",
     xlab = "FO", ylab = "FE")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

3.2 Probabilidad Rango 1

li <- 5000
ls <- 20000
prob_r1 <- pgamma(ls, shape_g, rate_g) -
  pgamma(li, shape_g, rate_g)


media_r1 <- shape_g / rate_g
sd_r1 <- sqrt(shape_g) / rate_g

cat("\n--- DATOS MAX LIBERACIÓN (RANGO 1) ---\n")
## 
## --- DATOS MAX LIBERACIÓN (RANGO 1) ---
cat("Modelo: Gamma\n")
## Modelo: Gamma
cat("Media:", round(media_r1, 2), "\n")
## Media: 5041.42
cat("Desviación estándar:", round(sd_r1, 2), "\n")
## Desviación estándar: 8060.4
cat("Probabilidad entre 5000 y 20000 galones:",
    round(prob_r1 * 100, 2), "%\n")
## Probabilidad entre 5000 y 20000 galones: 24.25 %

3.3 Conclusión - Rango 1

Para valores de máxima liberación de hasta 40 000 galones, la variable se ajusta a un modelo Gamma, con una media de 5 041.42 galones y una desviación estándar de 8 060.40 galones. La probabilidad de que el volumen liberado se encuentre entre 5 000 y 20 000 galones es del 24.25 %, lo que evidencia una alta dispersión en los volúmenes de derrame dentro de este rango. Este comportamiento refleja la variabilidad inherente a eventos de liberación moderada de hidrocarburos.

3.4 Rango 2: Modelo Log-Normal (> 40.000)

r2 <- max_lib[max_lib > 40000]
fit_ln <- fitdist(r2, "lnorm")
ml_r2 <- fit_ln$estimate["meanlog"]
sl_r2 <- fit_ln$estimate["sdlog"]

hist_r2 <- hist(r2, breaks = nclass.Sturges(r2), plot = FALSE)

plot(hist_r2, freq = FALSE, col = "steelblue",
     main = "Max Liberación – Rango 2 (Log-normal)",
     xlab = "Galones")
curve(dlnorm(x, ml_r2, sl_r2), add = TRUE, col = "red", lwd = 3)

Fo <- hist_r2$counts
breaks <- hist_r2$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  plnorm(breaks[i+1], ml_r2, sl_r2) -
    plnorm(breaks[i], ml_r2, sl_r2))

Fe <- P * length(r2)

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Max Liberación (Rango 2)",
     xlab = "FO", ylab = "FE")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

3.5 Probabilidad Rango 2

li <- 50000
ls <- 100000
prob_r2 <- plnorm(ls, ml_r2, sl_r2) -
  plnorm(li, ml_r2, sl_r2)

media_r2 <- exp(ml_r2 + 0.5 * sl_r2^2)
sd_r2 <- sqrt((exp(sl_r2^2) - 1) * exp(2 * ml_r2 + sl_r2^2))

cat("\n--- DATOS MAX LIBERACIÓN (RANGO 2) ---\n")
## 
## --- DATOS MAX LIBERACIÓN (RANGO 2) ---
cat("Modelo: Log-normal\n")
## Modelo: Log-normal
cat("Media:", round(media_r2, 2), "\n")
## Media: 2162817
cat("Desviación estándar:", round(sd_r2, 2), "\n")
## Desviación estándar: 10786748
cat("Probabilidad entre 50000 y 100000 galones:",
    round(prob_r2 * 100, 2), "%\n")
## Probabilidad entre 50000 y 100000 galones: 9.35 %

3.6 Conclusión - Rango 2

Para valores superiores a 40 000 galones, la variable se modela mediante una distribución log-normal, con una media de 2 162 817 galones y una desviación estándar de 10 786 748 galones, lo que indica una dispersión extremadamente alta. La probabilidad de que el volumen liberado se sitúe entre 50 000 y 100 000 galones es únicamente del 9.35 %, evidenciando la presencia de eventos extremos de gran magnitud que influyen significativamente en la media.

4 Respuesta Actual

resp <- na.omit(datos$Respuesta_actual_galones)
resp <- resp[resp > 0]

4.1 Rango 1: Gamma (≤ 10.000)

rr1 <- resp[resp <= 10000]
fit_gr <- fitdist(rr1, "gamma", method = "mme")
shape_gr <- fit_gr$estimate["shape"]
rate_gr <- fit_gr$estimate["rate"]

hist_rr1 <- hist(rr1, breaks = nclass.Sturges(rr1), plot = FALSE)

plot(hist_rr1, freq = FALSE, col = "lightgreen",
     main = "Respuesta Actual – Rango 1 (Gamma)",
     xlab = "Galones")
curve(dgamma(x, shape_gr, rate_gr), add = TRUE, col = "darkgreen", lwd = 3)

Fo <- hist_rr1$counts
breaks <- hist_rr1$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  pgamma(breaks[i+1], shape_gr, rate_gr) -
    pgamma(breaks[i], shape_gr, rate_gr))

Fe <- P * length(rr1)

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Respuesta Actual (Rango 1)",
     xlab = "FO", ylab = "FE")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

4.2 Probabilidad Rango 1

li <- 2000
ls <- 8000
prob_rr1 <- pgamma(ls, shape_gr, rate_gr) -
  pgamma(li, shape_gr, rate_gr)

media_rr1 <- shape_gr / rate_gr
sd_rr1 <- sqrt(shape_gr) / rate_gr

cat("\n--- DATOS RESPUESTA ACTUAL (RANGO 1) ---\n")
## 
## --- DATOS RESPUESTA ACTUAL (RANGO 1) ---
cat("Modelo: Gamma\n")
## Modelo: Gamma
cat("Media:", round(media_rr1, 2), "\n")
## Media: 2009.5
cat("Desviación estándar:", round(sd_rr1, 2), "\n")
## Desviación estándar: 2510.08
cat("Probabilidad entre 2000 y 8000 galones:",
    round(prob_rr1 * 100, 2), "%\n")
## Probabilidad entre 2000 y 8000 galones: 30.24 %

4.3 Conclusión - Rango 1

Para valores de respuesta actual hasta 10 000 galones, la variable sigue un modelo Gamma, con una media de 2 009.50 galones y una desviación estándar de 2 510.08 galones. La probabilidad de que la respuesta se encuentre entre 2 000 y 8 000 galones es del 30.24 %, lo que indica que una proporción importante de los eventos presenta respuestas operativas de magnitud moderada. ## Rango 2: Log-Normal (> 10.000)

rr2 <- resp[resp > 10000]
fit_ln_r <- fitdist(rr2, "lnorm")
ml_rr2 <- fit_ln_r$estimate["meanlog"]
sl_rr2 <- fit_ln_r$estimate["sdlog"]

hist_rr2 <- hist(rr2, breaks = nclass.Sturges(rr2), plot = FALSE)

plot(hist_rr2, freq = FALSE, col = "gray",
     main = "Respuesta Actual – Rango 2 (Log-normal)",
     xlab = "Galones")
curve(dlnorm(x, ml_rr2, sl_rr2), add = TRUE, col = "red", lwd = 3)

Fo <- hist_rr2$counts
breaks <- hist_rr2$breaks
h <- length(Fo)

P <- sapply(1:h, function(i)
  plnorm(breaks[i+1], ml_rr2, sl_rr2) -
    plnorm(breaks[i], ml_rr2, sl_rr2))

Fe <- P * length(rr2)

plot(Fo, Fe, pch = 19, col = "blue3",
     main = "FO vs FE – Respuesta Actual (Rango 2)",
     xlab = "FO", ylab = "FE")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)

4.4 Probabilidad Rango 2

li <- 20000
ls <- 60000
prob_rr2 <- plnorm(ls, ml_rr2, sl_rr2) -
  plnorm(li, ml_rr2, sl_rr2)

media_rr2 <- exp(ml_rr2 + 0.5 * sl_rr2^2)
sd_rr2 <- sqrt((exp(sl_rr2^2) - 1) * exp(2 * ml_rr2 + sl_rr2^2))

cat("\n--- DATOS RESPUESTA ACTUAL (RANGO 2) ---\n")
## 
## --- DATOS RESPUESTA ACTUAL (RANGO 2) ---
cat("Modelo: Log-normal\n")
## Modelo: Log-normal
cat("Media:", round(media_rr2, 2), "\n")
## Media: 1075465
cat("Desviación estándar:", round(sd_rr2, 2), "\n")
## Desviación estándar: 8669338
cat("Probabilidad entre 20000 y 60000 galones:",
    round(prob_rr2 * 100, 2), "%\n")
## Probabilidad entre 20000 y 60000 galones: 17.16 %

4.5 Conclusión - Rango 2

Para valores superiores a 10 000 galones, la respuesta actual se ajusta a una distribución log-normal, con una media de 1 075 465 galones y una desviación estándar de 8 669 338 galones. La probabilidad de que la respuesta se sitúe entre 20 000 y 60 000 galones es del 17.16 %, lo que evidencia una alta variabilidad y la presencia de respuestas excepcionales asociadas a derrames de gran magnitud.