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

2.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

3.Selecció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 principal, aunque —como se detalla en la Sección 6— presenta también un segundo repunte que exige un tratamiento en 2 zonas.

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

4.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: 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

5.Gráfica de Distribución de Frecuencia

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

Antes de proponer cualquier modelo, se visualiza también el histograma general de la Latitud, usando los mismos intervalos calculados en la Sección 4.

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.Tratamiento de Datos

El histograma general muestra un único pico de concentración muy marcado en el Hemisferio Norte, con una caída sostenida hacia la izquierda (Latitudes negativas / Hemisferio Sur) y otra hacia la derecha (Latitudes por encima de 60°). Se toma como punto de corte el Ecuador (0°): hacia la izquierda (Hemisferio Sur, [-60°, 0°)) la densidad crece levemente a medida que se acerca al Ecuador, así que se ajusta con una Exponencial; desde el Ecuador hacia la derecha (Hemisferio Norte, [0°, 80°]) se concentra el pico principal y el segundo repunte cercano a 50°-55° dentro de una misma forma acampanada, así que se ajusta con una Normal.

# Ajusta un modelo (Exponencial o Normal) a los datos de una zona.
# Exponencial:
#   pico = "izquierda": la densidad decae de izquierda a derecha (estándar).
#   pico = "derecha":   la densidad decae de derecha a izquierda (el pico
#                       está pegado al límite superior de la zona).
#   En ambos casos se desplaza la variable a dominio positivo antes de
#   ajustar (fitdistr lo exige); el desplazamiento se deshace solo para
#   graficar.
# Normal: no requiere ningún desplazamiento.
ajustar_zona <- function(datos, modelo, pico = "izquierda") {
  if (modelo == "exponential") {
    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, modelo = modelo)
  } else {
    f <- fitdistr(datos, "normal")
    list(fit = f, offset = 0, pico = pico, 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
  if (res$modelo == "exponential") {
    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)
  } else {
    hi_teo <- dnorm(mc, mean = res$fit$estimate["mean"], sd = res$fit$estimate["sd"]) * 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)
  if (res$modelo == "exponential") {
    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"])
    etiqueta <- "Curva Exponencial"
  } else {
    ys       <- dnorm(xs, mean = res$fit$estimate["mean"], sd = res$fit$estimate["sd"])
    etiqueta <- "Curva Normal"
  }

  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", 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 - 0))]
zona2_inf <- zona1_sup
zona2_sup <- cortes_int[length(cortes_int)]

cat("Zona 1: [", zona1_inf, ",", zona1_sup, ")\n")
## Zona 1: [ -60 , 0 )
cat("Zona 2: [", zona2_inf, ",", zona2_sup, "]\n")
## Zona 2: [ 0 , 80 ]
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 [ -60 , 0 ): 636 obs
cat("Zona 2 [", zona2_inf, ",", zona2_sup, "]:", length(z2), "obs\n")
## Zona 2 [ 0 , 80 ]: 6901 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 )

A continuación se visualiza el histograma con el corte de zona aplicado:

lim_min_z <- zona1_inf
lim_max_z <- zona2_sup

brks_zon <- cortes_int[cortes_int >= lim_min_z & cortes_int <= lim_max_z]
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", "#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 = zona1_sup, col = "black", lty = 2, lwd = 2)

legend("topright",
       legend = c("Zona 1 (Exponencial)", "Zona 2 (Normal)"),
       fill = c("#82E0AA", "#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)

7.Conjetura

  • Zona 1 [-60° , 0°): la densidad crece hacia el pico principal, así que se conjetura un modelo Exponencial con el pico en el extremo derecho (decae de derecha a izquierda).
  • Zona 2 [0° , 80°]: agrupa el pico principal y el segundo repunte cercano a 50°-55° dentro de una misma forma acampanada, así que se conjetura un modelo Normal.

8.Parámetros

res1 <- ajustar_zona(z1, modelo = "exponential", pico = "derecha")
res2 <- ajustar_zona(z2, modelo = "normal")

cat("Zona 1 | Exponencial (rate =", round(res1$fit$estimate["rate"], 5), ")\n")
## Zona 1 | Exponencial (rate = 0.04464 )
cat("Zona 2 | Normal (media =", round(res2$fit$estimate["mean"], 3),
    ", sd =", round(res2$fit$estimate["sd"], 3), ")\n")
## Zona 2 | Normal (media = 37.199 , sd = 15.966 )

9.Sobreposición de la Realidad con el Modelo

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), "] — Normal"),
          "#AED6F1", zona2_inf, zona2_sup, k_zona = 4)

10.Test de Bondad

10.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 - Exponencial (%):", r_pearson1, "\n")
## Pearson Zona 1 - Exponencial (%): 96.98
cat("Pearson Zona 2 - Normal (%):", r_pearson2, "\n")
## Pearson Zona 2 - Normal (%): 99.19

10.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 = 1)  # Exponencial: rate
chi2 <- calc_chi(fr2, n_param = 2)  # Normal: media y sd

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.1055 | 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.0106 | Crítico: 3.8415 | ¿Aceptado?: TRUE

A continuación se presenta la 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("Exponencial", "Normal"),
  `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: Exponencial | Zona 2: Normal"
  ) %>%
  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: Exponencial | Zona 2: Normal
Zona Modelo Test Pearson (%) Chi Cuadrado Umbral de Aceptación Resultado
Zona 1 [-60, 0) Exponencial 96.98 0.1055 5.9915 Modelo Aceptado
Zona 2 [0, 80] Normal 99.19 0.0106 3.8415 Modelo Aceptado
Autor: Grupo 5

11.Cálculo de Probabilidades

A partir del modelo Normal 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

pnorm_zona2 <- function(v) pnorm(v, mean = res2$fit$estimate["mean"], sd = res2$fit$estimate["sd"])

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

cat("P(X <", round(q1_z2, 2), ") =", round(p_a * 100, 2), "%\n")
## P(X < 20 ) = 14.07 %
cat("P(", round(q1_z2, 2), "<= X <", round(q3_z2, 2), ") =", round(p_b * 100, 2), "%\n")
## P( 20 <= X < 60 ) = 78.27 %
cat("P(X >=", round(q3_z2, 2), ") =", round(p_c * 100, 2), "%\n")
## P(X >= 60 ) = 7.66 %

12.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("Desviación estándar:", round(sigma, 3), "\n")
## Desviación estándar: 22.829
cat("Intervalo de confianza al 95%: (", round(ic_inf, 3), "° ,", round(ic_sup, 3), "° )\n")
## Intervalo de confianza al 95%: ( 31.738 ° , 32.769 ° )

13.Conclusión

El comportamiento de la Latitud se explica con un modelo Exponencial en la Zona 1 [-60°, 0°) de parámetro rate = 0.04464 (n = 636), y con un modelo Normal en la Zona 2 [0°, 80°] de parámetros media = 37.199° y sd = 15.966° (n = 6901), con correlaciones de Pearson de 96.98% y 99.19% respectivamente, ambas dentro del umbral aceptado por el Test Chi-Cuadrado (Sección 10).

Podemos afirmar con un 95% de confianza que la media aritmética real de la Latitud se encuentra entre 31.738° y 32.769°, con una desviación estándar de 22.829°.

Adicionalmente, bajo el modelo Normal de la Zona 2, la probabilidad de que un yacimiento se ubique en el tramo inicial es del 14.07%, en el tramo central del 78.27%, y en el tramo final del 7.66%.