1.Carga de Datos y Librerías

library(readxl)
library(dplyr)
library(tidyr)
library(gt)
library(ggplot2)
library(scales)
library(forcats)
library(MASS)   # Fundamental para el ajuste de distribuciones (fitdistr)
datos <- read_excel("dataset_mundial_petro.xlsx")
cat("Número de registros:", nrow(datos), "\n")
## Número de registros: 49212
cat("Número de variables:", ncol(datos), "\n")
## Número de variables: 32

2.Extracción de Variable Aleatoria

La variable Longitude indica la coordenada geográfica de longitud de cada yacimiento. Es una variable cuantitativa continua, cuyos valores van desde coordenadas negativas (Hemisferio Occidental) hasta positivas (Hemisferio Oriental). Por su naturaleza bimodal ligada al hemisferio, se le aplica un tratamiento analítico segmentado.

n_total <- nrow(datos)
Variable <- na.omit(as.numeric(datos$Longitude))
n <- length(Variable)

cat("Número de registros (n):", n, "\n")
## Número de registros (n): 7537
cat("Valor mínimo:", round(min(Variable), 3), "\n")
## Valor mínimo: -152.129
cat("Valor máximo:", round(max(Variable), 3), "\n")
## Valor máximo: 174.361
n_occidental <- sum(Variable <= 0)
n_oriental   <- sum(Variable > 0)

cat("Registros en Hemisferio Occidental (<= 0°):", n_occidental,
    "(", round(n_occidental / n * 100, 2), "%)\n")
## Registros en Hemisferio Occidental (<= 0°): 5290 ( 70.19 %)
cat("Registros en Hemisferio Oriental (> 0°):", n_oriental,
    "(", round(n_oriental / n * 100, 2), "%)\n")
## Registros en Hemisferio Oriental (> 0°): 2247 ( 29.81 %)

3.Tabla de Distribución de Frecuencias

Se aplica la Regla de Sturges para determinar el número óptimo de intervalos de clase, ajustando los límites a múltiplos de 10 para facilitar la lectura.

BASE <- 10
min_int <- floor(min(Variable) / BASE) * BASE
max_int <- ceiling(max(Variable) / BASE) * BASE
k_sug   <- floor(1 + 3.322 * log10(n))

Rango_int    <- max_int - min_int
Amplitud_int <- ceiling((Rango_int / k_sug) / 10) * 10
if (Amplitud_int == 0) Amplitud_int <- 10

cortes_int <- seq(from = min_int, by = Amplitud_int, length.out = k_sug + 1)
if (max(cortes_int) < max(Variable)) cortes_int <- c(cortes_int, max(cortes_int) + Amplitud_int)

while (length(cortes_int) > 2 && cortes_int[length(cortes_int) - 1] >= max(Variable)) {
  cortes_int <- cortes_int[-length(cortes_int)]
}

k <- length(cortes_int) - 1
inter_int <- cut(Variable, breaks = cortes_int, include.lowest = TRUE, right = FALSE)
ni_int    <- as.vector(table(inter_int))

conteo <- data.frame(
  Li = cortes_int[1:k],
  Ls = cortes_int[2:(k + 1)],
  MC = (cortes_int[1:k] + cortes_int[2:(k + 1)]) / 2,
  ni = ni_int
) %>%
  mutate(
    hi     = round(ni / n, 4),
    hi_pct = round(ni / n * 100, 2)
  )

cat("Total de intervalos de clase:", k, "\n")
## Total de intervalos de clase: 12
cat("Amplitud de clase:", Amplitud_int, "\n")
## Amplitud de clase: 30
cat("Intervalo más frecuente: [", conteo$Li[which.max(conteo$ni)], ",",
    conteo$Ls[which.max(conteo$ni)], ") con", max(conteo$ni), "registros\n")
## Intervalo más frecuente: [ -130 , -100 ) con 2689 registros
fila_total <- tibble(
  Li = NA, Ls = NA, MC = NA,
  ni     = sum(conteo$ni),
  hi     = sum(conteo$hi),
  hi_pct = sum(conteo$hi_pct)
)

tdf_final <- bind_rows(conteo, fila_total)

tdf_final %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N° 1**"),
    subtitle = md("Distribución de yacimientos según Longitud Geográfica")
  ) %>%
  cols_label(
    Li     = "Lím. Inf (°)",
    Ls     = "Lím. Sup (°)",
    MC     = "Marca de Clase",
    ni     = "Frecuencia (ni)",
    hi     = "Proporción (hi)",
    hi_pct = "Porcentaje (hi%)"
  ) %>%
  fmt_number(columns = MC,     decimals = 1) %>%
  fmt_number(columns = hi,     decimals = 4) %>%
  fmt_number(columns = hi_pct, decimals = 2) %>%
  sub_missing(columns = c(Li, Ls, MC), missing_text = "TOTAL") %>%
  tab_source_note(source_note = "Autor: Grupo 5") %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.border.top.color            = "black",
    table.border.bottom.color         = "black",
    table.border.top.style            = "solid",
    table.border.bottom.style         = "solid",
    column_labels.font.weight         = "bold",
    column_labels.border.top.color    = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color       = "black",
    heading.border.bottom.width       = px(2),
    table_body.hlines.color           = "grey",
    table_body.border.bottom.color    = "black"
  )
Tabla N° 1
Distribución de yacimientos según Longitud Geográfica
Lím. Inf (°) Lím. Sup (°) Marca de Clase Frecuencia (ni) Proporción (hi) Porcentaje (hi%)
-160 -130 −145.0 37 0.0049 0.49
-130 -100 −115.0 2689 0.3568 35.68
-100 -70 −85.0 2018 0.2677 26.77
-70 -40 −55.0 403 0.0535 5.35
-40 -10 −25.0 51 0.0068 0.68
-10 20 5.0 1130 0.1499 14.99
20 50 35.0 431 0.0572 5.72
50 80 65.0 384 0.0509 5.09
80 110 95.0 183 0.0243 2.43
110 140 125.0 178 0.0236 2.36
140 170 155.0 26 0.0034 0.34
170 200 185.0 7 0.0009 0.09
TOTAL TOTAL TOTAL 7537 0.9999 99.99
Autor: Grupo 5

4.Análisis Gráfico

colores <- colorRampPalette(c("#2E86C1", "#AED6F1"))(k)

conteo <- conteo %>% mutate(densidad = round(hi / Amplitud_int, 5))

ggplot(conteo, aes(x = MC, y = densidad, fill = factor(MC))) +
  geom_col(width = Amplitud_int, color = "white") +
  geom_text(aes(label = round(densidad, 4)), vjust = -0.4, size = 2.8, fontface = "bold") +
  scale_fill_manual(values = colores) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title   = "Gráfica N°1: Histograma — Densidad de Probabilidad de Longitud",
    x       = "Longitud (°)",
    y       = "Densidad de Probabilidad",
    caption = paste0("n = ", format(n, big.mark = ","), " | Fuente: GOGET")
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold"),
        axis.title = element_text(face = "bold"))

5.Histograma General

Antes de proponer cualquier modelo, se visualiza el histograma general de la Longitud, usando los mismos intervalos calculados en la Sección 3.

grises <- gray(seq(0.35, 0.85, length.out = k))
h_gen  <- hist(Variable, breaks = cortes_int, plot = FALSE)
h_gen$density <- conteo$hi

par(mar = c(5, 6, 4, 2))
plot(h_gen, col = grises, border = "black", freq = FALSE,
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = cortes_int, labels = round(cortes_int, 0), las = 2, cex.axis = 0.8)
mtext("Densidad de Probabilidad", side = 2, line = 4.2, cex = 0.9)
mtext("Longitud (grados)", side = 1, line = 3.5, cex = 0.9)
mtext("Histograma General — Longitud de yacimientos de petróleo y gas (GOGET)",
      side = 3, line = 1.5, cex = 0.95, font = 2)

6.Conjetura de Modelo Exponencial por Zonas

El histograma general muestra un pico de concentración de yacimientos justo en torno a los -100°, con una caída marcada tanto hacia la izquierda (Longitudes más negativas) como hacia la derecha (Longitudes crecientes) — un comportamiento que en cada tramo se describe bien con una distribución Exponencial, aunque orientada en direcciones opuestas a cada lado del pico. Se reutilizan los mismos intervalos de la Regla de Sturges calculados en la Sección 3 (Tabla N°1), agrupados en 2 zonas:

  • Intervalo 1 [-160° , -130°): concentra un número de registros demasiado bajo para aportar información relevante al ajuste, por lo que se omite del análisis inferencial formal (ver nota en la Sección 6.1).
  • Zona 1 [-130° , -100°) (Intervalo 2 de Sturges): la densidad crece a medida que nos acercamos al límite superior (el pico general está justo ahí), así que se ajusta un modelo Exponencial con el pico en el extremo derecho (decae de derecha a izquierda).
  • Zona 2 [-100° , 200°] (Intervalos 3 a 12): la densidad decae de forma natural desde su extremo izquierdo, así que se ajusta un modelo Exponencial estándar (decae de izquierda a derecha).
intervalo1_inf <- cortes_int[1]
intervalo1_sup <- cortes_int[2]
zona1_inf      <- cortes_int[2]
zona1_sup      <- cortes_int[3]
zona2_inf      <- cortes_int[3]
zona2_sup      <- cortes_int[k + 1]

n_omitido <- sum(Variable >= intervalo1_inf & Variable < intervalo1_sup)

cat("Intervalo 1 omitido [", intervalo1_inf, ",", intervalo1_sup, "):", n_omitido,
    "registros (", round(n_omitido / n * 100, 2), "% del total)\n")
## Intervalo 1 omitido [ -160 , -130 ): 37 registros ( 0.49 % del total)
# Ajusta un modelo Exponencial a los datos de una zona.
# pico = "izquierda": la densidad decae de izquierda a derecha (estandar).
# pico = "derecha":   la densidad decae de derecha a izquierda (el pico
#                     esta pegado al limite superior de la zona).
# En ambos casos se desplaza la variable a dominio positivo antes de
# ajustar; el desplazamiento se deshace solo para graficar.
ajustar_zona <- function(datos, pico = "izquierda") {
  if (pico == "izquierda") {
    offset <- floor(min(datos)) - 1
    d_pos  <- datos - offset
  } else {
    offset <- ceiling(max(datos)) + 1
    d_pos  <- offset - datos
  }
  f <- fitdistr(d_pos, "exponential")
  list(fit = f, offset = offset, pico = pico)
}

# Frecuencia observada (hi_obs) y esperada (hi_teo) de una zona, en una
# malla de k_zona intervalos propios.
frecuencias_zona <- function(datos, res, lim_min, lim_max, k_zona = 4) {
  brks   <- seq(lim_min, lim_max, length.out = k_zona + 1)
  hi_obs <- hist(datos, breaks = brks, plot = FALSE)$counts / length(datos)
  mc     <- (head(brks, -1) + tail(brks, -1)) / 2
  mc_pos <- if (res$pico == "izquierda") mc - res$offset else res$offset - mc
  mc_pos[mc_pos <= 0] <- 1e-9
  hi_teo <- dexp(mc_pos, rate = res$fit$estimate["rate"]) * diff(brks)
  hi_teo <- hi_teo / sum(hi_teo)
  list(hi_obs = hi_obs, hi_teo = hi_teo)
}

# Grafica el histograma de una zona con su curva Exponencial superpuesta.
plot_zona <- function(datos, res, titulo, col_fill, lim_min, lim_max, k_zona = 4) {
  brks <- seq(lim_min, lim_max, length.out = k_zona + 1)
  h    <- hist(datos, breaks = brks, plot = FALSE)
  h$density <- (h$counts / length(datos)) / diff(brks)

  xs     <- seq(lim_min, lim_max, length.out = 500)
  xs_pos <- if (res$pico == "izquierda") xs - res$offset else res$offset - xs
  xs_pos[xs_pos <= 0] <- 1e-9
  ys <- dexp(xs_pos, rate = res$fit$estimate["rate"])

  par(mar = c(5, 6, 4, 2))
  plot(h, col = col_fill, border = "black", freq = FALSE,
       main = "", xlab = "", ylab = "", las = 1)
  lines(xs, ys, col = "#C0392B", lwd = 2.5)
  mtext("Densidad de Probabilidad", side = 2, line = 4.2, cex = 0.9)
  mtext("Longitud (grados)", side = 1, line = 3.5, cex = 0.9)
  mtext(titulo, side = 3, line = 1, cex = 0.9, font = 2)
  legend("topright",
         legend = c("Histograma", "Curva Exponencial"),
         fill = c(col_fill, NA), border = c("black", NA),
         lty = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.8)
}

6.1 Histograma con Cortes por Zona

lim_min_z <- intervalo1_inf
lim_max_z <- zona2_sup

brks_zon <- cortes_int
x_zonif  <- Variable[Variable >= lim_min_z & Variable <= lim_max_z]
n_zonif  <- length(x_zonif)

h_zon <- hist(x_zonif, breaks = brks_zon, plot = FALSE)
h_zon$density <- (h_zon$counts / n_zonif) / diff(brks_zon)

mc_zon <- (head(brks_zon, -1) + tail(brks_zon, -1)) / 2
colores_zona <- ifelse(mc_zon < zona1_inf, "gray85",
                 ifelse(mc_zon < zona1_sup, "gray45", "#AED6F1"))

par(mar = c(5, 6, 6, 2))
plot(h_zon, col = colores_zona, border = "black", freq = FALSE,
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = round(brks_zon, 0), labels = round(brks_zon, 0), las = 2, cex.axis = 0.75)
abline(v = c(zona1_inf, zona1_sup), col = "black", lty = 2, lwd = 2)

legend("topright",
       legend = c("Intervalo 1 (omitido)", "Zona 1 (Exponencial)", "Zona 2 (Exponencial)"),
       fill = c("gray85", "gray45", "#AED6F1"), border = "black",
       bty = "n", cex = 0.8)

mtext("Densidad de Probabilidad", side = 2, line = 4.2, cex = 0.9)
mtext("Longitud (grados)", side = 1, line = 3.8, cex = 0.9)
mtext("Cortes por Zona — Longitud de yacimientos de petróleo y gas (GOGET)",
      side = 3, line = 3.5, cex = 0.95, font = 2)

Nota: el Intervalo 1 [-160° , -130°) se omite del ajuste inferencial formal por concentrar apenas 37 registros (0.49% del total), un volumen insuficiente para un ajuste representativo.

6.2 Histogramas Individuales por Zona

z1 <- Variable[Variable >= zona1_inf & Variable < zona1_sup]
z2 <- Variable[Variable >= zona2_inf & Variable <= zona2_sup]

res1 <- ajustar_zona(z1, pico = "derecha")
res2 <- ajustar_zona(z2, pico = "izquierda")

cat("Zona 1 [", zona1_inf, ",", zona1_sup, "):", length(z1), "obs | Exponencial (rate =",
    round(res1$fit$estimate["rate"], 5), ")\n")
## Zona 1 [ -130 , -100 ): 2689 obs | Exponencial (rate = 0.1159 )
cat("Zona 2 [", zona2_inf, ",", zona2_sup, "]:", length(z2), "obs | Exponencial (rate =",
    round(res2$fit$estimate["rate"], 5), ")\n")
## Zona 2 [ -100 , 200 ]: 4811 obs | Exponencial (rate = 0.01304 )
plot_zona(z1, res1, paste0("Zona 1 [", round(zona1_inf, 1), " a ", round(zona1_sup, 1), "] — Exponencial"),
          "gray45", zona1_inf, zona1_sup, k_zona = 4)

plot_zona(z2, res2, paste0("Zona 2 [", round(zona2_inf, 1), " a ", round(zona2_sup, 1), "] — Exponencial"),
          "#AED6F1", zona2_inf, zona2_sup, k_zona = 7)

7.Test de bondad

7.1 Test de Pearson

fr1 <- frecuencias_zona(z1, res1, zona1_inf, zona1_sup, k_zona = 4)
fr2 <- frecuencias_zona(z2, res2, zona2_inf, zona2_sup, k_zona = 7)

r_pearson1 <- round(cor(fr1$hi_obs, fr1$hi_teo) * 100, 2)
r_pearson2 <- round(cor(fr2$hi_obs, fr2$hi_teo) * 100, 2)

cat("Pearson Zona 1 - Exponencial (%):", r_pearson1, "\n")
## Pearson Zona 1 - Exponencial (%): 97.9
cat("Pearson Zona 2 - Exponencial (%):", r_pearson2, "\n")
## Pearson Zona 2 - Exponencial (%): 79.59

7.2 Test Chi-Cuadrado

calc_chi <- function(fr) {
  chi2_calc <- sum(((fr$hi_obs - fr$hi_teo)^2) / fr$hi_teo)
  gl        <- length(fr$hi_obs) - 2
  chi2_crit <- qchisq(0.95, gl)
  list(chi2_calc = chi2_calc, chi2_crit = chi2_crit)
}

chi1 <- calc_chi(fr1); chi2 <- calc_chi(fr2)

cat("Chi-Cuadrado Zona 1:", round(chi1$chi2_calc, 4), "| Crítico:", round(chi1$chi2_crit, 4),
    "| ¿Aceptado?:", chi1$chi2_calc < chi1$chi2_crit, "\n")
## Chi-Cuadrado Zona 1: 0.0715 | Crítico: 5.9915 | ¿Aceptado?: TRUE
cat("Chi-Cuadrado Zona 2:", round(chi2$chi2_calc, 4), "| Crítico:", round(chi2$chi2_crit, 4),
    "| ¿Aceptado?:", chi2$chi2_calc < chi2$chi2_crit, "\n")
## Chi-Cuadrado Zona 2: 0.3474 | Crítico: 11.0705 | ¿Aceptado?: TRUE

Tabla resumen del test

resumen_ajuste <- data.frame(
  Zona = c(paste0("Intervalo 1 [", round(intervalo1_inf, 1), ", ", round(intervalo1_sup, 1), ") — omitido"),
           paste0("Zona 1 [", round(zona1_inf, 1), ", ", round(zona1_sup, 1), ")"),
           paste0("Zona 2 [", round(zona2_inf, 1), ", ", round(zona2_sup, 1), "]")),
  Modelo                  = c("No aplica", "Exponencial", "Exponencial"),
  `Test Pearson (%)`      = c(NA, r_pearson1, r_pearson2),
  `Chi Cuadrado`          = c(NA, round(chi1$chi2_calc, 4), round(chi2$chi2_calc, 4)),
  `Umbral de Aceptación`  = c(NA, round(chi1$chi2_crit, 4), round(chi2$chi2_crit, 4)),
  check.names = FALSE
) %>%
  mutate(Resultado = ifelse(is.na(`Chi Cuadrado`), "Omitido por irrelevancia de datos",
                      ifelse(`Chi Cuadrado` < `Umbral de Aceptación`,
                             "Modelo Aceptado", "Modelo Rechazado")))

resumen_ajuste %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N°2: Resumen del Test de Bondad — Modelo Exponencial por Zonas**"),
    subtitle = "2 zonas ajustadas con modelo Exponencial; el Intervalo 1 se omitio por irrelevancia de datos"
  ) %>%
  sub_missing(columns = everything(), missing_text = "—") %>%
  tab_source_note(source_note = "Autor: Grupo 5") %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.border.top.color            = "black",
    table.border.bottom.color         = "black",
    column_labels.border.top.color    = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    table_body.border.bottom.color    = "black"
  )
Tabla N°2: Resumen del Test de Bondad — Modelo Exponencial por Zonas
2 zonas ajustadas con modelo Exponencial; el Intervalo 1 se omitio por irrelevancia de datos
Zona Modelo Test Pearson (%) Chi Cuadrado Umbral de Aceptación Resultado
Intervalo 1 [-160, -130) — omitido No aplica — — — Omitido por irrelevancia de datos
Zona 1 [-130, -100) Exponencial 97.90 0.0715 5.9915 Modelo Aceptado
Zona 2 [-100, 200] Exponencial 79.59 0.3474 11.0705 Modelo Aceptado
Autor: Grupo 5

8.Cálculo de probabilidades

A partir del modelo Exponencial validado para la Zona 2, se calcula la probabilidad teórica de que un yacimiento se ubique en el tramo inicial, central y final de esa zona, usando los cuartiles 25% y 75% de su rango como puntos de corte.

q1_z2 <- zona2_inf + (zona2_sup - zona2_inf) * 0.25
q3_z2 <- zona2_inf + (zona2_sup - zona2_inf) * 0.75

pexp_zona2 <- function(v) pexp(v - res2$offset, rate = res2$fit$estimate["rate"])

p_a <- pexp_zona2(q1_z2)
p_b <- pexp_zona2(q3_z2) - pexp_zona2(q1_z2)
p_c <- 1 - pexp_zona2(q3_z2)

cat("P(X <", round(q1_z2, 2), ") =", round(p_a * 100, 2), "%\n")
## P(X < -25 ) = 62.88 %
cat("P(", round(q1_z2, 2), "<= X <", round(q3_z2, 2), ") =", round(p_b * 100, 2), "%\n")
## P( -25 <= X < 125 ) = 31.87 %
cat("P(X >=", round(q3_z2, 2), ") =", round(p_c * 100, 2), "%\n")
## P(X >= 125 ) = 5.25 %
cat("Registros omitidos del análisis por irrelevancia de datos (Intervalo 1):", n_omitido, "\n")
## Registros omitidos del análisis por irrelevancia de datos (Intervalo 1): 37

9.Intervalo de confianza

Por el Teorema del Límite Central (TLC), se calcula el intervalo de confianza al 95% para la media poblacional de la Longitud a nivel mundial (Z = 1.96), usando el total de registros válidos.

x_bar   <- mean(Variable, na.rm = TRUE)
sigma   <- sd(Variable, na.rm = TRUE)
z       <- qnorm(0.975)
margen  <- z * (sigma / sqrt(n))

ic_inf <- x_bar - margen
ic_sup <- x_bar + margen

cat("Media muestral (centro de masa longitudinal):", round(x_bar, 3), "\n")
## Media muestral (centro de masa longitudinal): -54.653
cat("Intervalo de confianza al 95%: (", round(ic_inf, 3), "° ,", round(ic_sup, 3), "° )\n")
## Intervalo de confianza al 95%: ( -56.186 ° , -53.12 ° )

10.Conclusión

El histograma general evidenció una concentración marcadamente decreciente en la Longitud, propia de una distribución Exponencial. El Intervalo 1 [-160°, -130°) se omitió del ajuste inferencial formal por concentrar apenas 37 registros (0.49%), un volumen insuficiente para un ajuste representativo. Con el resto de los datos se definieron 2 zonas, ambas modeladas con una distribución Exponencial: la Zona 1 (n=2689, rate = 0.1159) obtuvo una correlación de Pearson del 97.9%, y la Zona 2 (n=4811, rate = 0.01304) del 79.59%. Bajo el modelo Exponencial de la Zona 2, la probabilidad de que un yacimiento se ubique en el tramo inicial es del 62.88%, en el tramo central del 31.87%, y en el tramo final del 5.25%. Finalmente, el intervalo de confianza al 95% para la media poblacional de la Longitud a nivel mundial se ubica entre -56.186° y -53.12°.