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 Latitude indica la coordenada geográfica de latitud de cada yacimiento. Es una variable cuantitativa continua, cuyos valores van desde coordenadas negativas (Hemisferio Sur) hasta positivas (Hemisferio Norte). A diferencia de la Longitud (que mostró dos picos de concentración), la Latitud concentra sus registros en un único pico, por lo que su tratamiento analítico se hace en 2 zonas en vez de 3.

n_total <- nrow(datos)
Variable <- na.omit(as.numeric(datos$Latitude))
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: -53.971
cat("Valor máximo:", round(max(Variable), 3), "\n")
## Valor máximo: 73.434
n_sur   <- sum(Variable <= 0)
n_norte <- sum(Variable > 0)

cat("Registros en Hemisferio Sur (<= 0°):", n_sur,
    "(", round(n_sur / n * 100, 2), "%)\n")
## Registros en Hemisferio Sur (<= 0°): 636 ( 8.44 %)
cat("Registros en Hemisferio Norte (> 0°):", n_norte,
    "(", round(n_norte / n * 100, 2), "%)\n")
## Registros en Hemisferio Norte (> 0°): 6901 ( 91.56 %)

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: 7
cat("Amplitud de clase:", Amplitud_int, "\n")
## Amplitud de clase: 20
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: [ 20 , 40 ) con 2971 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 Latitud 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 Latitud Geográfica
Lím. Inf (°) Lím. Sup (°) Marca de Clase Frecuencia (ni) Proporción (hi) Porcentaje (hi%)
-60 -40 −50.0 79 0.0105 1.05
-40 -20 −30.0 248 0.0329 3.29
-20 0 −10.0 309 0.0410 4.10
0 20 10.0 956 0.1268 12.68
20 40 30.0 2971 0.3942 39.42
40 60 50.0 2666 0.3537 35.37
60 80 70.0 308 0.0409 4.09
TOTAL TOTAL TOTAL 7537 1.0000 100.00
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 Latitud",
    x       = "Latitud (°)",
    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 Latitud, 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("Latitud (grados)", side = 1, line = 3.5, cex = 0.9)
mtext("Histograma General — Latitud 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 en torno a los [20°, 40°), con una caída marcada hacia la izquierda (Latitudes negativas / Hemisferio Sur). Sin embargo, al examinar más de cerca el tramo derecho, aparece un segundo repunte real de concentración entre 50° y 55° (probablemente por clústeres geográficos como Rusia y el Mar del Norte) — el mismo patrón bimodal que se observó en Longitud. Por eso, igual que allá, se necesitan 3 zonas, usando el rango completo de la variable sin omitir ningún tramo:

zona1_inf      <- cortes_int[1]
zona1_sup      <- cortes_int[which.min(abs(cortes_int - 40))]
zona2_inf      <- zona1_sup
zona2_sup      <- 47
zona3_inf      <- 47
zona3_sup      <- cortes_int[which.min(abs(cortes_int - 80))]
  • Zona 1 [-60° , 40°): la densidad crece hacia el primer pico general, así que se ajusta un modelo Exponencial con el pico en el extremo derecho (decae de derecha a izquierda).
  • Zona 2 [40° , 47°): justo después del primer pico, la densidad decae desde su extremo izquierdo hasta el valle en 47°, así que se ajusta un modelo Exponencial estándar (decae de izquierda a derecha).
  • Zona 3 [47° , 80°]: tras el segundo repunte en torno a 50°-55°, la densidad vuelve a decaer desde ese extremo izquierdo, así que también se ajusta un modelo Exponencial estándar (decae de izquierda a derecha).
# 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("Latitud (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 <- zona1_inf
lim_max_z <- zona3_sup

brks_zon <- sort(unique(c(cortes_int[cortes_int >= lim_min_z & cortes_int <= lim_max_z], zona2_sup)))
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_sup, "#82E0AA",
                 ifelse(mc_zon < zona2_sup, "#F5B041", "#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, zona2_sup, zona3_sup), col = "black", lty = 2, lwd = 2)

legend("topright",
       legend = c("Zona 1 (Exponencial)", "Zona 2 (Exponencial)", "Zona 3 (Exponencial)"),
       fill = c("#82E0AA", "#F5B041", "#AED6F1"), border = "black",
       bty = "n", cex = 0.8)

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

6.2 Histogramas Individuales por Zona

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

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

cat("Zona 1 [", zona1_inf, ",", zona1_sup, "):", length(z1), "obs | Exponencial (rate =",
    round(res1$fit$estimate["rate"], 5), ")\n")
## Zona 1 [ -60 , 40 ): 4563 obs | Exponencial (rate = 0.04571 )
cat("Zona 2 [", zona2_inf, ",", zona2_sup, "):", length(z2), "obs | Exponencial (rate =",
    round(res2$fit$estimate["rate"], 5), ")\n")
## Zona 2 [ 40 , 47 ): 429 obs | Exponencial (rate = 0.26689 )
cat("Zona 3 [", zona3_inf, ",", zona3_sup, "]:", length(z3), "obs | Exponencial (rate =",
    round(res3$fit$estimate["rate"], 5), ")\n")
## Zona 3 [ 47 , 80 ]: 2545 obs | Exponencial (rate = 0.12458 )
plot_zona(z1, res1, paste0("Zona 1 [", round(zona1_inf, 1), " a ", round(zona1_sup, 1), "] — Exponencial"),
          "#82E0AA", zona1_inf, zona1_sup, k_zona = 4)

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

plot_zona(z3, res3, paste0("Zona 3 [", round(zona3_inf, 1), " a ", round(zona3_sup, 1), "] — Exponencial"),
          "#AED6F1", zona3_inf, zona3_sup, k_zona = 4)

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 = 4)
fr3 <- frecuencias_zona(z3, res3, zona3_inf, zona3_sup, k_zona = 4)

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

cat("Pearson Zona 1 - Exponencial (%):", r_pearson1, "\n")
## Pearson Zona 1 - Exponencial (%): 99.76
cat("Pearson Zona 2 - Exponencial (%):", r_pearson2, "\n")
## Pearson Zona 2 - Exponencial (%): 97.52
cat("Pearson Zona 3 - Exponencial (%):", r_pearson3, "\n")
## Pearson Zona 3 - Exponencial (%): 99.43

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); chi3 <- calc_chi(fr3)

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.0339 | 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.0223 | Crítico: 5.9915 | ¿Aceptado?: TRUE
cat("Chi-Cuadrado Zona 3:", round(chi3$chi2_calc, 4), "| Crítico:", round(chi3$chi2_crit, 4),
    "| ¿Aceptado?:", chi3$chi2_calc < chi3$chi2_crit, "\n")
## Chi-Cuadrado Zona 3: 0.0454 | Crítico: 5.9915 | ¿Aceptado?: TRUE

Tabla resumen del test

resumen_ajuste <- data.frame(
  Zona = c(paste0("Zona 1 [", round(zona1_inf, 1), ", ", round(zona1_sup, 1), ")"),
           paste0("Zona 2 [", round(zona2_inf, 1), ", ", round(zona2_sup, 1), ")"),
           paste0("Zona 3 [", round(zona3_inf, 1), ", ", round(zona3_sup, 1), "]")),
  Modelo                  = c("Exponencial", "Exponencial", "Exponencial"),
  `Test Pearson (%)`      = c(r_pearson1, r_pearson2, r_pearson3),
  `Chi Cuadrado`          = c(round(chi1$chi2_calc, 4), round(chi2$chi2_calc, 4), round(chi3$chi2_calc, 4)),
  `Umbral de Aceptación`  = c(round(chi1$chi2_crit, 4), round(chi2$chi2_crit, 4), round(chi3$chi2_crit, 4)),
  check.names = FALSE
) %>%
  mutate(Resultado = 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 = "3 zonas ajustadas con modelo Exponencial (rango completo de la variable, ver Seccion 6)"
  ) %>%
  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
3 zonas ajustadas con modelo Exponencial (rango completo de la variable, ver Seccion 6)
Zona Modelo Test Pearson (%) Chi Cuadrado Umbral de Aceptación Resultado
Zona 1 [-60, 40) Exponencial 99.76 0.0339 5.9915 Modelo Aceptado
Zona 2 [40, 47) Exponencial 97.52 0.0223 5.9915 Modelo Aceptado
Zona 3 [47, 80] Exponencial 99.43 0.0454 5.9915 Modelo Aceptado
Autor: Grupo 5

8.Cálculo de probabilidades

A partir del modelo Exponencial validado para la Zona 3, 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_z3 <- zona3_inf + (zona3_sup - zona3_inf) * 0.25
q3_z3 <- zona3_inf + (zona3_sup - zona3_inf) * 0.75

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

p_a <- pexp_zona3(q1_z3)
p_b <- pexp_zona3(q3_z3) - pexp_zona3(q1_z3)
p_c <- 1 - pexp_zona3(q3_z3)

cat("P(X <", round(q1_z3, 2), ") =", round(p_a * 100, 2), "%\n")
## P(X < 55.25 ) = 68.41 %
cat("P(", round(q1_z3, 2), "<= X <", round(q3_z3, 2), ") =", round(p_b * 100, 2), "%\n")
## P( 55.25 <= X < 71.75 ) = 27.55 %
cat("P(X >=", round(q3_z3, 2), ") =", round(p_c * 100, 2), "%\n")
## P(X >= 71.75 ) = 4.04 %

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 Latitud 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 latitudinal):", round(x_bar, 3), "\n")
## Media muestral (centro de masa latitudinal): 32.254
cat("Intervalo de confianza al 95%: (", round(ic_inf, 3), "° ,", round(ic_sup, 3), "° )\n")
## Intervalo de confianza al 95%: ( 31.738 ° , 32.769 ° )

10.Conclusión

El histograma general evidenció una concentración marcadamente decreciente en la Latitud, propia de una distribución Exponencial. Al igual que en Longitud, se detectó un segundo repunte real de concentración (entre 50° y 55°), por lo que se aprovechó el rango completo de la variable (-60° a 80°), sin omitir ningún tramo, dividido en 3 zonas, todas modeladas con una distribución Exponencial: la Zona 1 (n=4563, rate = 0.04571) obtuvo una correlación de Pearson del 99.76%, la Zona 2 (n=429, rate = 0.26689) del 97.52%, y la Zona 3 (n=2545, rate = 0.12458) del 99.43%. Bajo el modelo Exponencial de la Zona 3, la probabilidad de que un yacimiento se ubique en el tramo inicial es del 68.41%, en el tramo central del 27.55%, y en el tramo final del 4.04%. Finalmente, el intervalo de confianza al 95% para la media poblacional de la Latitud a nivel mundial se ubica entre 31.738° y 32.769°.