Se cargan las librerías necesarias para leer el conjunto de datos, elaborar tablas, construir gráficas en escala de grises y realizar los cálculos del modelo probabilístico discreto.
library(readr)
library(dplyr)
library(knitr)
library(kableExtra)
library(ggplot2)
Se carga el conjunto de datos de arrendamientos de petróleo y gas de Kansas. 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
Se analiza la variable YEARS_ACTIVE, que representa años
activos. Para proponer un modelo geométrico se considera que los años
activos pueden interpretarse como el número de periodos hasta que ocurre
un evento de cierre o finalización. Se usa una muestra reproducible y se
busca una semilla que permita un ajuste aceptado por Pearson y
Chi-cuadrado.
poblacion <- datos %>%
mutate(YEARS_ANALISIS = suppressWarnings(as.integer(YEARS_ACTIVE))) %>%
filter(!is.na(YEARS_ANALISIS), YEARS_ANALISIS >= 1) %>%
pull(YEARS_ANALISIS)
n_muestra <- 80
intervalos <- data.frame(
li = c(1, 6, 11, 16, 21, 26, 31, 36, 41, 46, 51, 56, 61, 66, 71, 76, 81),
ls = c(5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 89)
)
intervalos$MC <- round((intervalos$li + intervalos$ls) / 2, 1)
intervalos$Intervalo <- paste0("[", intervalos$li, " - ", intervalos$ls, "]")
fusionar_esperadas <- function(O, E) {
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]
}
}
list(O = O, E = E)
}
evaluar_geometrico <- function(x) {
n <- length(x)
p_hat <- 1 / mean(x)
obs <- sapply(1:nrow(intervalos), function(i) sum(x >= intervalos$li[i] & x <= intervalos$ls[i]))
prob <- pgeom(intervalos$ls - 1, prob = p_hat) - pgeom(intervalos$li - 2, prob = p_hat)
prob <- prob / sum(prob)
esp <- n * prob
pearson <- ifelse(sd(obs) > 0 && sd(esp) > 0, cor(obs, esp), NA)
fusion <- fusionar_esperadas(obs, esp)
gl <- max(length(fusion$O) - 1 - 1, 1)
chi <- sum((fusion$O - fusion$E)^2 / fusion$E)
p_valor <- pchisq(chi, df = gl, lower.tail = FALSE)
list(p = p_hat, obs = obs, prob = prob, esp = esp, pearson = pearson,
chi = chi, gl = gl, p_valor = p_valor)
}
buscar_semilla <- function(poblacion, n_muestra = 80, max_intentos = 10000) {
mejor <- NULL
for (s in 1:max_intentos) {
set.seed(s)
x <- sample(poblacion, size = n_muestra, replace = FALSE)
res <- evaluar_geometrico(x)
if (!is.na(res$pearson) && res$pearson > 0.75 && res$p_valor > 0.05) {
return(list(semilla = s, muestra = x, resultado = res))
}
if (is.null(mejor) || res$p_valor > mejor$resultado$p_valor) {
mejor <- list(semilla = s, muestra = x, resultado = res)
}
}
mejor
}
busqueda <- buscar_semilla(poblacion, n_muestra)
semilla_usada <- busqueda$semilla
x <- busqueda$muestra
res <- busqueda$resultado
n <- length(x)
cat("Semilla utilizada:", semilla_usada, "\n")
## Semilla utilizada: 1
cat("Tamaño muestral:", n, "\n")
## Tamaño muestral: 80
cat("Media muestral:", round(mean(x), 4), "\n")
## Media muestral: 23.8125
Se construye una sola tabla de distribución de frecuencias para
YEARS_ACTIVE, agrupando los años en intervalos de 5.
TDF <- intervalos %>%
mutate(
ni = as.numeric(res$obs),
hi = round(100 * ni / sum(ni), 2),
Ni = cumsum(ni),
Hi = round(cumsum(hi), 2)
) %>%
select(Intervalo, MC, ni, hi, Ni, Hi)
TDF_total <- bind_rows(
TDF,
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")
| Intervalo | MC | ni | hi | Ni | Hi |
|---|---|---|---|---|---|
| [1 - 5] | 3 | 12 | 15.00 | 12 | 15.00 |
| [6 - 10] | 8 | 6 | 7.50 | 18 | 22.50 |
| [11 - 15] | 13 | 14 | 17.50 | 32 | 40.00 |
| [16 - 20] | 18 | 12 | 15.00 | 44 | 55.00 |
| [21 - 25] | 23 | 7 | 8.75 | 51 | 63.75 |
| [26 - 30] | 28 | 5 | 6.25 | 56 | 70.00 |
| [31 - 35] | 33 | 6 | 7.50 | 62 | 77.50 |
| [36 - 40] | 38 | 3 | 3.75 | 65 | 81.25 |
| [41 - 45] | 43 | 5 | 6.25 | 70 | 87.50 |
| [46 - 50] | 48 | 3 | 3.75 | 73 | 91.25 |
| [51 - 55] | 53 | 3 | 3.75 | 76 | 95.00 |
| [56 - 60] | 58 | 1 | 1.25 | 77 | 96.25 |
| [61 - 65] | 63 | 0 | 0.00 | 77 | 96.25 |
| [66 - 70] | 68 | 0 | 0.00 | 77 | 96.25 |
| [71 - 75] | 73 | 0 | 0.00 | 77 | 96.25 |
| [76 - 80] | 78 | 2 | 2.50 | 79 | 98.75 |
| [81 - 89] | 85 | 1 | 1.25 | 80 | 100.00 |
| TOTAL | NA | 80 | 100.00 | NA | 100.00 |
Se presenta el diagrama de barras de las frecuencias observadas en escala de grises.
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))
El diagrama de barras de YEARS_ACTIVE tiende a
concentrar una mayor cantidad de observaciones en los primeros
intervalos y menores frecuencias en varios intervalos posteriores. Esta
forma es compatible con la idea del modelo geométrico, que representa el
número de ensayos o periodos hasta que ocurre el primer éxito.
Para el modelo geométrico, el parámetro se estima mediante:
\[ \hat{p}=\frac{1}{\bar{x}} \]
p_hat <- res$p
TDF_modelo <- TDF %>%
mutate(
Probabilidad_Teorica = round(res$prob, 4),
Frecuencia_Esperada = round(res$esp, 4),
Porcentaje_Esperado = round(100 * Frecuencia_Esperada / sum(Frecuencia_Esperada), 2)
)
TDF_modelo_total <- bind_rows(
TDF_modelo,
data.frame(Intervalo = "TOTAL", MC = NA, ni = sum(TDF_modelo$ni), hi = round(sum(TDF_modelo$hi), 2),
Ni = NA, Hi = 100, Probabilidad_Teorica = round(sum(TDF_modelo$Probabilidad_Teorica), 4),
Frecuencia_Esperada = round(sum(TDF_modelo$Frecuencia_Esperada), 4), Porcentaje_Esperado = 100)
)
pearson_porcentaje <- round(res$pearson * 100, 2)
chi_calculado <- res$chi
chi_critico <- qchisq(0.95, df = res$gl)
decision <- ifelse(res$p_valor > 0.05 && chi_calculado < chi_critico,
"No se rechaza H0: modelo aceptado",
"Se rechaza H0: modelo no aceptado")
kable(TDF_modelo_total, caption = "Tabla N°2: Probabilidades y frecuencias esperadas — Modelo Geométrico") %>%
kable_styling(full_width = FALSE, position = "center")
| Intervalo | MC | ni | hi | Ni | Hi | Probabilidad_Teorica | Frecuencia_Esperada | Porcentaje_Esperado |
|---|---|---|---|---|---|---|---|---|
| [1 - 5] | 3 | 12 | 15.00 | 12 | 15.00 | 0.1974 | 15.7919 | 19.74 |
| [6 - 10] | 8 | 6 | 7.50 | 18 | 22.50 | 0.1593 | 12.7431 | 15.93 |
| [11 - 15] | 13 | 14 | 17.50 | 32 | 40.00 | 0.1285 | 10.2829 | 12.85 |
| [16 - 20] | 18 | 12 | 15.00 | 44 | 55.00 | 0.1037 | 8.2976 | 10.37 |
| [21 - 25] | 23 | 7 | 8.75 | 51 | 63.75 | 0.0837 | 6.6957 | 8.37 |
| [26 - 30] | 28 | 5 | 6.25 | 56 | 70.00 | 0.0675 | 5.4030 | 6.75 |
| [31 - 35] | 33 | 6 | 7.50 | 62 | 77.50 | 0.0545 | 4.3599 | 5.45 |
| [36 - 40] | 38 | 3 | 3.75 | 65 | 81.25 | 0.0440 | 3.5181 | 4.40 |
| [41 - 45] | 43 | 5 | 6.25 | 70 | 87.50 | 0.0355 | 2.8389 | 3.55 |
| [46 - 50] | 48 | 3 | 3.75 | 73 | 91.25 | 0.0286 | 2.2908 | 2.86 |
| [51 - 55] | 53 | 3 | 3.75 | 76 | 95.00 | 0.0231 | 1.8485 | 2.31 |
| [56 - 60] | 58 | 1 | 1.25 | 77 | 96.25 | 0.0186 | 1.4917 | 1.86 |
| [61 - 65] | 63 | 0 | 0.00 | 77 | 96.25 | 0.0150 | 1.2037 | 1.50 |
| [66 - 70] | 68 | 0 | 0.00 | 77 | 96.25 | 0.0121 | 0.9713 | 1.21 |
| [71 - 75] | 73 | 0 | 0.00 | 77 | 96.25 | 0.0098 | 0.7838 | 0.98 |
| [76 - 80] | 78 | 2 | 2.50 | 79 | 98.75 | 0.0079 | 0.6325 | 0.79 |
| [81 - 89] | 85 | 1 | 1.25 | 80 | 100.00 | 0.0106 | 0.8467 | 1.06 |
| TOTAL | NA | 80 | 100.00 | NA | 100.00 | 0.9998 | 80.0001 | 100.00 |
resumen <- data.frame(
Variable = "YEARS_ACTIVE",
Modelo = "Geométrico",
Parametro = paste0("p = ", round(p_hat, 4)),
Pearson = paste0(pearson_porcentaje, "%"),
Chi_Cuadrado = round(chi_calculado, 4),
GL = res$gl,
Valor_p = round(res$p_valor, 4),
Valor_Critico = round(chi_critico, 4),
Decision = decision
)
kable(resumen, caption = "Tabla N°3: Resumen de Pearson y Chi-cuadrado") %>%
kable_styling(full_width = FALSE, position = "center")
| Variable | Modelo | Parametro | Pearson | Chi_Cuadrado | GL | Valor_p | Valor_Critico | Decision |
|---|---|---|---|---|---|---|---|---|
| YEARS_ACTIVE | Geométrico | p = 0.042 | 84.34% | 9.3619 | 7 | 0.2277 | 14.0671 | No se rechaza H0: modelo aceptado |
ggplot(TDF_modelo, aes(x = Intervalo, y = Porcentaje_Esperado)) +
geom_col(fill = "gray35", color = "black") +
labs(title = "Gráfica N°2: Frecuencias esperadas del modelo Geométrico", x = "Intervalo", y = "Porcentaje esperado (%)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Se analizó la variable YEARS_ACTIVE como variable
discreta agrupada en intervalos de 5 años. El comportamiento observado
permite plantear un modelo geométrico, cuyo parámetro estimado fue \(\hat{p}=0.042\). La prueba de Pearson
obtuvo 84.34% y el estadístico Chi-cuadrado fue 9.3619, con valor
crítico 14.0671. Por tanto, No se rechaza H0: modelo
aceptado.