Estadística Inferencial — Modelo de Probabilidad Discreto aplicado al Inicio de Perforación (NY)
Eduardo Loayza
2026-07-16
Introducción y Metodología
Se retoma la variable Fecha de Inicio de Perforación (Date Spudded) 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 se inició 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.
1. Carga de Librerías
2. Carga de Datos
# ==========================================================================
# 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.csvcat("Columnas detectadas:", ncol(Datos_Brutos), "| Registros totales cargados:", nrow(Datos_Brutos), "\n")## Columnas detectadas: 52 | Registros totales cargados: 473903. Extracción y Preparación de la Variable
Datos_Brutos <- Datos_Brutos %>%
mutate(Fecha_Spud = as.Date(Date.Spudded, 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_Spud),
Fecha_Spud >= as.Date("2020-01-01"),
Fecha_Spud <= 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_Spud, name = "pozos_dia")
rango_fechas <- seq(as.Date("2020-01-01"), as.Date("2025-12-31"), by = "day")
Serie_Diaria <- data.frame(Fecha_Spud = rango_fechas) %>%
left_join(conteo_diario, by = "Fecha_Spud") %>%
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_Spud, 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 inicio de perforación por SEMANA\n")## Variable analizada: N° de pozos con inicio 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.0796## Varianza: 1.3515## Razón varianza/media: 1.252## Asimetría: 1.5923 | Curtosis: 8.93534. Tabla de Distribución de Frecuencias
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 inicio 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 inicio 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 | 121 | 0.3854 | 0.3854 |
| 1 | 97 | 0.3089 | 0.6943 |
| 2 | 60 | 0.1911 | 0.8854 |
| 3 | 27 | 0.0860 | 0.9713 |
| 4 | 8 | 0.0255 | 0.9968 |
| 9 | 1 | 0.0032 | 1.0000 |
| Fuente: NYS DEC — Oil, Gas & Other Regulated Wells. Elaboración: EDUARDO. | |||
5. Gráfico de Barras — Frecuencia Relativa Observada
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 iniciados por semana",
x = "N° de pozos iniciados 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)
)
6. Justificación del Modelo de Probabilidad
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 (spuds) 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 (spuds) 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.0796.
6.1 Frecuencias Observadas vs. Esperadas (Modelo Poisson)
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 121 0.339725332 106.673754
## 2 1 97 0.366773527 115.166888
## 3 2 60 0.197987621 62.168113
## 4 3 27 0.071250322 22.372601
## 5 4 8 0.019230780 6.038465
## 6 5 o más 1 0.005032418 1.580179plot_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 vs. 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)
)
6.2 Función de Probabilidad Acumulada (CDF): Empírica vs. Teórica
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 vs. Poisson",
subtitle = paste0("X ~ Poisson(λ̂ = ", round(lambda_hat, 4), ")"),
x = "N° de pozos iniciados 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"
)
7. Test de Pearson y Chi-Cuadrado
Se contrasta:
- H0: X se distribuye según una Poisson(λ̂) — el modelo propuesto es adecuado.
- Ha: X no se distribuye según una Poisson(λ̂) — el modelo propuesto no es adecuado.
\[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): 6.6726## 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)8. Tabla Resumen de Validación del Modelo
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),
"** vs. **X² crítico (α = 0.05, gl = ", gl, ") = ", round(chi2_critico, 4), "** | ",
"**Decisión:** ", decision
)))| VALIDACIÓN DEL MODELO POISSON | ||||
| λ̂ = 1.0796 · n = 314 semanas · 2020-2025 | ||||
| Clase (x) | FOi | Pi (teórica) | FEi | Aporte a X² |
|---|---|---|---|---|
| 0 | 121 | 0.3397 | 106.6738 | 1.9240 |
| 1 | 97 | 0.3668 | 115.1669 | 2.8657 |
| 2 | 60 | 0.1980 | 62.1681 | 0.0756 |
| 3 | 27 | 0.0713 | 22.3726 | 0.9571 |
| 4 | 8 | 0.0192 | 6.0385 | 0.6372 |
| 5 o más | 1 | 0.0050 | 1.5802 | 0.2130 |
| X² calculado = 6.6726 vs. X² crítico (α = 0.05, gl = 4) = 9.4877 | Decisión: No se rechaza H0 (el modelo Poisson es un ajuste adecuado) | ||||
9. Conclusiones
Conclusiones
- Modelo: de las distribuciones evaluadas (Bernoulli, Binomial, Geométrica, Poisson), Poisson es la que corresponde a X (conteo de eventos independientes en un intervalo fijo).
- Estimación: λ̂ = 1.0796 pozos/semana (n = 314 semanas, 2019-12-29 a 2025-12-28).
- Prueba de Pearson: X² calc = 6.6726 vs. X² crítico = 9.4877 (α = 0.05, gl = 4) → el modelo Poisson se acepta.
- Sobre-dispersión: varianza (1.3515) = 1.25 × la media, la tasa no fue perfectamente constante en el periodo. Una Binomial Negativa podría ajustar mejor esta variabilidad.
- En la práctica: ~1.1 pozos/semana en promedio; P(X = 0) ≈ 33.97%; P(X ≥ 4) ≈ 2.43%.