1. Cargar Librerías

library(readr); library(dplyr); library(gt); library(MASS)
cat("Librerías cargadas correctamente.\n")

## Librerías cargadas correctamente.

2. Cargar Datos

ruta_csv <- file.choose()
datos    <- read_csv(ruta_csv, show_col_types = FALSE)
cat("Archivo:", basename(ruta_csv), "| Filas:", nrow(datos), "\n")

## Archivo: oil_and_gas_leases_data (2).csv | Filas: 47757

3. Conteo

# AVG_PRODUCTION = CUMULATIVE_PRODUCTION / YEARS_ACTIVE, es decir, la
# producción ANUAL promedio de cada pozo. Es una variable extremadamente
# asimétrica (rango de 0.75 a más de 800,000): se recorta el 1% superior
# de valores atípicos (criterio estándar para variables muy sesgadas),
# ya que sin ese recorte el histograma queda degenerado (casi todos los
# datos caen en un solo intervalo) y ningún ajuste es estadísticamente
# válido, aunque el Pearson calculado sobre el rango completo pudiera
# verse alto de forma engañosa.
prod_raw <- datos %>%
  mutate(PA = suppressWarnings(as.numeric(AVG_PRODUCTION))) %>%
  filter(!is.na(PA), PA > 0) %>%
  pull(PA)

cap_p99 <- quantile(prod_raw, 0.99, na.rm = TRUE)
x_raw   <- prod_raw[prod_raw <= cap_p99]

n_conteo  <- length(x_raw)
k_sturges <- ceiling(1 + 3.322 * log10(n_conteo))

cat("Observaciones válidas (sin recorte):", length(prod_raw), "\n")

## Observaciones válidas (sin recorte): 47757

cat("Percentil 99 (tope de recorte):", round(cap_p99, 2), "\n")

## Percentil 99 (tope de recorte): 49354.18

cat("Observaciones tras recorte de outliers:", n_conteo, "\n")

## Observaciones tras recorte de outliers: 47279

cat("Clases (Regla de Sturges):", k_sturges, "\n")

## Clases (Regla de Sturges): 17

cat("Mínimo:", round(min(x_raw), 4), "| Máximo:", round(max(x_raw), 4), "\n")

## Mínimo: 0.75 | Máximo: 49276.21

4. Tabla de Distribución de Frecuencias

intervalos_cut <- cut(x, breaks=breaks_vec, right=FALSE, include.lowest=TRUE)
freq_abs       <- as.integer(table(intervalos_cut))
hi_dec         <- freq_abs / n
Ni_asc         <- cumsum(freq_abs); Hi_asc  <- cumsum(hi_dec)
Ni_desc        <- n - c(0, head(Ni_asc,-1)); Hi_desc <- 1 - c(0, head(Hi_asc,-1))
etiq           <- paste0("[",format(round(lim_inf,2),big.mark=",")," - ",format(round(lim_sup,2),big.mark=","),")")
etiq[k]        <- paste0("[",format(round(lim_inf[k],2),big.mark=",")," - ",format(round(lim_sup[k],2),big.mark=","),"]")

bind_rows(
  data.frame(Intervalo=etiq, MC=round(mc,2), ni=freq_abs,
             hi_pct=round(hi_dec*100,2), hi_real=round(hi_dec,4),
             Ni_a=Ni_asc, Hi_a=round(Hi_asc,4),
             Ni_d=Ni_desc, Hi_d=round(Hi_desc,4), stringsAsFactors=FALSE),
  data.frame(Intervalo="TOTAL", MC=NA_real_, ni=sum(freq_abs),
             hi_pct=round(sum(hi_dec)*100,2), hi_real=round(sum(hi_dec),4),
             Ni_a=max(Ni_asc), Hi_a=round(max(Hi_asc),4),
             Ni_d=max(Ni_desc), Hi_d=round(max(Hi_desc),4), stringsAsFactors=FALSE)
) %>%
  gt() %>%
  tab_header(title=md("**Tabla N\u00b01: Distribución de Frecuencias**"),
             subtitle=md(paste0("*Producción Anual Promedio, Kansas (n=",format(n,big.mark=","),")*"))) %>%
  cols_label(Intervalo=md("**Intervalo**"), MC=md("**MC**"), ni=md("**ni**"),
             hi_pct=md("**hi%**"), hi_real=md("**hi**"),
             Ni_a=md("**Ni (asc)**"), Hi_a=md("**Hi (asc)**"),
             Ni_d=md("**Ni (desc)**"), Hi_d=md("**Hi (desc)**")) %>%
  tab_style(style=list(cell_fill(color="#2C2C2C"),cell_text(color="white",weight="bold")),
            locations=cells_column_labels()) %>%
  tab_style(style=cell_fill(color="#F5F5F5"), locations=cells_body(rows=seq(1,k+1,by=2))) %>%
  tab_style(style=list(cell_fill(color="#D6D6D6"),cell_text(weight="bold")),
            locations=cells_body(rows=Intervalo=="TOTAL")) %>%
  fmt_missing(columns=everything(), missing_text="-") %>%
  tab_source_note(source_note=md("*Autor: Fernando Almeida*")) %>%
  tab_options(table.width=pct(100), table.font.size=px(13), data_row.padding=px(6))

Intervalo	MC	ni	hi%	hi	Ni (asc)	Hi (asc)	Ni (desc)	Hi (desc)
Tabla N°1: Distribución de Frecuencias
Producción Anual Promedio, Kansas (n=47,279)
[ 0.75 - 4,928.30)	2464.52	31483	66.59	0.6659	31483	0.6659	47279	1.0000
[ 4,928.30 - 9,855.84)	7392.07	7018	14.84	0.1484	38501	0.8143	15796	0.3341
[ 9,855.84 - 14,783.39)	12319.62	3291	6.96	0.0696	41792	0.8839	8778	0.1857
[14,783.39 - 19,710.93)	17247.16	2067	4.37	0.0437	43859	0.9277	5487	0.1161
[19,710.93 - 24,638.48)	22174.71	1489	3.15	0.0315	45348	0.9592	3420	0.0723
[24,638.48 - 29,566.03)	27102.25	978	2.07	0.0207	46326	0.9798	1931	0.0408
[29,566.03 - 34,493.57)	32029.80	475	1.00	0.0100	46801	0.9899	953	0.0202
[34,493.57 - 39,421.12)	36957.35	222	0.47	0.0047	47023	0.9946	478	0.0101
[39,421.12 - 44,348.66)	41884.89	143	0.30	0.0030	47166	0.9976	256	0.0054
[44,348.66 - 49,276.21]	46812.44	113	0.24	0.0024	47279	1.0000	113	0.0024
TOTAL	-	47279	100.00	1.0000	47279	1.0000	47279	1.0000
Autor: Fernando Almeida

5. Gráfico Histograma General

grises        <- gray(seq(0.25, 0.80, length.out=k))
h_obj         <- hist(x, breaks=breaks_vec, plot=FALSE)
h_obj$density <- hi_dec

par(mar=c(5,7,6,2))
plot(h_obj, col=grises, border="black", freq=FALSE,
     main="", xlab="", ylab="", las=1, xaxt="n")
axis(1, at=breaks_vec, labels=format(round(breaks_vec,0),big.mark=","), las=2, cex.axis=0.75)
mtext("Densidad de Probabilidad", side=2, line=5, cex=1)
mtext("Producción Anual Promedio (bbl/año)", side=1, line=4.5, cex=1)
mtext("Histograma General - Producción Anual, Kansas, EE.UU.",
      side=3, line=3, cex=0.95, font=2)

6. Conjetura y Ajuste del Modelo

El histograma de AVG_PRODUCTION (producción anual promedio por pozo) muestra una fuerte asimetría a la derecha, típica de variables de producción petrolera: muchos pozos con producción baja o moderada y unos pocos con producción muy alta. Este comportamiento es característico de una distribución Log-Normal, por lo que se conjetura ese modelo, contrastándolo también contra Normal y Exponencial (ambos desplazados al mínimo cuando corresponde).

set.seed(42)

lbl <- c(normal = "Normal", lognormal = "Log-Normal", exponential = "Exponencial")

k_adaptativo <- function(n) {
  max(3, min(k, ceiling(1 + 3.322 * log10(n))))
}

calc_pearson_zona <- function(datos, modelo, fit, offset, lim_min, lim_max) {
  k_zona <- k_adaptativo(length(datos))
  brks   <- seq(lim_min, lim_max, length.out = k_zona + 1)
  hi_obs <- hist(datos, breaks = brks, plot = FALSE)$counts / length(datos)
  mc_z   <- (head(brks,-1) + tail(brks,-1)) / 2

  if (modelo == "lognormal") {
    mc_pos <- mc_z - offset; mc_pos[mc_pos <= 0] <- 1e-9
    hi_teo <- dlnorm(mc_pos, meanlog = fit$estimate["meanlog"],
                     sdlog = fit$estimate["sdlog"]) * diff(brks)
  } else if (modelo == "exponential") {
    mc_pos <- mc_z - offset; mc_pos[mc_pos <= 0] <- 1e-9
    hi_teo <- dexp(mc_pos, rate = fit$estimate["rate"]) * diff(brks)
  } else {
    hi_teo <- dnorm(mc_z, mean = fit$estimate["mean"], sd = fit$estimate["sd"]) * diff(brks)
  }
  hi_teo <- hi_teo / sum(hi_teo)
  round(cor(hi_obs, hi_teo) * 100, 2)
}

evaluar_modelo <- function(datos, n_samp = 350, lim_min, lim_max) {
  set.seed(42)
  n_s    <- min(n_samp, length(datos))
  samp   <- sample(datos, size = n_s, replace = FALSE)
  offset <- min(datos) - 0.001
  d_pos  <- datos - offset
  s_pos  <- samp  - offset

  resultados <- list()

  tryCatch({
    f  <- fitdistr(datos, "normal")
    ks <- ks.test(samp, "pnorm",
                  mean = f$estimate["mean"], sd = f$estimate["sd"])
    resultados[["normal"]] <- list(
      fit = f, pval = ks$p.value,
      pearson = calc_pearson_zona(datos, "normal", f, 0, lim_min, lim_max),
      offset = 0)
  }, error = function(e){})

  tryCatch({
    f  <- fitdistr(d_pos, "lognormal")
    ks <- ks.test(s_pos, "plnorm",
                  meanlog = f$estimate["meanlog"], sdlog = f$estimate["sdlog"])
    resultados[["lognormal"]] <- list(
      fit = f, pval = ks$p.value,
      pearson = calc_pearson_zona(datos, "lognormal", f, offset, lim_min, lim_max),
      offset = offset)
  }, error = function(e){})

  tryCatch({
    f  <- fitdistr(d_pos, "exponential")
    ks <- ks.test(s_pos, "pexp", rate = f$estimate["rate"])
    resultados[["exponential"]] <- list(
      fit = f, pval = ks$p.value,
      pearson = calc_pearson_zona(datos, "exponential", f, offset, lim_min, lim_max),
      offset = offset)
  }, error = function(e){})

  mejor <- NULL
  mejor_pearson <- -Inf
  for (nombre in names(resultados)) {
    valor_pearson <- resultados[[nombre]]$pearson
    if (!is.na(valor_pearson) && valor_pearson > mejor_pearson) {
      mejor_pearson <- valor_pearson
      mejor <- nombre
    }
  }
  c(resultados[[mejor]], list(modelo = mejor))
}

validar <- function(datos, res, lim_min, lim_max) {
  set.seed(42)
  samp   <- sample(datos, size = min(50, length(datos)), replace = FALSE)
  offset <- res$offset

  pearson <- calc_pearson_zona(datos, res$modelo, res$fit, offset, lim_min, lim_max)

  if (res$modelo == "lognormal") {
    ks_res <- ks.test(samp - offset, "plnorm",
                      meanlog = res$fit$estimate["meanlog"],
                      sdlog   = res$fit$estimate["sdlog"])
  } else if (res$modelo == "exponential") {
    ks_res <- ks.test(samp - offset, "pexp", rate = res$fit$estimate["rate"])
  } else {
    ks_res <- ks.test(samp, "pnorm",
                      mean = res$fit$estimate["mean"],
                      sd   = res$fit$estimate["sd"])
  }

  pval_ks <- round(ks_res$p.value, 4)

  if (pearson > 70) {
    resultado_val <- "APROBADO"
  } else {
    resultado_val <- "RECHAZADO"
  }

  list(pearson = pearson, pval = pval_ks, val = resultado_val)
}

r_prod <- evaluar_modelo(x, lim_min = x_min, lim_max = x_max)
v_prod <- validar(x, r_prod, x_min, x_max)

cat("Modelo seleccionado:", lbl[r_prod$modelo], "\n")

## Modelo seleccionado: Log-Normal

cat("Parámetros:\n"); print(round(r_prod$fit$estimate, 6))

## Parámetros:

##  meanlog    sdlog 
## 7.754001 1.530154

cat("Pearson R%:", v_prod$pearson, "| KS p-valor:", v_prod$pval, "|", v_prod$val, "\n")

## Pearson R%: 99.94 | KS p-valor: 0.1649 | APROBADO

6.1 Histograma con Curva Teórica Ajustada

offset_prod <- r_prod$offset
xs          <- seq(x_min, x_max, length.out = 500)

if (r_prod$modelo == "lognormal") {
  xs_pos <- xs - offset_prod; xs_pos[xs_pos <= 0] <- 1e-9
  ys <- dlnorm(xs_pos, meanlog = r_prod$fit$estimate["meanlog"],
               sdlog = r_prod$fit$estimate["sdlog"])
} else if (r_prod$modelo == "exponential") {
  xs_pos <- xs - offset_prod; xs_pos[xs_pos <= 0] <- 1e-9
  ys <- dexp(xs_pos, rate = r_prod$fit$estimate["rate"])
} else {
  ys <- dnorm(xs, mean = r_prod$fit$estimate["mean"], sd = r_prod$fit$estimate["sd"])
}

par(mar=c(5,7,6,2))
plot(h_obj, col=grises, border="black", freq=FALSE,
     main="", xlab="", ylab="", las=1, xaxt="n",
     ylim=range(c(hi_dec, ys), na.rm=TRUE))
axis(1, at=breaks_vec, labels=format(round(breaks_vec,0),big.mark=","), las=2, cex.axis=0.75)
lines(xs, ys, col="black", lwd=2.5)
mtext("Densidad de Probabilidad", side=2, line=5, cex=1)
mtext("Producción Anual Promedio (bbl/año)", side=1, line=4.5, cex=1)
mtext(paste0("Ajuste ", lbl[r_prod$modelo], " - Producción Anual, Kansas, EE.UU."),
      side=3, line=3, cex=0.95, font=2)
legend("topright",
       legend = c("Histograma", paste0("Curva ", lbl[r_prod$modelo])),
       fill   = c("gray55", NA), border = c("black", NA),
       lty    = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.9)

7. Tabla de Validación

tabla_resumen <- data.frame(
  Zona       = "Producción Anual (recorte P99)",
  Modelo     = as.character(lbl[r_prod$modelo]),
  Pearson    = v_prod$pearson,
  KS_p       = v_prod$pval,
  Validacion = v_prod$val,
  stringsAsFactors = FALSE
)

tabla_resumen %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N\u00b02: Resumen de Validación**"),
    subtitle = md("*Pearson (R%) y Kolmogorov-Smirnov (p-valor)*")
  ) %>%
  cols_label(Zona=md("**Zona**"), Modelo=md("**Modelo Aplicado**"),
             Pearson=md("**Pearson (R %)**"), KS_p=md("**K-S (p-valor)**"),
             Validacion=md("**Validación**")) %>%
  tab_style(style=list(cell_fill(color="#2C2C2C"),cell_text(color="white",weight="bold")),
            locations=cells_column_labels()) %>%
  tab_style(style=list(cell_fill(color="#2C2C2C"),cell_text(color="white",weight="bold")),
            locations=cells_title(groups="title")) %>%
  tab_style(style=cell_text(color="darkgreen",weight="bold"),
            locations=cells_body(columns=Validacion, rows=Validacion=="APROBADO")) %>%
  tab_style(style=cell_text(color="darkred",weight="bold"),
            locations=cells_body(columns=Validacion, rows=Validacion=="RECHAZADO")) %>%
  tab_style(style=cell_borders(sides="bottom",color="#E0E0E0",weight=px(1)),
            locations=cells_body(rows=everything())) %>%
  cols_align(align="center", columns=c(Pearson,KS_p,Validacion)) %>%
  cols_align(align="left",   columns=c(Zona,Modelo)) %>%
  fmt_number(columns=Pearson, decimals=2) %>%
  fmt_number(columns=KS_p,    decimals=4) %>%
  tab_source_note(source_note=md("*Autor: Fernando Almeida*")) %>%
  tab_options(table.width=pct(90), table.font.size=px(13),
              heading.title.font.size=px(16), heading.subtitle.font.size=px(12),
              data_row.padding=px(6),
              column_labels.border.top.width=px(2),
              column_labels.border.bottom.width=px(2),
              table_body.border.bottom.width=px(2),
              table.border.top.style="hidden", table.border.bottom.style="hidden")

Zona	Modelo Aplicado	Pearson (R %)	K-S (p-valor)	Validación
Tabla N°2: Resumen de Validación
Pearson (R%) y Kolmogorov-Smirnov (p-valor)
Producción Anual (recorte P99)	Log-Normal	99.94	0.1649	APROBADO
Autor: Fernando Almeida

8. Intervalo de Confianza

El Intervalo de Confianza representa el puente fundamental entre el modelo empírico observado (Log-Normal) y la estimación poblacional. Aunque la distribución original de la Producción Anual presenta fuerte asimetría, el TLC garantiza que la distribución de las medias muestrales tenderá a la normalidad debido al volumen masivo de datos (\(n=\) 47,279).

Los postulados de confianza empírica sugieren:

\[P(\bar{x} - E < \mu < \bar{x} + E) \approx 68\%\]

\[P(\bar{x} - 2E < \mu < \bar{x} + 2E) \approx 95\%\]

\[P(\bar{x} - 3E < \mu < \bar{x} + 3E) \approx 99\%\]

Donde el Margen de Error (\(E\)) se define como:

\[E = \frac{\sigma}{\sqrt{n}}\]

media_muestral <- mean(x)
desv_est       <- sd(x)
n_total        <- length(x)
z_95           <- 1.96
error_est      <- desv_est / sqrt(n_total)
margen         <- z_95 * error_est
lim_inf_ic     <- media_muestral - margen
lim_sup_ic     <- media_muestral + margen

data.frame(
  Parametro      = "Producción Anual Promedio Kansas (bbl/año)",
  Lim_Inferior   = round(lim_inf_ic, 2),
  Media_Muestral = round(media_muestral, 2),
  Lim_Superior   = round(lim_sup_ic, 2),
  Error_Estandar = paste0("+/- ", round(margen, 2)),
  Confianza      = "95% (Z=1.96)",
  stringsAsFactors = FALSE
) %>%
  gt() %>%
  tab_header(
    title    = md("**TABLA N\u00b0 3: ESTIMACIÓN DE LA MEDIA POBLACIONAL**"),
    subtitle = md("*Inferencia Estadística para la Variable Producción Anual*")
  ) %>%
  cols_label(
    Parametro      = md("**Parámetro**"),
    Lim_Inferior   = md("**Lim_Inferior**"),
    Media_Muestral = md("**Media_Muestral**"),
    Lim_Superior   = md("**Lim_Superior**"),
    Error_Estandar = md("**Error_Estándar**"),
    Confianza      = md("**Confianza**")
  ) %>%
  tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
            locations = cells_column_labels()) %>%
  tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
            locations = cells_title(groups = "title")) %>%
  tab_style(style = list(cell_fill(color = "#C8E6C9"), cell_text(weight = "bold")),
            locations = cells_body(columns = Media_Muestral)) %>%
  tab_style(style = cell_borders(sides = "bottom", color = "#E0E0E0", weight = px(1)),
            locations = cells_body(rows = everything())) %>%
  cols_align(align = "center", columns = c(Lim_Inferior, Media_Muestral, Lim_Superior, Error_Estandar, Confianza)) %>%
  cols_align(align = "left",   columns = Parametro) %>%
  tab_source_note(source_note = md("*Autor: Fernando Almeida*")) %>%
  tab_options(table.width = pct(100), table.font.size = px(13),
              heading.title.font.size = px(16), heading.subtitle.font.size = px(12),
              data_row.padding = px(6),
              column_labels.border.top.width = px(2),
              column_labels.border.bottom.width = px(2),
              table_body.border.bottom.width = px(2),
              table.border.top.style = "hidden", table.border.bottom.style = "hidden")

Parámetro	Lim_Inferior	Media_Muestral	Lim_Superior	Error_Estándar	Confianza
TABLA N° 3: ESTIMACIÓN DE LA MEDIA POBLACIONAL
Inferencia Estadística para la Variable Producción Anual
Producción Anual Promedio Kansas (bbl/año)	5646.03	5714.8	5783.57	+/- 68.77	95% (Z=1.96)
Autor: Fernando Almeida

9. Conclusiones

Se trabajó con la variable AVG_PRODUCTION (producción anual promedio por pozo, calculada como CUMULATIVE_PRODUCTION / YEARS_ACTIVE), recortando el 1% superior de valores atípicos (tope en 49,354 bbl/año) para evitar un histograma degenerado por la fuerte cola derecha de la variable. Se evaluaron los modelos Normal, Log-Normal y Exponencial sobre el rango resultante (0.75 a 49,276 bbl/año), y se seleccionó el de mayor correlación de Pearson: Log-Normal (meanlog = 7.7540, sdlog = 1.5302). La prueba de Pearson obtuvo 99.94% y el p-valor de Kolmogorov-Smirnov fue 0.1649, por lo que el modelo quedó APROBADO.

Autor: Fernando Almeida — Análisis Estadístico, Kansas Hydrocarbon Leases Dataset

Producción Anual — Ajuste de Modelo

Fernando Almeida

2026-07-21