Variables Continuas
Cargar los datos y librerías
setwd("D:/Data")
datos <- read.csv("derrames_globales_.csv", header = TRUE, sep = ";", dec =".")
library(dplyr)## Warning: package 'dplyr' was built under R version 4.5.2##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(janitor)## Warning: package 'janitor' was built under R version 4.5.2##
## Adjuntando el paquete: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(e1071)## Warning: package 'e1071' was built under R version 4.5.2df <- clean_names(datos)
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" ...1 Variable Latitud
latitud_raw <- gsub(",", ".", df$latitud)
latitud <- suppressWarnings(as.numeric(latitud_raw))
latitud <- na.omit(latitud)
# Regla de Sturges
n <- length(latitud)
k <- ceiling(1 + 3.322 * log10(n))
R <- max(latitud) - min(latitud)
A <- R / k
# Límites y Marca de Clase (MC)
liminf <- seq(from = min(latitud), by = A, length.out = k)
limsup <- liminf + A
breaks_lat <- c(liminf, max(limsup) + 0.01)
MC <- (liminf + limsup) / 2
rango_lat <- cut(latitud, breaks = breaks_lat, right = FALSE, include.lowest = TRUE)
ni <- as.numeric(table(rango_lat))
hi_perc <- (ni / n) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Consolidada
Tabla_Lat_Final <- data.frame(
liminf = round(liminf, 2),
limsup = round(limsup, 2),
MC = round(MC, 2),
ni = ni,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_Lat_Final)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 -78.00 -66.49 -72.25 2 0.06 2 3550 0.06 100.00
## 2 -66.49 -54.99 -60.74 1 0.03 3 3548 0.08 99.94
## 3 -54.99 -43.48 -49.23 2 0.06 5 3547 0.14 99.92
## 4 -43.48 -31.97 -37.73 4 0.11 9 3545 0.25 99.86
## 5 -31.97 -20.47 -26.22 7 0.20 16 3541 0.45 99.75
## 6 -20.47 -8.96 -14.71 10 0.28 26 3534 0.73 99.55
## 7 -8.96 2.55 -3.21 4 0.11 30 3524 0.85 99.27
## 8 2.55 14.05 8.30 30 0.85 60 3520 1.69 99.15
## 9 14.05 25.56 19.81 240 6.76 300 3490 8.45 98.31
## 10 25.56 37.06 31.31 1626 45.80 1926 3250 54.25 91.55
## 11 37.06 48.57 42.82 1189 33.49 3115 1624 87.75 45.75
## 12 48.57 60.08 54.32 330 9.30 3445 435 97.04 12.25
## 13 60.08 71.58 65.83 105 2.96 3550 105 100.00 2.961.1 Gráfica 1: Distribución de la latitud (Sturges)
Se usa Sturges para definir la amplitud de cada intervalo pero se prefiere usar el segundo histograma en la presentación final porque resulta altamente beneficioso para una lectura rápida. Aunque el primer gráfico refleja la matemática exacta de la regla de Sturges, el eje horizontal queda saturado de decimales, lo que hace que la gráfica sea un poco confusa e incómoda de interpretar. En cambio, el segundo histograma es mucho más limpio; al usar el ajuste automático de R con números enteros y redondos, la visualización se vuelve más clara y directa para el lector, manteniendo intacta la forma real de los datos y el tamaño de las barras.
ymax_lat <- max(Tabla_Lat_Final$ni) * 1.2
hist(latitud, breaks = breaks_lat, freq = TRUE,
main = "Histograma de la Latitud de Derrames Petroleros",
xlab = "Latitud (grados)", ylab = "Frecuencia",
col = terrain.colors(length(breaks_lat) - 1), border = "black", las = 1,
xlim = c(min(breaks_lat), max(breaks_lat)),
ylim = c(0, ymax_lat),
xaxt = "n") ## Warning in plot.histogram(r, freq = freq1, col = col, border = border, angle =
## angle, : the AREAS in the plot are wrong -- rather use 'freq = FALSE'axis(side = 1, at = breaks_lat, labels = round(breaks_lat, 2), las = 2, cex.axis = 0.8)
1.2 Gráfica 2: Distribución de la latitud (R)
# Extraer los cálculos automáticos de R
h_auto <- hist(latitud, plot = FALSE)
# R ya nos da los breaks limpios y el conteo exacto
breaks_r <- h_auto$breaks
ni_r <- h_auto$counts
ymax_lat_r <- max(ni_r) * 1.2
hist(latitud, breaks = breaks_r, freq = TRUE,
main = "Histograma de la Latitud (Ajuste Automático R)",
xlab = "Latitud (grados)", ylab = "Frecuencia",
col = terrain.colors(length(breaks_r) - 1), border = "black", las = 1,
xlim = c(min(breaks_r), max(breaks_r)),
ylim = c(0, ymax_lat_r))
# Nueva Tabla
n <- length(latitud)
MC <- h_auto$mids
hi_perc <- (ni_r / n) * 100
Niasc <- cumsum(ni_r)
Nidsc <- rev(cumsum(rev(ni_r)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Completa
Tabla_Lat_R <- data.frame(
liminf = breaks_r[-length(breaks_r)],
limsup = breaks_r[-1],
MC = MC,
ni = ni_r,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_Lat_R)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 -80 -70 -75 2 0.06 2 3550 0.06 100.00
## 2 -70 -60 -65 1 0.03 3 3548 0.08 99.94
## 3 -60 -50 -55 2 0.06 5 3547 0.14 99.92
## 4 -50 -40 -45 0 0.00 5 3545 0.14 99.86
## 5 -40 -30 -35 5 0.14 10 3545 0.28 99.86
## 6 -30 -20 -25 7 0.20 17 3540 0.48 99.72
## 7 -20 -10 -15 9 0.25 26 3533 0.73 99.52
## 8 -10 0 -5 3 0.08 29 3524 0.82 99.27
## 9 0 10 5 11 0.31 40 3521 1.13 99.18
## 10 10 20 15 136 3.83 176 3510 4.96 98.87
## 11 20 30 25 1261 35.52 1437 3374 40.48 95.04
## 12 30 40 35 803 22.62 2240 2113 63.10 59.52
## 13 40 50 45 906 25.52 3146 1310 88.62 36.90
## 14 50 60 55 297 8.37 3443 404 96.99 11.38
## 15 60 70 65 97 2.73 3540 107 99.72 3.01
## 16 70 80 75 10 0.28 3550 10 100.00 0.281.3 Gráfica 3: Boxplot de latitud
boxplot(latitud, horizontal = TRUE, col = "beige",
main = "Diagrama de Caja de la Latitud", xlab = "Latitud (grados)")
1.4 Gráfica 4: Ojivas ascendentes y descendentes
ymax_ojiva_lat <- max(Tabla_Lat_R$Niasc) * 1.1
plot(Tabla_Lat_R$MC, Tabla_Lat_R$Niasc, type = "b", pch = 19, col = "black",
ylim = c(0, ymax_ojiva_lat), main = "Ojivas de Latitud (Ascendente y Descendente)",
xlab = "Latitud", ylab = "Frecuencia acumulada")
lines(Tabla_Lat_R$MC, Tabla_Lat_R$Nidsc, col = "blue", type = "b", pch = 19)
grid()
legend("topleft", legend = c("Ascendente", "Descendente"), col = c("black", "blue"),
lty = 1, pch = 19, cex = 0.8)
1.5 Indicadores estadísticos
# Función para moda
get_mode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
#Cálculos de todos los indicadores principales
media_lat <- mean(latitud)
mediana_lat <- median(latitud)
moda_lat <- get_mode(latitud)
varianza_lat <- var(latitud)
desv_estandar_lat <- sd(latitud)
coef_variacion_lat <- (desv_estandar_lat / media_lat) * 100
asimetria_lat <- skewness(latitud)
curtosis_lat <- kurtosis(latitud)
#Creación de la tabla consolidada
indicadores_lat <- data.frame(
Indicador = c("Moda (datos crudos)",
"Mediana",
"Media (x̄)",
"Desviación Estándar (σ)",
"Varianza (σ²)",
"Coef. Variación (%)",
"Asimetría",
"Curtosis"),
Valor = c(round(moda_lat, 2),
round(mediana_lat, 2),
round(media_lat, 2),
round(desv_estandar_lat, 2),
round(varianza_lat, 2),
round(coef_variacion_lat, 2),
round(asimetria_lat, 2),
round(curtosis_lat, 2))
)
# Impresión de la tabla
cat("Tabla de Indicadores: Latitud\n\n")## Tabla de Indicadores: Latitudprint(indicadores_lat)## Indicador Valor
## 1 Moda (datos crudos) 36.88
## 2 Mediana 34.78
## 3 Media (x̄) 36.32
## 4 Desviación Estándar (σ) 12.28
## 5 Varianza (σ²) 150.71
## 6 Coef. Variación (%) 33.80
## 7 Asimetría -1.07
## 8 Curtosis 9.14Conclusiones
El comportamiento de la variable latitud fluctúa entre -78° y 71.58°. Sus valores están en torno a una media de 36.32°, con una desviación estándar de 12.28, mostrando una distribución heterogénea. Sus valores se concentran fuertemente en el hemisferio norte, específicamente entre las latitudes de 25° y 50°. Los datos presentan un sesgo negativo debido a la presencia de valores atípicos que se extienden hacia el sur, entre -75° y -35°. Por lo tanto, el análisis del comportamiento de esta variable es altamente beneficioso, ya que nos permite identificar con claridad las zonas geográficas globales de mayor riesgo y ocurrencia de incidentes petroleros.
2 Variable Longitud
longitud_raw <- gsub(",", ".", df$longuitud)
longitud <- suppressWarnings(as.numeric(longitud_raw))
longitud <- na.omit(longitud)
# Limpieza de datos: Filtramos valores geográficamente imposibles
longitud <- longitud[longitud >= -180 & longitud <= 180]
# Regla de Sturges
n <- length(longitud)
k <- ceiling(1 + 3.322 * log10(n))
rango_min <- min(longitud)
rango_max <- max(longitud)
# Límites y Marca de Clase
breaks_long <- seq(from = rango_min, to = rango_max, length.out = k + 1)
liminf <- breaks_long[-(k + 1)]
limsup <- breaks_long[-1]
MC <- (liminf + limsup) / 2
# Agrupación y Frecuencias
rango_long <- cut(longitud, breaks = breaks_long, right = FALSE, include.lowest = TRUE)
ni <- as.numeric(table(rango_long))
hi_perc <- (ni / n) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Consolidada
Tabla_Long_Final <- data.frame(
liminf = round(liminf, 2),
limsup = round(limsup, 2),
MC = round(MC, 2),
ni = ni,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_Long_Final)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 -179.55 -151.99 -165.77 292 8.23 292 3548 8.23 100.00
## 2 -151.99 -124.44 -138.21 222 6.26 514 3256 14.49 91.77
## 3 -124.44 -96.88 -110.66 582 16.40 1096 3034 30.89 85.51
## 4 -96.88 -69.32 -83.10 2246 63.30 3342 2452 94.19 69.11
## 5 -69.32 -41.76 -55.54 124 3.49 3466 206 97.69 5.81
## 6 -41.76 -14.21 -27.99 3 0.08 3469 82 97.77 2.31
## 7 -14.21 13.35 -0.43 18 0.51 3487 79 98.28 2.23
## 8 13.35 40.91 27.13 9 0.25 3496 61 98.53 1.72
## 9 40.91 68.46 54.68 11 0.31 3507 52 98.84 1.47
## 10 68.46 96.02 82.24 1 0.03 3508 41 98.87 1.16
## 11 96.02 123.58 109.80 6 0.17 3514 40 99.04 1.13
## 12 123.58 151.13 137.36 24 0.68 3538 34 99.72 0.96
## 13 151.13 178.69 164.91 10 0.28 3548 10 100.00 0.282.1 Gráfica 1: Distribución de la longitud (Sturges)
ymax_long <- max(ni) * 1.2
hist(longitud, breaks = breaks_long, freq = TRUE,
main = "Histograma de la Longitud (Regla de Sturges)",
xlab = "Longitud (grados)", ylab = "Frecuencia",
col = terrain.colors(length(breaks_long) - 1), border = "black", las = 1,
xlim = c(min(breaks_long), max(breaks_long)),
ylim = c(0, ymax_long),
xaxt = "n")
axis(side = 1, at = breaks_long, labels = round(breaks_long, 2), las = 2, cex.axis = 0.8)
2.2 Gráfica 2: Distribución de la longitud (R)
# Extraer los cálculos automáticos de R
h_auto_long <- hist(longitud, plot = FALSE)
# R ya nos da los breaks limpios y el conteo exacto
breaks_r_long <- h_auto_long$breaks
ni_r_long <- h_auto_long$counts
ymax_long_r <- max(ni_r_long) * 1.2
hist(longitud, breaks = breaks_r_long, freq = TRUE,
main = "Histograma de la Longitud (Ajuste Automático R)",
xlab = "Longitud (grados)", ylab = "Frecuencia",
col = terrain.colors(length(breaks_r_long) - 1), border = "black", las = 1,
xlim = c(min(breaks_r_long), max(breaks_r_long)),
ylim = c(0, ymax_long_r))
# Nueva Tabla
n <- length(longitud)
MC <- h_auto_long$mids
hi_perc <- (ni_r_long / n) * 100
Niasc <- cumsum(ni_r_long)
Nidsc <- rev(cumsum(rev(ni_r_long)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
#Tabla completa
Tabla_Long_R <- data.frame(
liminf = breaks_r_long[-length(breaks_r_long)],
limsup = breaks_r_long[-1],
MC = MC,
ni = ni_r_long,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_Long_R)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 -180 -160 -170 162 4.57 162 3548 4.57 100.00
## 2 -160 -140 -150 207 5.83 369 3386 10.40 95.43
## 3 -140 -120 -130 609 17.16 978 3179 27.56 89.60
## 4 -120 -100 -110 95 2.68 1073 2570 30.24 72.44
## 5 -100 -80 -90 1519 42.81 2592 2475 73.06 69.76
## 6 -80 -60 -70 869 24.49 3461 956 97.55 26.94
## 7 -60 -40 -50 5 0.14 3466 87 97.69 2.45
## 8 -40 -20 -30 1 0.03 3467 82 97.72 2.31
## 9 -20 0 -10 14 0.39 3481 81 98.11 2.28
## 10 0 20 10 7 0.20 3488 67 98.31 1.89
## 11 20 40 30 8 0.23 3496 60 98.53 1.69
## 12 40 60 50 10 0.28 3506 52 98.82 1.47
## 13 60 80 70 1 0.03 3507 42 98.84 1.18
## 14 80 100 90 1 0.03 3508 41 98.87 1.16
## 15 100 120 110 2 0.06 3510 40 98.93 1.13
## 16 120 140 130 9 0.25 3519 38 99.18 1.07
## 17 140 160 150 20 0.56 3539 29 99.75 0.82
## 18 160 180 170 9 0.25 3548 9 100.00 0.252.3 Gráfica 3: Boxplot de longitud
boxplot(longitud, horizontal = TRUE, col = "beige",
main = "Diagrama de Caja de la Longitud", xlab = "Longitud (grados)")
2.4 Gráfica 4: Ojiva ascendentes y descendentes
ymax_ojiva_long <- max(Tabla_Long_R$Niasc) * 1.1
plot(Tabla_Long_R$MC, Tabla_Long_R$Niasc, type = "b", pch = 19, col = "black",
ylim = c(0, ymax_ojiva_long), main = "Ojivas de Longitud (Ascendente y Descendente)",
xlab = "Longitud", ylab = "Frecuencia acumulada")
lines(Tabla_Long_R$MC, Tabla_Long_R$Nidsc, col = "blue", type = "b", pch = 19)
grid()
legend("topleft", legend = c("Ascendente", "Descendente"), col = c("black", "blue"),
lty = 1, pch = 19, cex = 0.8)
2.5 Indicadores estadísticos
# Cálculos de todos los indicadores principales
media_long <- mean(longitud)
mediana_long <- median(longitud)
moda_long <- get_mode(longitud) # Usa la función get_mode creada en la sección de Latitud
varianza_long <- var(longitud)
desv_estandar_long <- sd(longitud)
coef_variacion_long <- (desv_estandar_long / media_long) * 100
asimetria_long <- skewness(longitud)
curtosis_long <- kurtosis(longitud)
# Creación de la tabla consolidada
indicadores_long <- data.frame(
Indicador = c("Moda (datos crudos)",
"Mediana",
"Media (x̄)",
"Desviación Estándar (σ)",
"Varianza (σ²)",
"Coef. Variación (%)",
"Asimetría",
"Curtosis"),
Valor = c(round(moda_long, 2),
round(mediana_long, 2),
round(media_long, 2),
round(desv_estandar_long, 2),
round(varianza_long, 2),
round(coef_variacion_long, 2),
round(asimetria_long, 2),
round(curtosis_long, 2))
)
# Impresión de la tabla
cat("Tabla de Indicadores: Longitud\n\n")## Tabla de Indicadores: Longitudprint(indicadores_long)## Indicador Valor
## 1 Moda (datos crudos) -76.33
## 2 Mediana -90.06
## 3 Media (x̄) -95.57
## 4 Desviación Estándar (σ) 39.77
## 5 Varianza (σ²) 1581.30
## 6 Coef. Variación (%) -41.61
## 7 Asimetría 2.42
## 8 Curtosis 14.83Conclusiones
El comportamiento de la variable longitud fluctúa entre -179.55° y 178.69°. Sus valores están en torno a una media de -95.57°, con una desviación estándar de 39.77, mostrando una distribución altamente dispersa. Sus valores se concentran fuertemente en el hemisferio occidental en torno a las longitudes que enmarcan el continente americano. Los datos presentan un sesgo positivo debido a la presencia de valores atípicos que se extienden hacia el hemisferio oriental, entre 5° y 165°. Por lo tanto, el análisis del comportamiento de esta variable es altamente beneficioso, ya que, en conjunto con la latitud, nos permite identificar con precisión las coordenadas globales de mayor riesgo y ocurrencia de incidentes petroleros.
3 Variable Máxima Liberación de Petróleo
max_lib_raw <- as.numeric(as.character(df$maximo_liberacion_galones))
max_lib <- na.omit(max_lib_raw)
# Limpieza de datos
max_lib <- max_lib[max_lib >= 0]
# Regla de Sturges y Amplitud
n <- length(max_lib)
k <- ceiling(1 + 3.322 * log10(n))
R <- max(max_lib) - min(max_lib)
A <- R / k
# Límites y Marca de Clase (MC)
liminf <- seq(from = min(max_lib), by = A, length.out = k)
limsup <- liminf + A
# Ajuste del último límite para asegurar la inclusión del valor máximo
breaks_max_lib <- c(liminf, max(limsup) + 0.1)
MC <- (liminf + limsup) / 2
# Agrupación y Frecuencias
rango_max_lib <- cut(max_lib, breaks = breaks_max_lib, right = FALSE, include.lowest = TRUE)
ni <- as.numeric(table(rango_max_lib))
hi_perc <- (ni / n) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Consolidada
Tabla_MaxLib_Final <- data.frame(
liminf = round(liminf, 2),
limsup = round(limsup, 2),
MC = round(MC, 2),
ni = ni,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_MaxLib_Final)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 0 25846155 12923077 2158 99.54 2158 2168 99.54 100.00
## 2 25846155 51692309 38769232 7 0.32 2165 10 99.86 0.46
## 3 51692309 77538464 64615386 1 0.05 2166 3 99.91 0.14
## 4 77538464 103384618 90461541 0 0.00 2166 2 99.91 0.09
## 5 103384618 129230773 116307695 0 0.00 2166 2 99.91 0.09
## 6 129230773 155076927 142153850 0 0.00 2166 2 99.91 0.09
## 7 155076927 180923082 168000005 0 0.00 2166 2 99.91 0.09
## 8 180923082 206769236 193846159 1 0.05 2167 2 99.95 0.09
## 9 206769236 232615391 219692314 0 0.00 2167 1 99.95 0.05
## 10 232615391 258461545 245538468 0 0.00 2167 1 99.95 0.05
## 11 258461545 284307700 271384623 0 0.00 2167 1 99.95 0.05
## 12 284307700 310153854 297230777 0 0.00 2167 1 99.95 0.05
## 13 310153854 336000009 323076932 1 0.05 2168 1 100.00 0.053.1 Gráfica 1: Distribución de Máxima liberación (Sturges)
ymax_maxlib <- max(Tabla_MaxLib_Final$Niasc) * 1.1
old_par <- par(mar = c(9, 5, 4, 2))
hist(max_lib, breaks = breaks_max_lib, freq = TRUE,
main = "Histograma de la Máxima Lib. de Petróleo (Sturges)",
xlab = "", ylab = "Frecuencia",
col = terrain.colors(length(breaks_max_lib) - 1), border = "black", las = 1,
ylim = c(0, ymax_maxlib),
xaxt = "n")
axis(side = 1, at = breaks_max_lib, labels = round(breaks_max_lib, 2), las = 2, cex.axis = 0.8)
mtext("Máxima liberación (galones)", side = 1, line = 7)
par(old_par)3.2 Gráfica 2: Distribución de Máxima liberación (R)
# Extraer los cálculos automáticos de R
h_auto_maxlib <- hist(max_lib, plot = FALSE)
# R ya nos da los breaks limpios y el conteo exacto
breaks_r_maxlib <- h_auto_maxlib$breaks
ni_r_maxlib <- h_auto_maxlib$counts
ymax_maxlib_r <- max(ni_r_maxlib) * 1.2
hist(max_lib, breaks = breaks_r_maxlib, freq = TRUE,
main = "Histograma de la Máxima Liberación (Ajuste Automático R)",
xlab = "Máxima liberación (galones)", ylab = "Frecuencia",
col = terrain.colors(length(breaks_r_maxlib) - 1), border = "black", las = 1,
xlim = c(min(breaks_r_maxlib), max(breaks_r_maxlib)),
ylim = c(0, ymax_maxlib_r))
# Nueva Tabla
n <- length(max_lib)
MC <- h_auto_maxlib$mids
hi_perc <- (ni_r_maxlib / n) * 100
Niasc <- cumsum(ni_r_maxlib)
Nidsc <- rev(cumsum(rev(ni_r_maxlib)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Completa
Tabla_MaxLib_R <- data.frame(
liminf = breaks_r_maxlib[-length(breaks_r_maxlib)],
limsup = breaks_r_maxlib[-1],
MC = MC,
ni = ni_r_maxlib,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_MaxLib_R)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 0.0e+00 2.0e+07 1.0e+07 2151 99.22 2151 2168 99.22 100.00
## 2 2.0e+07 4.0e+07 3.0e+07 13 0.60 2164 17 99.82 0.78
## 3 4.0e+07 6.0e+07 5.0e+07 1 0.05 2165 4 99.86 0.18
## 4 6.0e+07 8.0e+07 7.0e+07 1 0.05 2166 3 99.91 0.14
## 5 8.0e+07 1.0e+08 9.0e+07 0 0.00 2166 2 99.91 0.09
## 6 1.0e+08 1.2e+08 1.1e+08 0 0.00 2166 2 99.91 0.09
## 7 1.2e+08 1.4e+08 1.3e+08 0 0.00 2166 2 99.91 0.09
## 8 1.4e+08 1.6e+08 1.5e+08 0 0.00 2166 2 99.91 0.09
## 9 1.6e+08 1.8e+08 1.7e+08 0 0.00 2166 2 99.91 0.09
## 10 1.8e+08 2.0e+08 1.9e+08 0 0.00 2166 2 99.91 0.09
## 11 2.0e+08 2.2e+08 2.1e+08 1 0.05 2167 2 99.95 0.09
## 12 2.2e+08 2.4e+08 2.3e+08 0 0.00 2167 1 99.95 0.05
## 13 2.4e+08 2.6e+08 2.5e+08 0 0.00 2167 1 99.95 0.05
## 14 2.6e+08 2.8e+08 2.7e+08 0 0.00 2167 1 99.95 0.05
## 15 2.8e+08 3.0e+08 2.9e+08 0 0.00 2167 1 99.95 0.05
## 16 3.0e+08 3.2e+08 3.1e+08 0 0.00 2167 1 99.95 0.05
## 17 3.2e+08 3.4e+08 3.3e+08 1 0.05 2168 1 100.00 0.053.3 Gráfica 3: Boxplot de Máxima liberación
boxplot(max_lib, horizontal = TRUE, col = "beige",
main = "Diagrama de Caja de la Máxima Liberación",
xlab = "Máxima liberación (galones)")
3.4 Gráfica 4: Ojiva ascendentes y descendentes
ymax_ojiva_maxlib <- max(Tabla_MaxLib_R$Niasc) * 1.1
plot(Tabla_MaxLib_R$MC, Tabla_MaxLib_R$Niasc, type = "b", pch = 19, col = "black",
ylim = c(0, ymax_ojiva_maxlib), main = "Ojivas de Máxima Liberación (Ascendente y Descendente)",
xlab = "Máxima liberación (galones)", ylab = "Frecuencia acumulada")
lines(Tabla_MaxLib_R$MC, Tabla_MaxLib_R$Nidsc, col = "blue", type = "b", pch = 19)
grid()
legend("right", legend = c("Ascendente", "Descendente"), col = c("black", "blue"),
lty = 1, pch = 19, cex = 0.8)
3.5 Indicadores estadísticos
# Cálculos de los indicadores principales
media_maxlib <- mean(max_lib)
mediana_maxlib <- median(max_lib)
moda_maxlib <- get_mode(max_lib)
varianza_maxlib <- var(max_lib)
desv_estandar_maxlib <- sd(max_lib)
coef_variacion_maxlib <- (desv_estandar_maxlib / media_maxlib) * 100
asimetria_maxlib <- skewness(max_lib)
curtosis_maxlib <- kurtosis(max_lib)
# Creación de la tabla consolidada
indicadores_maxlib <- data.frame(
Indicador = c("Moda (datos crudos)",
"Mediana",
"Media (x̄)",
"Desviación Estándar (σ)",
"Varianza (σ²)",
"Coef. Variación (%)",
"Asimetría",
"Curtosis"),
Valor = c(round(moda_maxlib, 2),
round(mediana_maxlib, 2),
round(media_maxlib, 2),
round(desv_estandar_maxlib, 2),
round(varianza_maxlib, 2),
round(coef_variacion_maxlib, 2),
round(asimetria_maxlib, 2),
round(curtosis_maxlib, 2))
)
# Impresión de la tabla
cat("Tabla de Indicadores: Máxima Liberación\n\n")## Tabla de Indicadores: Máxima Liberaciónprint(indicadores_maxlib)## Indicador Valor
## 1 Moda (datos crudos) 1.000000e+03
## 2 Mediana 3.150000e+03
## 3 Media (x̄) 7.896504e+05
## 4 Desviación Estándar (σ) 9.010100e+06
## 5 Varianza (σ²) 8.118191e+13
## 6 Coef. Variación (%) 1.141020e+03
## 7 Asimetría 2.958000e+01
## 8 Curtosis 1.004690e+03Conclusiones
El comportamiento de la variable máxima liberación de petróleo fluctúa entre 0 y 336000009 galones. Sus valores están en torno a una mediana de 3150 galones, con una desviación estándar de 9010100 galones, mostrando una distribución extremadamente dispersa y heterogénea. Sus valores se concentran fuertemente en los rangos más bajos. Los datos presentan un pronunciado sesgo positivo debido a la presencia de valores atípicos extremos que se extienden hacia la derecha, representando derrames masivos de entre 50000000 y 300000000 galones. Por lo tanto, el análisis del comportamiento de esta variable es altamente beneficioso, ya que nos permite identificar que, aunque la inmensa mayoría de los incidentes son de pequeña magnitud, el riesgo global está fuertemente impactado por unos pocos eventos de proporciones catastróficas.
4 Variable Respuesta actual
Dado que la variable Respuesta actual tiene una cola larga hacia la derecha, casi todos los datos se concentran en el intervalo inicial de 0 a 40.000, lo que oculta los detalles. Al aplicar nuevamente la regla de Sturges solo a este bloque denso, logramos ver cómo se distribuyen realmente las respuestas más comunes, evitamos que los valores atípicos distorsionen nuestros gráficos y sacamos a la luz la variación que estaba oculta en el grueso de los datos.
# Extracción y limpieza de la variable
Resp_actual_raw <- as.numeric(as.character(datos$Respuesta_actual_galones))## Warning: NAs introducidos por coerciónResp_actual <- na.omit(Resp_actual_raw)
# Definir rango fijo (0 - 40,000)
Resp_actual_rango1 <- Resp_actual[Resp_actual >= 0 & Resp_actual <= 40000]
# Regla de Sturges y Amplitud sobre los datos filtrados
n <- length(Resp_actual_rango1)
k <- round(1 + 3.322 * log10(n))
R <- max(Resp_actual_rango1) - min(Resp_actual_rango1)
A <- R / k
# Límites y Marca de Clase (MC)
liminf <- seq(from = min(Resp_actual_rango1), by = A, length.out = k)
limsup <- liminf + A
# Ajuste del último límite para asegurar la inclusión del valor máximo
breaks_resp_actual <- c(liminf, max(limsup) + 0.1)
MC <- (liminf + limsup) / 2
# Agrupación y Frecuencias
rango_resp_actual <- cut(Resp_actual_rango1, breaks = breaks_resp_actual, right = FALSE, include.lowest = TRUE)
ni <- as.numeric(table(rango_resp_actual))
hi_perc <- (ni / n) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Consolidada
Tabla_RespActual_Final <- data.frame(
liminf = round(liminf, 2),
limsup = round(limsup, 2),
MC = round(MC, 2),
ni = ni,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_RespActual_Final)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 0.00 3333.33 1666.67 1134 74.90 1134 1514 74.90 100.00
## 2 3333.33 6666.67 5000.00 138 9.11 1272 380 84.02 25.10
## 3 6666.67 10000.00 8333.33 59 3.90 1331 242 87.91 15.98
## 4 10000.00 13333.33 11666.67 59 3.90 1390 183 91.81 12.09
## 5 13333.33 16666.67 15000.00 22 1.45 1412 124 93.26 8.19
## 6 16666.67 20000.00 18333.33 14 0.92 1426 102 94.19 6.74
## 7 20000.00 23333.33 21666.67 20 1.32 1446 88 95.51 5.81
## 8 23333.33 26666.67 25000.00 21 1.39 1467 68 96.90 4.49
## 9 26666.67 30000.00 28333.33 7 0.46 1474 47 97.36 3.10
## 10 30000.00 33333.33 31666.67 18 1.19 1492 40 98.55 2.64
## 11 33333.33 36666.67 35000.00 12 0.79 1504 22 99.34 1.45
## 12 36666.67 40000.00 38333.33 10 0.66 1514 10 100.00 0.664.1 Gráfica 1: Distribución de respuesta actual (Struges)
ymax_resp <- max(Tabla_RespActual_Final$Niasc) * 1.1
old_par <- par(mar = c(8, 5, 4, 2))
hist(Resp_actual_rango1, breaks = breaks_resp_actual, freq = TRUE,
main = "Histograma de la Respuesta Actual (Rango 1 - Sturges)",
xlab = "", ylab = "Frecuencia",
col = "purple2", border = "black", las = 1,
ylim = c(0, ymax_resp),
xaxt = "n")## Warning in plot.histogram(r, freq = freq1, col = col, border = border, angle =
## angle, : the AREAS in the plot are wrong -- rather use 'freq = FALSE'axis(side = 1, at = breaks_resp_actual, labels = round(breaks_resp_actual, 2), las = 2, cex.axis = 0.8)
mtext("Cantidad recolectada (galones)", side = 1, line = 6)
par(old_par)4.2 Gráfica 2: Distribución de respuesta actual (R)
# Extraer los cálculos automáticos de R
h_auto_resp <- hist(Resp_actual_rango1, plot = FALSE)
# R ya nos da los breaks limpios y el conteo exacto
breaks_r_resp <- h_auto_resp$breaks
ni_r_resp <- h_auto_resp$counts
ymax_resp_r <- max(ni_r_resp) * 1.2
hist(Resp_actual_rango1, breaks = breaks_r_resp, freq = TRUE,
main = "Histograma de la Respuesta Actual (Ajuste Automático R)",
xlab = "Cantidad recolectada (galones)", ylab = "Frecuencia",
col = "purple2", border = "black", las = 1,
xlim = c(min(breaks_r_resp), max(breaks_r_resp)),
ylim = c(0, ymax_resp_r))
# Nueva Tabla
n <- length(Resp_actual_rango1)
MC <- h_auto_resp$mids
hi_perc <- (ni_r_resp / n) * 100
Niasc <- cumsum(ni_r_resp)
Nidsc <- rev(cumsum(rev(ni_r_resp)))
Hiasc_perc <- cumsum(hi_perc)
Hidsc_perc <- rev(cumsum(rev(hi_perc)))
# Tabla Completa
Tabla_RespActual_R <- data.frame(
liminf = breaks_r_resp[-length(breaks_r_resp)],
limsup = breaks_r_resp[-1],
MC = MC,
ni = ni_r_resp,
hi_perc = round(hi_perc, 2),
Niasc = Niasc,
Nidsc = Nidsc,
Hiasc_perc = round(Hiasc_perc, 2),
Hidsc_perc = round(Hidsc_perc, 2)
)
print(Tabla_RespActual_R)## liminf limsup MC ni hi_perc Niasc Nidsc Hiasc_perc Hidsc_perc
## 1 0 5000 2500 1251 82.63 1251 1514 82.63 100.00
## 2 5000 10000 7500 102 6.74 1353 263 89.37 17.37
## 3 10000 15000 12500 51 3.37 1404 161 92.73 10.63
## 4 15000 20000 17500 31 2.05 1435 110 94.78 7.27
## 5 20000 25000 22500 23 1.52 1458 79 96.30 5.22
## 6 25000 30000 27500 29 1.92 1487 56 98.22 3.70
## 7 30000 35000 32500 16 1.06 1503 27 99.27 1.78
## 8 35000 40000 37500 11 0.73 1514 11 100.00 0.734.3 Gráfica 2: Boxplot de respuesta actual
boxplot(Resp_actual_rango1, horizontal = TRUE, col = "purple2",
main = "Diagrama de Caja de la Respuesta Actual",
xlab = "Cantidad recolectada (galones)")
4.4 Gráfica 3: Ojiva ascendentes y descendentes
ymax_ojiva_resp <- max(Tabla_RespActual_R$Niasc) * 1.1
plot(Tabla_RespActual_R$MC, Tabla_RespActual_R$Niasc, type = "b", pch = 19, col = "black",
ylim = c(0, ymax_ojiva_resp), main = "Ojivas de Respuesta Actual (Ascendente y Descendente)",
xlab = "Cantidad recolectada (galones)", ylab = "Frecuencia acumulada")
lines(Tabla_RespActual_R$MC, Tabla_RespActual_R$Nidsc, col = "blue", type = "b", pch = 19)
grid()
legend("right", legend = c("Ascendente", "Descendente"), col = c("black", "blue"),
lty = 1, pch = 19, cex = 0.8)
4.5 Indicadores estadísticos
get_mode <- function(v) {
v <- v[v != 0 & !is.na(v)]
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
moda_resp <- get_mode(Resp_actual_rango1)
# Cálculos de los indicadores principales
media_resp <- mean(Resp_actual_rango1)
mediana_resp <- median(Resp_actual_rango1)
varianza_resp <- var(Resp_actual_rango1)
desv_estandar_resp <- sd(Resp_actual_rango1)
coef_variacion_resp <- (desv_estandar_resp / media_resp) * 100
asimetria_resp <- skewness(Resp_actual_rango1)
curtosis_resp <- kurtosis(Resp_actual_rango1)
# Creación de la tabla consolidada
indicadores_resp <- data.frame(
Indicador = c("Moda (datos crudos)",
"Mediana",
"Media (x̄)",
"Desviación Estándar (σ)",
"Varianza (σ²)",
"Coef. Variación (%)",
"Asimetría",
"Curtosis"),
Valor = c(round(moda_resp, 2),
round(mediana_resp, 2),
round(media_resp, 2),
round(desv_estandar_resp, 2),
round(varianza_resp, 2),
round(coef_variacion_resp, 2),
round(asimetria_resp, 2),
round(curtosis_resp, 2))
)
# Impresión de la tabla
cat("Tabla de Indicadores: Respuesta Actual (Rango 1: 0–40,000 galones)\n\n")## Tabla de Indicadores: Respuesta Actual (Rango 1: 0–40,000 galones)print(indicadores_resp)## Indicador Valor
## 1 Moda (datos crudos) 1000.00
## 2 Mediana 500.00
## 3 Media (x̄) 3599.25
## 4 Desviación Estándar (σ) 7195.78
## 5 Varianza (σ²) 51779178.24
## 6 Coef. Variación (%) 199.92
## 7 Asimetría 2.87
## 8 Curtosis 8.31Conclusiones
El comportamiento de la variable respuesta actual en el rango 1 fluctúa entre 0 y 40000 galones. Sus valores están en torno a una mediana de 500 galones, con una desviación estándar de 7195.78 galones, mostrando una distribución altamente dispersa y heterogénea.Los datos presentan un claro sesgo positivo debido a la presencia de valores atípicos que se extienden hacia la derecha, representando respuestas o recuperaciones de mayor magnitud situadas entre 9000 y 40000 galones. Por lo cual, el análisis del comportamiento de esta variable es altamente beneficioso, ya que nos permite identificar que, aunque la operatividad habitual implica la recolección de volúmenes pequeños, los protocolos de respuesta deben estar dimensionados para gestionar eventos atípicos de contención masiva.