0.Carga de 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)

1.Carga de Datos

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.Selecció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, el modelo teórico a ajustar es una distribución Log-normal, aplicada de forma independiente a cada hemisferio.

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 Frecuencia

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.Gráfica de Distribución de Frecuencia

Con la Tabla N°1 ya construida, se grafica la densidad de probabilidad de cada intervalo para visualizar la forma general de la distribución. Esto es el paso previo indispensable antes de conjeturar qué modelo teórico podría ajustarse a los datos.

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"))

Se visualiza también el histograma general en formato base R, con los mismos intervalos de la Sección 3, para ubicar visualmente los picos de concentración.

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)

Se confirma un pico dominante en torno a los -100° (Hemisferio Occidental) y un segundo repunte en torno a los 0°-20° (Hemisferio Oriental), por lo que se trata la variable en 2 zonas independientes, como se desarrolla a continuación.

5.Modelo aplicado por Zona

Se reutilizan los cortes de la Regla de Sturges (Sección 3), agrupados en 2 zonas: Zona 1 (Hemisferio Occidental) y Zona 2 (Hemisferio Oriental).

# Ajusta un modelo (Log-normal o Exponencial) a los datos de una zona.
# Se desplaza la variable a dominio positivo (fitdistr exige soporte
# positivo); el desplazamiento se deshace solo para graficar.
ajustar_zona <- function(datos, modelo) {
  offset <- floor(min(datos)) - 1
  d_pos  <- datos - offset
  f <- fitdistr(d_pos, modelo)
  list(fit = f, offset = offset, modelo = modelo)
}

# 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 <- mc - res$offset
  mc_pos[mc_pos <= 0] <- 1e-9
  hi_teo <- if (res$modelo == "lognormal") {
    dlnorm(mc_pos, meanlog = res$fit$estimate["meanlog"],
           sdlog = res$fit$estimate["sdlog"]) * diff(brks)
  } else {
    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 teórica 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 <- xs - res$offset
  xs_pos[xs_pos <= 0] <- 1e-9
  if (res$modelo == "lognormal") {
    ys      <- dlnorm(xs_pos, meanlog = res$fit$estimate["meanlog"],
                       sdlog = res$fit$estimate["sdlog"])
    etiqueta <- "Curva Log-normal"
  } else {
    ys      <- dexp(xs_pos, rate = res$fit$estimate["rate"])
    etiqueta <- "Curva Exponencial"
  }

  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", etiqueta),
         fill = c(col_fill, NA), border = c("black", NA),
         lty = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.8)
}
zona1_inf <- cortes_int[1]
zona1_sup <- cortes_int[which.min(abs(cortes_int - (-10)))]
zona2_inf <- zona1_sup
zona2_sup <- cortes_int[length(cortes_int)]

cat("Zona 1 — Hemisferio Occidental [", zona1_inf, ",", zona1_sup, "):",
    sum(Variable >= zona1_inf & Variable < zona1_sup), "obs\n")
## Zona 1 — Hemisferio Occidental [ -160 , -10 ): 5198 obs
cat("Zona 2 — Hemisferio Oriental   [", zona2_inf, ",", zona2_sup, "]:",
    sum(Variable >= zona2_inf & Variable <= zona2_sup), "obs\n")
## Zona 2 — Hemisferio Oriental   [ -10 , 200 ]: 2339 obs
z1 <- Variable[Variable >= zona1_inf & Variable < zona1_sup]
z2 <- Variable[Variable >= zona2_inf & Variable <= zona2_sup]

cat("Zona 1 [", zona1_inf, ",", zona1_sup, "):", length(z1), "obs\n")
## Zona 1 [ -160 , -10 ): 5198 obs
cat("Zona 2 [", zona2_inf, ",", zona2_sup, "]:", length(z2), "obs\n")
## Zona 2 [ -10 , 200 ]: 2339 obs
cat("Verificación — n(Zona 1) + n(Zona 2):", length(z1) + length(z2), "(debe ser", n, ")\n")
## Verificación — n(Zona 1) + n(Zona 2): 7537 (debe ser 7537 )
h_zon <- hist(Variable, breaks = cortes_int, plot = FALSE)
h_zon$density <- conteo$hi / Amplitud_int

mc_zon <- (head(cortes_int, -1) + tail(cortes_int, -1)) / 2
colores_zona <- ifelse(mc_zon < zona2_inf, "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(cortes_int, 0), labels = round(cortes_int, 0), las = 2, cex.axis = 0.75)
abline(v = zona2_inf, col = "black", lty = 2, lwd = 2)

legend("topright",
       legend = c("Zona 1 — Log-normal", "Zona 2 — Exponencial"),
       fill = c("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)

6.Conjetura

  • Zona 1 — Hemisferio Occidental [-160° , -10°): la densidad sube hasta el pico cercano a -100° y luego decae, así que se conjetura un modelo Log-normal.
  • Zona 2 — Hemisferio Oriental [-10° , 200°]: la densidad decae de forma natural desde su extremo izquierdo, así que se conjetura un modelo Exponencial.

Se plantea entonces, para cada zona, la siguiente hipótesis:

  • H₀: la Longitud dentro de la zona se distribuye según el modelo conjeturado.
  • H₁: la Longitud dentro de la zona no se distribuye según el modelo conjeturado.

Esta hipótesis se contrasta formalmente en la Sección 9 (Test de Bondad).

7.Parámetros

res1 <- ajustar_zona(z1, modelo = "lognormal")
res2 <- ajustar_zona(z2, modelo = "exponential")

cat("Zona 1 | Log-normal (meanlog =", round(res1$fit$estimate["meanlog"], 5),
    ", sdlog =", round(res1$fit$estimate["sdlog"], 5), ")\n")
## Zona 1 | Log-normal (meanlog = 4.00067 , sdlog = 0.35003 )
cat("Zona 2 | Exponencial (rate =", round(res2$fit$estimate["rate"], 5), ")\n")
## Zona 2 | Exponencial (rate = 0.02063 )

8.Sobreposición de la Realidad con el Modelo

plot_zona(z1, res1, paste0("Zona 1 [", round(zona1_inf, 1), " a ", round(zona1_sup, 1), "] — Log-normal"),
          "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 = 4)

9.Test de Bondad

9.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)

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 - Log-normal (%):", r_pearson1, "\n")
## Pearson Zona 1 - Log-normal (%): 99.96
cat("Pearson Zona 2 - Exponencial (%):", r_pearson2, "\n")
## Pearson Zona 2 - Exponencial (%): 98.86

9.2 Test Chi-Cuadrado

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

chi1 <- calc_chi(fr1, n_param = 2)  # Log-normal: meanlog y sdlog
chi2 <- calc_chi(fr2, n_param = 1)  # Exponencial: rate

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.1397 | Crítico: 3.8415 | ¿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.071 | 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), "]")),
  Modelo                  = c("Log-normal", "Exponencial"),
  `Test Pearson (%)`      = c(r_pearson1, r_pearson2),
  `Chi Cuadrado`          = c(round(chi1$chi2_calc, 4), round(chi2$chi2_calc, 4)),
  `Umbral de Aceptación`  = c(round(chi1$chi2_crit, 4), round(chi2$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 por Zona**"),
    subtitle = "Zona 1: Log-normal | Zona 2: Exponencial"
  ) %>%
  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 por Zona
Zona 1: Log-normal | Zona 2: Exponencial
Zona Modelo Test Pearson (%) Chi Cuadrado Umbral de Aceptación Resultado
Zona 1 [-160, -10) Log-normal 99.96 0.1397 3.8415 Modelo Aceptado
Zona 2 [-10, 200] Exponencial 98.86 0.0710 5.9915 Modelo Aceptado
Autor: Grupo 5

10.Cálculo de Probabilidades

Bajo el modelo Exponencial de 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.

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 < 42.5 ) = 66.14 %
cat("P(", round(q1_z2, 2), "<= X <", round(q3_z2, 2), ") =", round(p_b * 100, 2), "%\n")
## P( 42.5 <= X < 147.5 ) = 29.98 %
cat("P(X >=", round(q3_z2, 2), ") =", round(p_c * 100, 2), "%\n")
## P(X >= 147.5 ) = 3.88 %

11.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 (Z = 1.96).

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:", round(x_bar, 3), "\n")
## Media muestral: -54.653
cat("Desviación estándar:", round(sigma, 3), "\n")
## Desviación estándar: 67.9
cat("Intervalo de confianza al 95%: (", round(ic_inf, 3), "° ,", round(ic_sup, 3), "° )\n")
## Intervalo de confianza al 95%: ( -56.186 ° , -53.12 ° )

12.Conclusión

El comportamiento de la Longitud se explica con un modelo Log-normal en la Zona 1 — Hemisferio Occidental (meanlog = 4.00067, sdlog = 0.35003) y un modelo Exponencial en la Zona 2 — Hemisferio Oriental (rate = 0.02063). Podemos afirmar con un 95% de confianza que la media aritmética real de la Longitud se encuentra entre -56.186° y -53.12°, con una desviación estándar de 67.9°.