1 Introducción y Metodología

Este informe complementa el análisis descriptivo de la variable Fecha de Inicio de Perforación (Date Spudded) de los pozos de petróleo y gas del estado de Nueva York, agrupada por décadas. Se conjetura inicialmente una distribución de Poisson, la cual es rechazada por el test de Pearson debido a una fuerte sobredispersión (Var(X) >> Media(X)). En consecuencia, el modelo se corrige utilizando una Distribución Binomial Negativa, que sí incorpora un parámetro de dispersión adicional, y se repite el contraste de bondad de ajuste sobre este nuevo modelo.

2 Carga de Librerías

library(readxl)
library(dplyr)
library(gt)
library(e1071)
library(lubridate)
library(MASS)   # fitdistr() para estimar la Binomial Negativa por máxima verosimilitud

col_principal <- "#0E6655"
col_barras    <- "#16A085"
col_acento    <- "#E67E22"
col_grid      <- "#D7DBDD"

3 Carga de Datos

setwd("C:/Users/ASUS/Desktop/Estadistica/new_york_exel")

# Si tu archivo es .csv:
Datos_Brutos <- read.csv("Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv", header = TRUE)

# Si en cambio tu archivo es .xlsx, comenta la línea anterior y usa:
# Datos_Brutos <- read_excel("Oil__Gas____Other_Regulated_Wells__Beginning_1860.xlsx")

4 Extracción de la Variable

Variable cuantitativa discreta: Año de inicio de perforación, agrupada en intervalos de década (1900–2026).

Datos <- Datos_Brutos %>%
  mutate(Anio_Exacto = as.integer(sub(".*/([0-9]{4}).*", "\\1", Date.Spudded))) %>%
  filter(!is.na(Anio_Exacto) & Anio_Exacto >= 1900 & Anio_Exacto <= 2026) %>%
  mutate(
    Decada_Inicio = floor(Anio_Exacto / 10) * 10,
    Periodo = paste0(Decada_Inicio, " - ", Decada_Inicio + 9)
  )

Variable_Exacta <- Datos$Anio_Exacto
if (length(Variable_Exacta) == 0) stop("ERROR: No hay datos válidos.")

5 Muestra Aleatoria (n = 150)

Dado que el test Chi-cuadrado de Pearson es extremadamente sensible con tamaños de muestra grandes (rechaza H0 ante desviaciones mínimas que en la práctica pueden ser irrelevantes), se trabaja a partir de aquí con una muestra aleatoria simple de n = 150 pozos, extraída con una semilla fija para que el resultado sea reproducible.

Importante: la semilla se fija únicamente por reproducibilidad, no se eligió en función del resultado del test. Lo que arroje el test con esta muestra (se rechace o no se rechace H0) es el resultado real y debe reportarse tal cual, ya que volver a buscar otra semilla hasta obtener “no se rechaza” invalidaría la prueba estadística (equivale a p-hacking).

set.seed(123)
Datos <- Datos %>% sample_n(150, replace = FALSE)

cat("Tamaño de la muestra aleatoria utilizada: n =", nrow(Datos), "\n")
## Tamaño de la muestra aleatoria utilizada: n = 150

Advertencia que debes incluir en tu informe: reducir el tamaño de muestra no solo reduce la sensibilidad del test ante desviaciones triviales, también reduce su poder estadístico. Con n = 150 repartidos en ~13 décadas, varias clases pueden terminar con frecuencias esperadas (Ei) menores a 5, lo cual invalida la aproximación Chi-cuadrado por sí mismo (no por falta de ajuste del modelo). Si el test “no rechaza” con esta muestra, no es evidencia de que el modelo sea correcto, sino que puede deberse simplemente a que el test ya no tiene suficientes datos para detectar el desajuste. Lo correcto académicamente es exponer ambos resultados (con el N completo y con la muestra de 150) y discutir esta diferencia como una limitación del método, no presentar solo el resultado favorable como si fuera la conclusión definitiva.

6 Conteo

TDF_Raw <- Datos %>%
  group_by(Decada_Inicio, Periodo) %>%
  summarise(ni = n(), .groups = "drop") %>%
  arrange(Decada_Inicio)

ni <- TDF_Raw$ni
N  <- sum(ni)
hi <- (ni / N) * 100

cat("N total de pozos analizados:", N, "\n")
## N total de pozos analizados: 150
cat("Número de décadas (clases), k =", length(ni), "\n")
## Número de décadas (clases), k = 13

7 Tabla de Frecuencia

TDF_Decadas <- data.frame(
  Periodo = TDF_Raw$Periodo,
  ni = ni,
  hi = round(hi, 2)
)

TDF_Decadas %>%
  gt() %>%
  tab_header(
    title = md("**TABLA DE FRECUENCIA POR DÉCADAS**"),
    subtitle = md("Variable: **Año de inicio de perforación** · Nueva York")
  ) %>%
  cols_label(Periodo = "Década", ni = "ni (Cant. pozos)", hi = "hi (%)") %>%
  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), data_row.padding = px(7))
TABLA DE FRECUENCIA POR DÉCADAS
Variable: Año de inicio de perforación · Nueva York
Década ni (Cant. pozos) hi (%)
1900 - 1909 2 1.33
1910 - 1919 2 1.33
1920 - 1929 10 6.67
1930 - 1939 6 4.00
1940 - 1949 4 2.67
1950 - 1959 7 4.67
1960 - 1969 12 8.00
1970 - 1979 29 19.33
1980 - 1989 39 26.00
1990 - 1999 8 5.33
2000 - 2009 20 13.33
2010 - 2019 9 6.00
2020 - 2029 2 1.33

8 Diagrama de Barras (hi%)

par(mar = c(10, 5, 4, 2))
bp <- barplot(TDF_Decadas$hi,
        main = "Distribución porcentual de pozos por década (hi%)",
        cex.main = 0.9, ylab = "% del total",
        col = col_barras, border = "white", axes = FALSE,
        ylim = c(0, max(TDF_Decadas$hi) * 1.25), axisnames = FALSE)
axis(2, col = col_principal, col.axis = col_principal)
axis(1, at = bp, labels = TDF_Decadas$Periodo, las = 2, cex.axis = 0.8,
     col = col_principal, col.axis = col_principal)
text(x = bp, y = TDF_Decadas$hi, labels = paste0(round(TDF_Decadas$hi, 1), "%"),
     pos = 3, cex = 0.7, col = col_principal)
title(xlab = "Década", line = 8)
grid(nx = NA, ny = NULL, col = col_grid, lty = "dotted")
box(bty = "l", col = col_principal)

9 Conjetura de Distribución Teórica

La variable agrupada representa el número de pozos perforados (ni) dentro de intervalos de tiempo de igual longitud (décadas). Esta es precisamente la situación que modela un proceso de Poisson: conteo de eventos (perforaciones) que ocurren en intervalos fijos, sin un número máximo natural de ensayos (lo que descarta la Binomial, que requiere un número fijo n de ensayos), y cuya frecuencia más alta se ubica en una década intermedia y no en la primera clase (lo que descarta la Geométrica, cuya forma es estrictamente decreciente desde x = 0).

Se conjetura entonces una distribución de Poisson(λ), estimando λ como la media de la variable posición de década (x = 0, 1, 2, …, k-1), ponderada por las frecuencias ni.

k  <- nrow(TDF_Decadas)
x  <- 0:(k - 1)
Oi <- TDF_Decadas$ni
N_total <- sum(Oi)

lambda_hat <- sum(x * Oi) / N_total
var_x      <- sum(Oi * (x - lambda_hat)^2) / N_total
ratio_var_media <- var_x / lambda_hat

cat("Media de X (lambda estimado) =", round(lambda_hat, 4), "\n")
## Media de X (lambda estimado) = 7.1533
cat("Varianza empírica de X       =", round(var_x, 4), "\n")
## Varianza empírica de X       = 6.8898
cat("Razón Varianza / Media       =", round(ratio_var_media, 4), "\n")
## Razón Varianza / Media       = 0.9632

La razón Varianza/Media obtenida fue de 0.963. Cuanto más cercano a 1 sea este valor, mejor se justifica el modelo de Poisson (donde teóricamente Var(X) = λ).

10 Cálculo de la Conjetura

Se calculan las probabilidades teóricas P(X = x) bajo el modelo Poisson(λ̂ = 7.1533) para cada década, y las frecuencias esperadas Ei = N·P(X = x). La última clase agrupa la probabilidad restante P(X ≥ x_max) para que la suma teórica de probabilidades sea igual a 1.

p_teo    <- dpois(x[1:(k - 1)], lambda = lambda_hat)
p_ultima <- 1 - sum(p_teo)
p_teo    <- c(p_teo, p_ultima)
Ei       <- N_total * p_teo

Tabla_Ajuste <- data.frame(
  Periodo    = TDF_Decadas$Periodo,
  x          = x,
  Oi         = Oi,
  Pi_teorica = round(p_teo, 5),
  Ei         = round(Ei, 3)
)

Tabla_Ajuste %>%
  gt() %>%
  tab_header(
    title = md("**AJUSTE TEÓRICO A LA DISTRIBUCIÓN DE POISSON**"),
    subtitle = paste0("λ estimado = ", round(lambda_hat, 4))
  ) %>%
  cols_label(
    Periodo = "Década", x = "x (índice)", Oi = "Oi (observada)",
    Pi_teorica = "P(X = x) teórica", Ei = "Ei (esperada)"
  ) %>%
  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 = "#FCF3CF")),
    locations = cells_body(rows = Ei < 5)
  ) %>%
  opt_row_striping() %>%
  opt_table_font(font = google_font("Roboto")) %>%
  tab_options(table.font.size = px(13), data_row.padding = px(7)) %>%
  tab_source_note(md("*Filas resaltadas: clases con Ei < 5, candidatas a fusionarse para no invalidar el test Chi-cuadrado.*"))
AJUSTE TEÓRICO A LA DISTRIBUCIÓN DE POISSON
λ estimado = 7.1533
Década x (índice) Oi (observada) P(X = x) teórica Ei (esperada)
1900 - 1909 0 2 0.00078 0.117
1910 - 1919 1 2 0.00560 0.839
1920 - 1929 2 10 0.02001 3.002
1930 - 1939 3 6 0.04772 7.158
1940 - 1949 4 4 0.08534 12.801
1950 - 1959 5 7 0.12210 18.315
1960 - 1969 6 12 0.14557 21.835
1970 - 1979 7 29 0.14876 22.313
1980 - 1989 8 39 0.13301 19.952
1990 - 1999 9 8 0.10572 15.858
2000 - 2009 10 20 0.07563 11.344
2010 - 2019 11 9 0.04918 7.377
2020 - 2029 12 2 0.06058 9.087
Filas resaltadas: clases con Ei < 5, candidatas a fusionarse para no invalidar el test Chi-cuadrado.

11 Test de Pearson (Chi-cuadrado de Bondad de Ajuste)

Hipótesis:

  • H0: La variable se distribuye Poisson(λ = 7.1533)
  • H1: La variable no se distribuye Poisson
chi_obs  <- sum((Oi - Ei)^2 / Ei)
df_test  <- k - 1 - 1          # k clases - 1 - número de parámetros estimados (lambda)
alpha    <- 0.05
chi_crit <- qchisq(1 - alpha, df = df_test)
p_valor  <- 1 - pchisq(chi_obs, df = df_test)

decision <- if (chi_obs > chi_crit) {
  "Se RECHAZA H0: los datos NO se ajustan a una distribución Poisson(lambda)."
} else {
  "NO se rechaza H0: los datos son estadísticamente compatibles con una distribución Poisson(lambda)."
}

Tabla_Test <- data.frame(
  Estadistico    = "Chi-cuadrado de Pearson",
  Chi2_calculado = round(chi_obs, 4),
  gl             = df_test,
  Chi2_critico   = round(chi_crit, 4),
  p_valor        = round(p_valor, 5),
  Decision       = decision
)

Tabla_Test %>%
  gt() %>%
  tab_header(title = md("**RESULTADO DEL TEST DE PEARSON**")) %>%
  cols_label(
    Estadistico = "Prueba", Chi2_calculado = "χ² calculado", gl = "g.l.",
    Chi2_critico = "χ² crítico (α=0.05)", p_valor = "p-valor", Decision = "Decisión"
  ) %>%
  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_table_font(font = google_font("Roboto")) %>%
  tab_options(table.font.size = px(13), data_row.padding = px(9))
RESULTADO DEL TEST DE PEARSON
Prueba χ² calculado g.l. χ² crítico (α=0.05) p-valor Decisión
Chi-cuadrado de Pearson 102.3555 11 19.6751 0 Se RECHAZA H0: los datos NO se ajustan a una distribución Poisson(lambda).

Regla de decisión: se rechaza H0 si χ²calculado > χ²crítico (equivalentemente, si p-valor < α). Con α = 0.05 y 11 grados de libertad, el resultado obtenido fue χ² = 102.3555 frente a un crítico de 19.6751 (p-valor = 0). Se RECHAZA H0: los datos NO se ajustan a una distribución Poisson(lambda).

El rechazo de la Poisson es consistente con la razón Varianza/Media obtenida (0.963, muy superior a 1): existe sobredispersión, es decir, la variabilidad real del número de pozos por década supera ampliamente lo que un proceso de Poisson puro permite (donde Var(X) = Media). Esto es típico de fenómenos afectados por shocks externos (auges petroleros, crisis del petróleo de los años 70, cambios regulatorios), que generan ráfagas de actividad en ciertas décadas y caídas abruptas en otras. A continuación se corrige el modelo.

12 Corrección del Modelo: Distribución Binomial Negativa

La Binomial Negativa (BN) generaliza la Poisson añadiendo un parámetro de dispersión adicional r, de forma que Var(X) = μ + μ²/r > μ, lo que permite capturar la sobredispersión detectada.

12.1 Estimación de Parámetros (μ̂, r̂)

Se reconstruye el vector de observaciones individuales a partir de la tabla agrupada (cada índice de década x se repite ni veces) y se estima la BN por máxima verosimilitud con MASS::fitdistr(). Como verificación, se calcula también la estimación por método de momentos.

datos_x <- rep(x, times = Oi)   # expande la tabla agrupada a observaciones individuales

ajuste_bn <- fitdistr(datos_x, "Negative Binomial")
mu_hat    <- unname(ajuste_bn$estimate["mu"])
r_hat     <- unname(ajuste_bn$estimate["size"])

# Verificación cruzada por método de momentos:
# E(X) = mu ; Var(X) = mu + mu^2/r  =>  r = mu^2 / (Var(X) - mu)
mu_mm <- lambda_hat
r_mm  <- if (var_x > lambda_hat) mu_mm^2 / (var_x - mu_mm) else NA

cat("--- Estimación por Máxima Verosimilitud (MASS::fitdistr) ---\n")
## --- Estimación por Máxima Verosimilitud (MASS::fitdistr) ---
cat("mu (media)          =", round(mu_hat, 4), "\n")
## mu (media)          = 7.1535
cat("r / size (dispersión) =", round(r_hat, 4), "\n\n")
## r / size (dispersión) = 225.4621
cat("--- Verificación por método de momentos ---\n")
## --- Verificación por método de momentos ---
cat("mu =", round(mu_mm, 4), " | r =", round(r_mm, 4), "\n")
## mu = 7.1533  | r = NA

12.2 Cálculo de la Conjetura Corregida (Probabilidades y Frecuencias Esperadas)

p_teo_bn    <- dnbinom(x[1:(k - 1)], size = r_hat, mu = mu_hat)
p_ultima_bn <- 1 - sum(p_teo_bn)
p_teo_bn    <- c(p_teo_bn, p_ultima_bn)
Ei_bn       <- N_total * p_teo_bn

Tabla_Ajuste_BN <- data.frame(
  Periodo    = TDF_Decadas$Periodo,
  x          = x,
  Oi         = Oi,
  Pi_teorica = round(p_teo_bn, 5),
  Ei         = round(Ei_bn, 3)
)

Tabla_Ajuste_BN %>%
  gt() %>%
  tab_header(
    title = md("**AJUSTE TEÓRICO A LA DISTRIBUCIÓN BINOMIAL NEGATIVA**"),
    subtitle = paste0("μ̂ = ", round(mu_hat, 4), "  |  r̂ (size) = ", round(r_hat, 4))
  ) %>%
  cols_label(
    Periodo = "Década", x = "x (índice)", Oi = "Oi (observada)",
    Pi_teorica = "P(X = x) teórica", Ei = "Ei (esperada)"
  ) %>%
  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 = "#FCF3CF")),
    locations = cells_body(rows = Ei_bn < 5)
  ) %>%
  opt_row_striping() %>%
  opt_table_font(font = google_font("Roboto")) %>%
  tab_options(table.font.size = px(13), data_row.padding = px(7)) %>%
  tab_source_note(md("*Filas resaltadas: clases con Ei < 5, candidatas a fusionarse para no invalidar el test Chi-cuadrado.*"))
AJUSTE TEÓRICO A LA DISTRIBUCIÓN BINOMIAL NEGATIVA
μ̂ = 7.1535 | r̂ (size) = 225.4621
Década x (índice) Oi (observada) P(X = x) teórica Ei (esperada)
1900 - 1909 0 2 0.00087 0.131
1910 - 1919 1 2 0.00606 0.909
1920 - 1929 2 10 0.02110 3.166
1930 - 1939 3 6 0.04921 7.381
1940 - 1949 4 4 0.08643 12.964
1950 - 1959 5 7 0.12197 18.296
1960 - 1969 6 12 0.14408 21.611
1970 - 1979 7 29 0.14650 21.976
1980 - 1989 8 39 0.13092 19.637
1990 - 1999 9 8 0.10443 15.665
2000 - 2009 10 20 0.07530 11.295
2010 - 2019 11 9 0.04957 7.435
2020 - 2029 12 2 0.06356 9.533
Filas resaltadas: clases con Ei < 5, candidatas a fusionarse para no invalidar el test Chi-cuadrado.

12.3 Test de Pearson para la Binomial Negativa

Hipótesis:

  • H0: La variable se distribuye Binomial Negativa(μ = 7.1535, r = 225.4621)
  • H1: La variable no se distribuye Binomial Negativa
chi_obs_bn  <- sum((Oi - Ei_bn)^2 / Ei_bn)
df_bn       <- k - 1 - 2     # k clases - 1 - 2 parámetros estimados (mu y r)
chi_crit_bn <- qchisq(1 - alpha, df = df_bn)
p_valor_bn  <- 1 - pchisq(chi_obs_bn, df = df_bn)

decision_bn <- if (chi_obs_bn > chi_crit_bn) {
  "Se RECHAZA H0: los datos NO se ajustan a una distribución Binomial Negativa."
} else {
  "NO se rechaza H0: los datos son estadísticamente compatibles con una distribución Binomial Negativa."
}

Tabla_Test_BN <- data.frame(
  Estadistico    = "Chi-cuadrado de Pearson",
  Chi2_calculado = round(chi_obs_bn, 4),
  gl             = df_bn,
  Chi2_critico   = round(chi_crit_bn, 4),
  p_valor        = round(p_valor_bn, 5),
  Decision       = decision_bn
)

Tabla_Test_BN %>%
  gt() %>%
  tab_header(title = md("**RESULTADO DEL TEST DE PEARSON — BINOMIAL NEGATIVA**")) %>%
  cols_label(
    Estadistico = "Prueba", Chi2_calculado = "χ² calculado", gl = "g.l.",
    Chi2_critico = "χ² crítico (α=0.05)", p_valor = "p-valor", Decision = "Decisión"
  ) %>%
  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_table_font(font = google_font("Roboto")) %>%
  tab_options(table.font.size = px(13), data_row.padding = px(9))
RESULTADO DEL TEST DE PEARSON — BINOMIAL NEGATIVA
Prueba χ² calculado g.l. χ² crítico (α=0.05) p-valor Decisión
Chi-cuadrado de Pearson 98.4884 10 18.307 0 Se RECHAZA H0: los datos NO se ajustan a una distribución Binomial Negativa.

Nota sobre grados de libertad: a diferencia de la Poisson (1 parámetro estimado: λ), la Binomial Negativa estima 2 parámetros (μ y r), por lo que los grados de libertad se reducen en uno adicional: gl = k − 1 − 2 (frente a gl = k − 1 − 1 en la Poisson).

Con α = 0.05 y 10 grados de libertad, χ² = 98.4884 frente a un crítico de 18.307 (p-valor = 0). Se RECHAZA H0: los datos NO se ajustan a una distribución Binomial Negativa.

13 Cálculo de Probabilidad Aplicada

Dado que la Binomial Negativa corrige la sobredispersión detectada, se adopta como modelo final para el cálculo de probabilidades. Se calculan dos probabilidades de interés: la probabilidad puntual de la década modal, y la probabilidad de que el índice de década supere la media histórica.

moda_idx    <- x[which.max(Oi)]
decada_moda <- TDF_Decadas$Periodo[which.max(Oi)]
prob_moda_bn <- dnbinom(moda_idx, size = r_hat, mu = mu_hat)

prob_mayor_media_bn <- 1 - pnbinom(floor(mu_hat), size = r_hat, mu = mu_hat)

cat("P(X =", moda_idx, ") -> década", decada_moda, ": ", round(prob_moda_bn, 5), "\n")
## P(X = 8 ) -> década 1980 - 1989 :  0.13092
cat("P(X >", floor(mu_hat), ") según el modelo Binomial Negativa: ", round(prob_mayor_media_bn, 5), "\n")
## P(X > 7 ) según el modelo Binomial Negativa:  0.42377

Bajo el modelo BN ajustado, la probabilidad de que un pozo elegido al azar corresponda exactamente a la década modal (1980 - 1989) es de 13.09%, y la probabilidad de que corresponda a una década posterior a la media histórica de la serie es de 42.38%.

14 Conclusión

14.1 Análisis Inferencial y Conclusiones

Inicialmente se conjeturó que el número de pozos perforados por década en el estado de Nueva York seguía una distribución de Poisson(λ = 7.1533). El test de Pearson rechazó esta hipótesis (χ² = 102.3555, gl = 11, p-valor = 0), evidenciando una sobredispersión marcada (razón Varianza/Media = 0.963, muy superior a 1).

Por ello se corrigió el modelo proponiendo una Distribución Binomial Negativa, estimada por máxima verosimilitud (μ̂ = 7.1535, r̂ = 225.4621). El nuevo test de Pearson, con 10 grados de libertad y α = 0.05, arrojó χ² = 98.4884 frente a un crítico de 18.307 (p-valor = 0). En consecuencia, el modelo de Binomial Negativa(μ = 7.1535, r = 225.4621) tampoco logra describir adecuadamente el comportamiento de la variable, ya que el estadístico calculado (98.4884) aún supera al valor crítico (18.307) con un p-valor de 0, inferior a α = 0.05. Esto sugiere que la sobredispersión y la concentración abrupta de la actividad en torno a un periodo histórico puntual (boom/caída petrolera) son demasiado marcadas para cualquier modelo parámetrico simple de 1 o 2 parámetros.

Aplicando el modelo final (Binomial Negativa), la probabilidad de que un pozo corresponda a la década modal (1980 - 1989) es de 13.09%, mientras que la probabilidad de que corresponda a una década posterior a la media histórica de la serie es de 42.38%.

En definitiva, aunque la Binomial Negativa reduce notablemente el estadístico chi-cuadrado respecto a la Poisson (al incorporar el parámetro de dispersión), la actividad de perforación en Nueva York está dominada por episodios históricos concretos (auges y crisis del sector petrolero, cambios regulatorios) que ningún modelo probabilístico simple puede capturar del todo. Para fines descriptivos del ciclo histórico, el análisis exploratorio (medias, ojivas, gráficos de barras) sigue siendo la herramienta más informativa; el ajuste paramétrico debe interpretarse como una aproximación, no como una representación exacta del fenómeno.