Se retoma la variable Década de Finalización de perforación (Completion Year/Month/Day) de los pozos de petróleo y gas de Nueva York, ahora desde la estadística inferencial: se propone y valida un modelo de probabilidad discreto mediante la prueba de bondad de ajuste de Pearson (Chi-cuadrado).
Variable redefinida:
X = N° de pozos cuya perforación fue completada por semana calendario, periodo 2020-01-01 – 2025-12-31.
Se usa escala semanal (Poisson requiere λ constante) y se acota a 2020-2025 por ser el tramo más reciente y estable, evitando heterogeneidad de periodos anteriores (auges/caídas históricas del sector).
# ==========================================================================
# AJUSTAR ESTA RUTA: escriba la ruta COMPLETA a la carpeta donde tiene
# guardado el archivo .csv en su computador. Ejemplo Windows:
# ruta_archivo <- "C:/Users/ASUS/Desktop/Estadistica/new_york_exel/Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv"
# Ejemplo si el .Rmd y el .csv están en la MISMA carpeta, no hace falta ruta:
# ruta_archivo <- "Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv"
# ==========================================================================
ruta_archivo <- "Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv"
if (!file.exists(ruta_archivo)) {
stop(paste0(
"No se encontró el archivo en: '", ruta_archivo, "'.\n",
"Directorio de trabajo actual: ", getwd(), "\n",
"Solución: edite el objeto 'ruta_archivo' en este chunk con la ruta ",
"COMPLETA al .csv en su computador, o mueva el .csv a la carpeta: ", getwd()
))
}
Datos_Brutos <- read.csv(
ruta_archivo,
header = TRUE, sep = ";", fileEncoding = "latin1"
)
# Detección automática de separador: si el archivo no se separó en columnas
# (algunas versiones descargadas del portal usan "," en vez de ";"), se
# vuelve a leer con el separador correcto.
if (ncol(Datos_Brutos) <= 1) {
Datos_Brutos <- read.csv(
ruta_archivo,
header = TRUE, sep = ",", fileEncoding = "latin1"
)
}
if (ncol(Datos_Brutos) <= 1) {
Datos_Brutos <- read.csv(ruta_archivo, header = TRUE, sep = ";")
}
if (ncol(Datos_Brutos) <= 1) {
Datos_Brutos <- read.csv(ruta_archivo, header = TRUE, sep = ",")
}
if (ncol(Datos_Brutos) <= 1) {
stop("No se pudo separar el archivo en columnas. Verifique el delimitador manualmente.")
}
cat("Archivo leído desde:", normalizePath(ruta_archivo), "\n")## Archivo leído desde: C:\Users\ASUS\Downloads\Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv
cat("Columnas detectadas:", ncol(Datos_Brutos), "| Registros totales cargados:", nrow(Datos_Brutos), "\n")## Columnas detectadas: 52 | Registros totales cargados: 47390
El dataset puede traer la fecha de finalización de dos formas según
la versión de descarga: (A) separada en tres columnas
Completion Year / Completion Month /
Completion Day, o (B) en un solo campo de fecha tipo
Date Well Completed. Se detectan las columnas por
patrón (no por nombre fijo) para que el script funcione sin
importar pequeñas diferencias de nombre (espacios extra, mayúsculas,
etc.).
nombres_disp <- names(Datos_Brutos)
col_y <- nombres_disp[grepl("completion", nombres_disp, ignore.case = TRUE) &
grepl("year", nombres_disp, ignore.case = TRUE)]
col_m <- nombres_disp[grepl("completion", nombres_disp, ignore.case = TRUE) &
grepl("month", nombres_disp, ignore.case = TRUE)]
col_d <- nombres_disp[grepl("completion", nombres_disp, ignore.case = TRUE) &
grepl("day", nombres_disp, ignore.case = TRUE) &
!grepl("decade", nombres_disp, ignore.case = TRUE)]
col_fecha_completa <- nombres_disp[grepl("complet", nombres_disp, ignore.case = TRUE) &
grepl("date", nombres_disp, ignore.case = TRUE)]
if (length(col_y) > 0 && length(col_m) > 0 && length(col_d) > 0) {
modo_fecha_fin <- "ymd"
cat("Modo detectado: columnas separadas Year/Month/Day\n")
cat(" Year :", col_y[1], "\n Month:", col_m[1], "\n Day :", col_d[1], "\n")
} else if (length(col_fecha_completa) > 0) {
modo_fecha_fin <- "fecha_texto"
cat("Modo detectado: campo de fecha único\n")
cat(" Columna:", col_fecha_completa[1], "\n")
} else {
stop(paste0(
"No se encontró ninguna columna de finalización (se buscaron patrones ",
"'Completion'+'Year'/'Month'/'Day' y 'Complet*'+'Date').\n",
"COLUMNAS ENCONTRADAS EN EL ARCHIVO:\n", paste(nombres_disp, collapse = " | ")
))
}## Modo detectado: campo de fecha único
## Columna: Date.Well.Completed
if (modo_fecha_fin == "ymd") {
Datos_Brutos <- Datos_Brutos %>%
mutate(
Anio_c = suppressWarnings(as.integer(.data[[col_y[1]]])),
Mes_c = suppressWarnings(as.integer(.data[[col_m[1]]])),
Dia_c = suppressWarnings(as.integer(.data[[col_d[1]]])),
Fecha_Completado = suppressWarnings(as.Date(
paste(Anio_c, Mes_c, Dia_c, sep = "-"), format = "%Y-%m-%d"
))
)
} else {
Datos_Brutos <- Datos_Brutos %>%
mutate(
Fecha_Completado = suppressWarnings(as.Date(
.data[[col_fecha_completa[1]]], format = "%m/%d/%Y"
))
)
}
# --- Filtro al periodo 2020-01-01 a 2025-12-31 (ver justificación arriba) ---
Datos_Validos <- Datos_Brutos %>%
filter(
!is.na(Fecha_Completado),
Fecha_Completado >= as.Date("2020-01-01"),
Fecha_Completado <= as.Date("2025-12-31")
)
# Conteo de pozos por día (paso intermedio, para no perder ningún día sin actividad)
conteo_diario <- Datos_Validos %>%
count(Fecha_Completado, name = "pozos_dia")
rango_fechas <- seq(as.Date("2020-01-01"), as.Date("2025-12-31"), by = "day")
Serie_Diaria <- data.frame(Fecha_Completado = rango_fechas) %>%
left_join(conteo_diario, by = "Fecha_Completado") %>%
mutate(pozos_dia = ifelse(is.na(pozos_dia), 0, pozos_dia))
# --- Agregación a nivel SEMANAL: esta es la variable X del análisis ---
Serie_Semanal <- Serie_Diaria %>%
mutate(Semana = floor_date(Fecha_Completado, unit = "week")) %>%
group_by(Semana) %>%
summarise(pozos_semana = sum(pozos_dia), .groups = "drop")
# Completar TODAS las semanas del rango, incluidas las que tuvieron 0 pozos
rango_semanas <- seq(min(Serie_Semanal$Semana), max(Serie_Semanal$Semana), by = "week")
Serie_Semanal <- data.frame(Semana = rango_semanas) %>%
left_join(Serie_Semanal, by = "Semana") %>%
mutate(pozos_semana = ifelse(is.na(pozos_semana), 0, pozos_semana))
# X y n representan la escala SEMANAL en todo el resto del documento
X <- Serie_Semanal$pozos_semana
n <- length(X)
if (n == 0) stop("ERROR: No hay datos válidos.")
lambda_hat <- mean(X)
var_X <- var(X)
cat("Variable analizada: N° de pozos con finalización de perforación por SEMANA\n")## Variable analizada: N° de pozos con finalización de perforación por SEMANA
## Periodo: 2019-12-29 a 2025-12-28
## Número de semanas observadas (n): 314
## Media (lambda estimado, x̄): 1.2006
## Varianza: 1.5539
## Razón varianza/media: 1.294
## Asimetría: 1.46 | Curtosis: 6.4458
tabla_FO <- as.data.frame(table(X))
names(tabla_FO) <- c("x", "FOi")
tabla_FO$x <- as.numeric(as.character(tabla_FO$x))
tabla_FO <- tabla_FO %>%
arrange(x) %>%
mutate(
fi = FOi / n,
Fi_asc = cumsum(fi)
)
tabla_FO %>%
gt() %>%
tab_header(
title = md("**DISTRIBUCIÓN DE FRECUENCIAS SEMANALES**"),
subtitle = md("Variable: **N° de pozos con finalización de perforación por semana (X)** · Nueva York · 2020-2025")
) %>%
fmt_number(columns = c(fi, Fi_asc), decimals = 4) %>%
cols_label(
x = "N° de pozos por semana (x)", FOi = "Frec. Observada (FOi)",
fi = "Frec. Relativa (fi)", Fi_asc = "Frec. Relativa Acum. (Fi)"
) %>%
cols_align(align = "center", columns = everything()) %>%
tab_style(
style = list(cell_fill(color = col_principal), cell_text(color = "white", weight = "bold")),
locations = cells_title()
) %>%
tab_style(
style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()
) %>%
opt_row_striping() %>%
opt_table_font(font = google_font("Roboto")) %>%
tab_options(
table.font.size = px(13),
heading.align = "left",
data_row.padding = px(6),
table.border.top.color = col_principal,
table.border.bottom.color = col_principal,
column_labels.border.bottom.color = col_principal
) %>%
tab_source_note(md("*Fuente: NYS DEC — Oil, Gas & Other Regulated Wells. Elaboración: EDUARDO.*"))| DISTRIBUCIÓN DE FRECUENCIAS SEMANALES | |||
| Variable: N° de pozos con finalización de perforación por semana (X) · Nueva York · 2020-2025 | |||
| N° de pozos por semana (x) | Frec. Observada (FOi) | Frec. Relativa (fi) | Frec. Relativa Acum. (Fi) |
|---|---|---|---|
| 0 | 107 | 0.3408 | 0.3408 |
| 1 | 104 | 0.3312 | 0.6720 |
| 2 | 61 | 0.1943 | 0.8662 |
| 3 | 28 | 0.0892 | 0.9554 |
| 4 | 8 | 0.0255 | 0.9809 |
| 5 | 3 | 0.0096 | 0.9904 |
| 6 | 2 | 0.0064 | 0.9968 |
| 8 | 1 | 0.0032 | 1.0000 |
| Fuente: NYS DEC — Oil, Gas & Other Regulated Wells. Elaboración: EDUARDO. | |||
Al tratarse de una variable discreta de conteo (no de intervalos de clase), se usa un gráfico de barras —el análogo discreto del histograma— para observar la forma de la distribución semanal.
ggplot(tabla_FO, aes(x = factor(x), y = fi * 100)) +
geom_col(fill = col_barras, width = 0.7) +
labs(
title = "Gráfico N°1: Frecuencia relativa observada — pozos finalizados por semana",
x = "N° de pozos finalizados por semana (x)", y = "Frecuencia relativa (%)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(color = col_principal, face = "bold", size = 12),
axis.title = element_text(color = col_principal),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 7)
)Se evalúan las cuatro distribuciones discretas candidatas frente a la naturaleza de X:
comparacion <- data.frame(
Distribución = c("Bernoulli", "Binomial", "Geométrica", "Poisson"),
`Qué mide` = c(
"Éxito/fracaso en un único ensayo",
"N° de éxitos en n ensayos fijos con prob. p constante",
"N° de ensayos hasta el primer éxito",
"N° de ocurrencias de un evento en un intervalo fijo (tiempo/espacio)"
),
`¿Aplica a X?` = c(
"No: X no es binaria (0/1)",
"No: no existe un n° fijo de \"ensayos\" por semana",
"No: X no mide espera hasta un éxito",
"Sí: X cuenta eventos (finalizaciones) por unidad de tiempo fija (1 semana)"
),
check.names = FALSE
)
comparacion %>%
gt() %>%
tab_header(title = md("**SELECCIÓN DEL MODELO DE PROBABILIDAD**")) %>%
cols_align(align = "left", columns = everything()) %>%
tab_style(
style = list(cell_fill(color = col_principal), cell_text(color = "white", weight = "bold")),
locations = cells_title()
) %>%
tab_style(
style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()
) %>%
tab_style(
style = list(cell_fill(color = "#D0ECE7"), cell_text(weight = "bold")),
locations = cells_body(rows = Distribución == "Poisson")
) %>%
opt_table_font(font = google_font("Roboto")) %>%
tab_options(table.font.size = px(12.5), heading.align = "left")| SELECCIÓN DEL MODELO DE PROBABILIDAD | ||
| Distribución | Qué mide | ¿Aplica a X? |
|---|---|---|
| Bernoulli | Éxito/fracaso en un único ensayo | No: X no es binaria (0/1) |
| Binomial | N° de éxitos en n ensayos fijos con prob. p constante | No: no existe un n° fijo de "ensayos" por semana |
| Geométrica | N° de ensayos hasta el primer éxito | No: X no mide espera hasta un éxito |
| Poisson | N° de ocurrencias de un evento en un intervalo fijo (tiempo/espacio) | Sí: X cuenta eventos (finalizaciones) por unidad de tiempo fija (1 semana) |
Conclusión: X cuenta eventos discretos e independientes en un intervalo fijo de tiempo a tasa aproximadamente constante — la definición de un proceso de Poisson. Se modela X ~ Poisson(λ̂), con λ̂ = x̄ = 1.2006.
Para la prueba Chi-cuadrado se agrupan las clases con FEi < 5 (regla de Cochran), tanto en la cola derecha como en la izquierda si corresponde.
construir_clases_poisson <- function(x, lambda) {
n_obs <- length(x)
# --- Cola izquierda: agrupa x = 0..(k_ini-1) si su FEi < 5 ---
k_ini <- 0
while (n_obs * dpois(k_ini, lambda) < 5 && k_ini <= 200) k_ini <- k_ini + 1
FOi <- c(); Pi <- c(); etiquetas <- c()
if (k_ini > 0) {
FOi <- c(FOi, sum(x < k_ini))
Pi <- c(Pi, ppois(k_ini - 1, lambda))
etiquetas <- c(etiquetas, paste0(k_ini - 1, " o menos"))
}
# --- Clases individuales mientras FEi >= 5 ---
k <- k_ini
repeat {
p_k <- dpois(k, lambda)
FE_k <- n_obs * p_k
if (FE_k < 5) break
FOi <- c(FOi, sum(x == k))
Pi <- c(Pi, p_k)
etiquetas <- c(etiquetas, as.character(k))
k <- k + 1
if (k > 200) break
}
# --- Cola derecha: agrupa x >= k ---
FOi <- c(FOi, sum(x >= k))
Pi <- c(Pi, 1 - ppois(k - 1, lambda))
etiquetas <- c(etiquetas, paste0(k, " o más"))
data.frame(clase = etiquetas, FOi = FOi, Pi = Pi, FEi = n_obs * Pi)
}
Tabla_Chi <- construir_clases_poisson(X, lambda_hat)
Tabla_Chi## clase FOi Pi FEi
## 1 0 107 0.301002430 94.514763
## 2 1 104 0.361394637 113.477916
## 3 2 61 0.216951876 68.122889
## 4 3 28 0.086826812 27.263619
## 5 4 8 0.026061870 8.183427
## 6 5 o más 6 0.007762376 2.437386
plot_df <- Tabla_Chi %>%
select(clase, FOi, FEi) %>%
pivot_longer(cols = c(FOi, FEi), names_to = "Tipo", values_to = "Frecuencia") %>%
mutate(clase = factor(clase, levels = Tabla_Chi$clase))
ggplot(plot_df, aes(x = clase, y = Frecuencia, fill = Tipo)) +
geom_col(position = position_dodge(width = 0.75), width = 0.65) +
scale_fill_manual(
values = c(FOi = col_barras, FEi = col_acento),
labels = c(FOi = "Observada (FOi)", FEi = "Esperada · Poisson (FEi)")
) +
labs(
title = "Gráfico N°2: Frecuencias observadas y esperadas — Modelo Poisson",
subtitle = paste0("λ̂ = ", round(lambda_hat, 4)),
x = "N° de pozos por semana", y = "Frecuencia (semanas)", fill = ""
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(color = col_principal, face = "bold", size = 12),
legend.position = "top",
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 7)
)El siguiente gráfico es el soporte visual central de la validación: compara la función de probabilidad acumulada empírica de los datos con la CDF teórica de la Poisson(λ̂). Cuanto más se superpongan ambas curvas escalonadas, mejor es el ajuste del modelo.
x_max_plot <- max(X)
x_seq <- 0:x_max_plot
cdf_empirica_fn <- ecdf(X)
cdf_emp_vals <- cdf_empirica_fn(x_seq)
cdf_teo_vals <- ppois(x_seq, lambda_hat)
df_cdf <- data.frame(
x = rep(x_seq, 2),
F = c(cdf_emp_vals, cdf_teo_vals),
Tipo = rep(c("Empírica (datos)", "Teórica (Poisson)"), each = length(x_seq))
)
ggplot(df_cdf, aes(x = x, y = F, color = Tipo)) +
geom_step(linewidth = 1) +
geom_point(size = 1.6) +
scale_color_manual(values = c(
"Empírica (datos)" = col_barras,
"Teórica (Poisson)" = col_acento
)) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Gráfico N°3: Función de Probabilidad Acumulada (CDF) — Empírica y Poisson",
subtitle = paste0("X ~ Poisson(λ̂ = ", round(lambda_hat, 4), ")"),
x = "N° de pozos finalizados por semana (x)", y = "P(X ≤ x)", color = ""
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(color = col_principal, face = "bold", size = 12),
legend.position = "top"
)Se contrasta:
\[X^2 = \sum \frac{(FO_i - FE_i)^2}{FE_i} \qquad FE_i = P_i \times n \qquad gl = k - 1 - m\]
donde k es el número de clases y m el número de parámetros estimados a partir de los datos (en este caso, m = 1, ya que λ se estimó como x̄).
Tabla_Chi <- Tabla_Chi %>%
mutate(Aporte_Chi2 = (FOi - FEi)^2 / FEi)
chi2_calculado <- sum(Tabla_Chi$Aporte_Chi2)
k_clases <- nrow(Tabla_Chi)
m_parametros <- 1
gl <- k_clases - 1 - m_parametros
alpha <- 0.05
chi2_critico <- qchisq(1 - alpha, df = gl)
decision <- if (chi2_calculado <= chi2_critico) {
"No se rechaza H0 (el modelo Poisson es un ajuste adecuado)"
} else {
"Se rechaza H0 (el modelo Poisson no se ajusta adecuadamente)"
}
cat("Chi-cuadrado calculado (X² calc):", round(chi2_calculado, 4), "\n")## Chi-cuadrado calculado (X² calc): 8.417
## Número de clases (k): 6 | Parámetros estimados (m): 1
## Grados de libertad (gl = k - 1 - m): 4
## Chi-cuadrado crítico (alpha = 0.05): 9.4877
## Decisión: No se rechaza H0 (el modelo Poisson es un ajuste adecuado)
fila_max_aporte <- which.max(Tabla_Chi$Aporte_Chi2)
Tabla_Chi %>%
gt() %>%
tab_header(
title = md("**VALIDACIÓN DEL MODELO POISSON**"),
subtitle = md(paste0("λ̂ = ", round(lambda_hat, 4), " · n = ", n, " semanas · 2020-2025"))
) %>%
fmt_number(columns = c(Pi, FEi, Aporte_Chi2), decimals = 4) %>%
cols_label(
clase = "Clase (x)", FOi = "FOi", Pi = "Pi (teórica)",
FEi = "FEi", Aporte_Chi2 = "Aporte a X²"
) %>%
cols_align(align = "center", columns = everything()) %>%
tab_style(
style = list(cell_fill(color = col_principal), cell_text(color = "white", weight = "bold")),
locations = cells_title()
) %>%
tab_style(
style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()
) %>%
tab_style(
style = list(cell_fill(color = "#FDEBD0"), cell_text(weight = "bold")),
locations = cells_body(rows = fila_max_aporte)
) %>%
opt_row_striping() %>%
opt_table_font(font = google_font("Roboto")) %>%
tab_options(
table.font.size = px(13),
heading.align = "left",
data_row.padding = px(7),
table.border.top.color = col_principal,
table.border.bottom.color = col_principal,
column_labels.border.bottom.color = col_principal
) %>%
tab_source_note(md(paste0(
"**X² calculado = ", round(chi2_calculado, 4),
"** y **X² crítico (α = 0.05, gl = ", gl, ") = ", round(chi2_critico, 4), "** | ",
"**Decisión:** ", decision
)))| VALIDACIÓN DEL MODELO POISSON | ||||
| λ̂ = 1.2006 · n = 314 semanas · 2020-2025 | ||||
| Clase (x) | FOi | Pi (teórica) | FEi | Aporte a X² |
|---|---|---|---|---|
| 0 | 107 | 0.3010 | 94.5148 | 1.6493 |
| 1 | 104 | 0.3614 | 113.4779 | 0.7916 |
| 2 | 61 | 0.2170 | 68.1229 | 0.7448 |
| 3 | 28 | 0.0868 | 27.2636 | 0.0199 |
| 4 | 8 | 0.0261 | 8.1834 | 0.0041 |
| 5 o más | 6 | 0.0078 | 2.4374 | 5.2073 |
| X² calculado = 8.417 y X² crítico (α = 0.05, gl = 4) = 9.4877 | Decisión: No se rechaza H0 (el modelo Poisson es un ajuste adecuado) | ||||