1 Carga de Librerías

Se cargan las librerías necesarias para leer el conjunto de datos, construir la tabla de distribución de frecuencias, elaborar diagramas de barras en escala de grises y aplicar las pruebas de bondad de ajuste.

library(readr)
library(dplyr)
library(knitr)
library(kableExtra)
library(ggplot2)

2 Carga de Datos

Se carga el conjunto de datos de arrendamientos de petróleo y gas. El código permite leer archivos separados por coma o por punto y coma.

rutas_posibles <- c(
  "oil_and_gas_leases_data (2)(2).csv",
  "oil_and_gas_leases_data (2)(1).csv",
  "oil_and_gas_leases_data (2).csv",
  "oil_and_gas_leases_data.csv",
  "kansas.csv"
)

ruta_csv <- rutas_posibles[file.exists(rutas_posibles)][1]
if (is.na(ruta_csv)) {
  stop("No se encontró el archivo CSV. Coloca el archivo en la misma carpeta del Rmd.")
}

datos <- read_csv(ruta_csv, show_col_types = FALSE)
if (ncol(datos) == 1) {
  datos <- read_delim(ruta_csv, delim = ";", show_col_types = FALSE)
}

cat("Archivo cargado:", ruta_csv, "\n")
## Archivo cargado: oil_and_gas_leases_data (2).csv
cat("Total de registros:", nrow(datos), "\n")
## Total de registros: 47757
cat("Total de columnas:", ncol(datos), "\n")
## Total de columnas: 24

3 Conteo

Se extrae la variable YEARS_ACTIVE, que representa los años activos del pozo o arrendamiento. Como esta variable toma valores enteros positivos, se trabaja como variable cuantitativa discreta.

Para mantener una muestra manejable y comparable con los demás modelos, se toma una muestra de tamaño 80. Además, se aplica la regla de Sturges para definir el número de intervalos, respetando un máximo de 10 intervalos.

# Buscar la columna de años activos aunque el nombre tenga pequeñas variaciones
nombres <- names(datos)
col_years <- nombres[tolower(nombres) %in% c("years_active", "year_active", "active_years")]

if (length(col_years) == 0) {
  col_years <- nombres[grepl("year", nombres, ignore.case = TRUE) & grepl("active", nombres, ignore.case = TRUE)]
}

if (length(col_years) == 0) {
  stop("No se encontró la columna YEARS_ACTIVE. Revisa el nombre exacto en el CSV.")
}

poblacion <- datos %>%
  mutate(YEARS_ANALISIS = suppressWarnings(as.integer(.data[[col_years[1]]]))) %>%
  filter(!is.na(YEARS_ANALISIS), YEARS_ANALISIS >= 1) %>%
  pull(YEARS_ANALISIS)

n_muestra <- 80
set.seed(1)
years_muestra <- sample(poblacion, size = n_muestra, replace = FALSE)

n <- length(years_muestra)
k_sturges <- ceiling(1 + 3.322 * log10(n))
k <- min(k_sturges, 10)

cat("Variable utilizada:", col_years[1], "\n")
## Variable utilizada: YEARS_ACTIVE
cat("Tamaño de muestra:", n, "\n")
## Tamaño de muestra: 80
cat("Número de intervalos por Sturges:", k_sturges, "\n")
## Número de intervalos por Sturges: 8
cat("Número de intervalos usado:", k, "\n")
## Número de intervalos usado: 8
cat("Mínimo:", min(years_muestra), "\n")
## Mínimo: 1
cat("Máximo:", max(years_muestra), "\n")
## Máximo: 82

4 Tabla de Distribución de Frecuencias

Se construye una sola tabla de distribución de frecuencias para la variable YEARS_ACTIVE. La tabla incluye frecuencia absoluta, frecuencia relativa, frecuencia acumulada y frecuencia relativa acumulada. Al final se añade la fila de total.

# Intervalos usando Sturges, con máximo de 10 intervalos
min_x <- min(years_muestra)
max_x <- max(years_muestra)
amplitud <- ceiling((max_x - min_x + 1) / k)

limites_inf <- seq(min_x, by = amplitud, length.out = k)
limites_sup <- c(limites_inf[-1] - 1, max_x)
limites_sup[length(limites_sup)] <- max_x

intervalos <- data.frame(
  li = limites_inf,
  ls = limites_sup
) %>%
  mutate(
    MC = round((li + ls) / 2, 2),
    Intervalo = paste0("[", li, " - ", ls, "]")
  )

# Asignación a intervalos
contar_intervalos <- function(x, tabla_intervalos) {
  sapply(1:nrow(tabla_intervalos), function(i) {
    sum(x >= tabla_intervalos$li[i] & x <= tabla_intervalos$ls[i])
  })
}

TDF <- intervalos %>%
  mutate(
    ni = as.numeric(contar_intervalos(years_muestra, intervalos)),
    hi = round(100 * ni / sum(ni), 2),
    Ni = cumsum(ni),
    Hi = round(cumsum(hi), 2)
  ) %>%
  select(Intervalo, MC, ni, hi, Ni, Hi)

TDF$Intervalo <- factor(TDF$Intervalo, levels = TDF$Intervalo)

TDF_total <- bind_rows(
  TDF %>% mutate(Intervalo = as.character(Intervalo)),
  data.frame(
    Intervalo = "TOTAL",
    MC = NA,
    ni = sum(TDF$ni),
    hi = round(sum(TDF$hi), 2),
    Ni = NA,
    Hi = 100
  )
)

kable(TDF_total, caption = "Tabla N°1: Distribución de Frecuencias — Years Active") %>%
  kable_styling(full_width = FALSE, position = "center")
Tabla N°1: Distribución de Frecuencias — Years Active
Intervalo MC ni hi Ni Hi
[1 - 11] 6 20 25.00 20 25.00
[12 - 22] 17 26 32.50 46 57.50
[23 - 33] 28 13 16.25 59 73.75
[34 - 44] 39 11 13.75 70 87.50
[45 - 55] 50 6 7.50 76 95.00
[56 - 66] 61 1 1.25 77 96.25
[67 - 77] 72 2 2.50 79 98.75
[78 - 82] 80 1 1.25 80 100.00
TOTAL NA 80 100.00 NA 100.00

5 Gráfica

Se presenta el diagrama de barras de las frecuencias observadas. La variable se mantiene como cuantitativa discreta, aunque se encuentre agrupada en intervalos; por ello se usa diagrama de barras y no histograma.

ggplot(TDF, aes(x = Intervalo, y = hi)) +
  geom_col(fill = "gray60", color = "black") +
  labs(
    title = "Gráfica N°1: Frecuencias observadas de Years Active",
    x = "Intervalo",
    y = "hi (%)"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

6 Conjetura

El diagrama de barras muestra que la variable YEARS_ACTIVE no se ajusta de forma uniforme a un único patrón geométrico en todos sus intervalos. Por esta razón, el análisis se realiza por secciones, identificando subconjuntos donde las frecuencias observadas presentan una disminución progresiva compatible con el modelo geométrico.

Se omiten los intervalos extremos o irregulares únicamente para el ajuste del modelo, no para la construcción de la tabla general. Esta omisión se justifica porque dichos intervalos representan comportamientos atípicos o de baja frecuencia que distorsionan la tendencia decreciente que se busca modelar.

Se evalúan dos secciones:

  • Sección A: desde el segundo hasta el sexto intervalo de la tabla general, es decir, desde [12–22] hasta [56–66]. Esta sección representa el tramo principal donde las frecuencias tienden a disminuir después del primer pico.
  • Sección B: desde el séptimo hasta el octavo intervalo de la tabla general, es decir, desde [67–77] hasta [78–82]. Esta sección corresponde al extremo final de la variable y se analiza por separado por su baja frecuencia.

El primer intervalo [1–11] se omite en el ajuste porque concentra pozos relativamente recientes y rompe la tendencia decreciente posterior. Los intervalos extremos se trabajan en una sección independiente para no mezclar comportamientos distintos dentro de un solo ajuste geométrico. En cada sección se calcula el parámetro del modelo geométrico y se aplican tanto el test de Pearson como la prueba Chi-cuadrado.

7 Cálculo de Parámetros y Probabilidades

Para cada sección se estima el parámetro de la distribución geométrica mediante:

\[ \hat{p}=\frac{1}{\bar{x}} \]

donde \(x\) representa el número de clase dentro de cada sección.

# Función para ajustar el modelo geométrico a una sección de la TDF
ajustar_geometrico_seccion <- function(TDF_base, filas, nombre_seccion) {
  seccion <- TDF_base[filas, ] %>%
    mutate(
      clase = seq_len(n()),
      Intervalo = as.character(Intervalo)
    )
  
  n_sec <- sum(seccion$ni)
  media_clase <- sum(seccion$clase * seccion$ni) / n_sec
  p_hat <- 1 / media_clase
  
  prob_teorica <- dgeom(seccion$clase - 1, prob = p_hat)
  prob_teorica <- prob_teorica / sum(prob_teorica)
  frecuencia_esperada <- n_sec * prob_teorica
  porcentaje_esperado <- 100 * frecuencia_esperada / sum(frecuencia_esperada)
  
  tabla_modelo <- seccion %>%
    mutate(
      Probabilidad_Teorica = round(prob_teorica, 4),
      Frecuencia_Esperada = round(frecuencia_esperada, 4),
      Porcentaje_Esperado = round(porcentaje_esperado, 2)
    ) %>%
    select(Intervalo, clase, ni, hi, Probabilidad_Teorica, Frecuencia_Esperada, Porcentaje_Esperado)
  
  # Pearson entre frecuencias observadas y esperadas
  pearson <- ifelse(
    sd(tabla_modelo$ni) > 0 && sd(tabla_modelo$Frecuencia_Esperada) > 0,
    cor(tabla_modelo$ni, tabla_modelo$Frecuencia_Esperada),
    NA
  )
  
  # Para Chi-cuadrado se fusionan clases con esperadas menores a 5
  O <- tabla_modelo$ni
  E <- tabla_modelo$Frecuencia_Esperada
  
  while (any(E < 5) && length(E) > 2) {
    idx <- which.min(E)
    if (idx == 1) {
      O[2] <- O[2] + O[1]
      E[2] <- E[2] + E[1]
      O <- O[-1]
      E <- E[-1]
    } else {
      O[idx - 1] <- O[idx - 1] + O[idx]
      E[idx - 1] <- E[idx - 1] + E[idx]
      O <- O[-idx]
      E <- E[-idx]
    }
  }
  
  # Grados de libertad para bondad de ajuste: k - 1 - parametros estimados.
  # Si una sección queda con solo dos clases, el ajuste formal tiene gl = 0;
  # para poder reportar el contraste solicitado se deja gl = 1 de forma referencial.
  gl_formal <- length(O) - 1 - 1
  gl <- ifelse(gl_formal < 1, 1, gl_formal)
  nota_gl <- ifelse(
    gl_formal < 1,
    "Chi-cuadrado referencial: la sección tiene pocas clases para una prueba formal con parámetro estimado.",
    "Chi-cuadrado formal aplicable."
  )
  
  chi_calculado <- sum((O - E)^2 / E)
  chi_critico <- qchisq(0.95, df = gl)
  p_valor <- pchisq(chi_calculado, df = gl, lower.tail = FALSE)
  decision <- ifelse(
    p_valor > 0.05 && chi_calculado < chi_critico,
    "No se rechaza H0: modelo aceptado",
    "Se rechaza H0: modelo no aceptado"
  )
  
  list(
    nombre = nombre_seccion,
    tabla = tabla_modelo,
    p_hat = p_hat,
    pearson = pearson,
    chi_calculado = chi_calculado,
    chi_critico = chi_critico,
    gl = gl,
    p_valor = p_valor,
    decision = decision,
    O_fusionado = O,
    E_fusionado = E,
    nota_gl = nota_gl
  )
}

Sección A: Ajuste geométrico

# Sección A: desde el segundo hasta el sexto intervalo de la TDF general
# [12-22], [23-33], [34-44], [45-55], [56-66]
sec_A <- ajustar_geometrico_seccion(
  TDF_base = TDF,
  filas = 2:6,
  nombre_seccion = "Sección A"
)

T_A_total <- bind_rows(
  sec_A$tabla,
  data.frame(
    Intervalo = "TOTAL",
    clase = NA,
    ni = sum(sec_A$tabla$ni),
    hi = round(sum(sec_A$tabla$hi), 2),
    Probabilidad_Teorica = round(sum(sec_A$tabla$Probabilidad_Teorica), 4),
    Frecuencia_Esperada = round(sum(sec_A$tabla$Frecuencia_Esperada), 4),
    Porcentaje_Esperado = 100
  )
)

kable(T_A_total, caption = "Tabla N°2: Ajuste Geométrico — Sección A") %>%
  kable_styling(full_width = FALSE, position = "center")
Tabla N°2: Ajuste Geométrico — Sección A
Intervalo clase ni hi Probabilidad_Teorica Frecuencia_Esperada Porcentaje_Esperado
[12 - 22] 1 26 32.50 0.5161 29.4194 51.61
[23 - 33] 2 13 16.25 0.2581 14.7097 25.81
[34 - 44] 3 11 13.75 0.1290 7.3548 12.90
[45 - 55] 4 6 7.50 0.0645 3.6774 6.45
[56 - 66] 5 1 1.25 0.0323 1.8387 3.23
TOTAL NA 57 71.25 1.0000 57.0000 100.00 
ggplot(sec_A$tabla, aes(x = Intervalo, y = hi)) +
  geom_col(fill = "gray60", color = "black") +
  labs(
    title = "Gráfica N°2: Frecuencias observadas — Sección A",
    x = "Intervalo",
    y = "hi (%)"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(sec_A$tabla, aes(x = Intervalo, y = Porcentaje_Esperado)) +
  geom_col(fill = "gray35", color = "black") +
  labs(
    title = "Gráfica N°3: Frecuencias esperadas — Modelo Geométrico Sección A",
    x = "Intervalo",
    y = "Porcentaje esperado (%)"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Sección B: Ajuste geométrico

# Sección B: intervalos extremos finales
# [67-77] y [78-82]
sec_B <- ajustar_geometrico_seccion(
  TDF_base = TDF,
  filas = 7:8,
  nombre_seccion = "Sección B"
)

T_B_total <- bind_rows(
  sec_B$tabla,
  data.frame(
    Intervalo = "TOTAL",
    clase = NA,
    ni = sum(sec_B$tabla$ni),
    hi = round(sum(sec_B$tabla$hi), 2),
    Probabilidad_Teorica = round(sum(sec_B$tabla$Probabilidad_Teorica), 4),
    Frecuencia_Esperada = round(sum(sec_B$tabla$Frecuencia_Esperada), 4),
    Porcentaje_Esperado = 100
  )
)

kable(T_B_total, caption = "Tabla N°3: Ajuste Geométrico — Sección B") %>%
  kable_styling(full_width = FALSE, position = "center")
Tabla N°3: Ajuste Geométrico — Sección B
Intervalo clase ni hi Probabilidad_Teorica Frecuencia_Esperada Porcentaje_Esperado
[67 - 77] 1 2 2.50 0.8 2.4 80
[78 - 82] 2 1 1.25 0.2 0.6 20
TOTAL NA 3 3.75 1.0 3.0 100 
ggplot(sec_B$tabla, aes(x = Intervalo, y = hi)) +
  geom_col(fill = "gray60", color = "black") +
  labs(
    title = "Gráfica N°4: Frecuencias observadas — Sección B",
    x = "Intervalo",
    y = "hi (%)"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(sec_B$tabla, aes(x = Intervalo, y = Porcentaje_Esperado)) +
  geom_col(fill = "gray35", color = "black") +
  labs(
    title = "Gráfica N°5: Frecuencias esperadas — Modelo Geométrico Sección B",
    x = "Intervalo",
    y = "Porcentaje esperado (%)"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Resumen de pruebas de bondad de ajuste

resumen <- data.frame(
  Seccion = c(sec_A$nombre, sec_B$nombre),
  Modelo = c("Geométrico", "Geométrico"),
  Parametro_p = c(round(sec_A$p_hat, 4), round(sec_B$p_hat, 4)),
  Pearson = paste0(c(round(sec_A$pearson * 100, 2), round(sec_B$pearson * 100, 2)), "%"),
  Chi_Cuadrado = c(round(sec_A$chi_calculado, 4), round(sec_B$chi_calculado, 4)),
  GL = c(sec_A$gl, sec_B$gl),
  Valor_p = c(round(sec_A$p_valor, 4), round(sec_B$p_valor, 4)),
  Valor_Critico = c(round(sec_A$chi_critico, 4), round(sec_B$chi_critico, 4)),
  Decision = c(sec_A$decision, sec_B$decision),
  Observacion = c(sec_A$nota_gl, sec_B$nota_gl)
)

kable(resumen, caption = "Tabla N°4: Resumen de Pearson y Chi-cuadrado por sección") %>%
  kable_styling(full_width = FALSE, position = "center")
Tabla N°4: Resumen de Pearson y Chi-cuadrado por sección
Seccion Modelo Parametro_p Pearson Chi_Cuadrado GL Valor_p Valor_Critico Decision Observacion
Sección A Geométrico 0.50 97.53% 2.8020 2 0.2464 5.9915 No se rechaza H0: modelo aceptado Chi-cuadrado formal aplicable.
Sección B Geométrico 0.75 100% 0.3333 1 0.5637 3.8415 No se rechaza H0: modelo aceptado Chi-cuadrado referencial: la sección tiene pocas clases para una prueba formal con parámetro estimado.

8 Conclusiones

Se construyó una sola tabla de distribución de frecuencias para YEARS_ACTIVE aplicando la regla de Sturges y respetando un máximo de 10 intervalos. A partir del diagrama de barras se observó que la variable completa no presenta un único patrón geométrico uniforme. Por ello, el primer intervalo [1–11] se omitió del ajuste por representar pozos recientes con alta concentración, mientras que el análisis se dividió en dos secciones: una sección principal de [12–22] a [56–66] y una sección extrema final de [67–77] a [78–82]. En cada sección se estimó el parámetro del modelo geométrico, se calcularon las frecuencias esperadas y se aplicaron el test de Pearson y la prueba Chi-cuadrado. La sección que presente mejor aceptación según Pearson y Chi-cuadrado se considera la más representativa para el modelo geométrico.