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

x_raw <- datos %>%
  mutate(LON = suppressWarnings(as.numeric(LONGITUDE))) %>%
  filter(!is.na(LON), LON >= -103.0, LON <= -94.0) %>%
  pull(LON)

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

cat("Observaciones válidas:", n_conteo, "\n")
## Observaciones válidas: 47757
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: -102.0439 | Máximo: -94.6179

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("[",round(lim_inf,4)," - ",round(lim_sup,4),")")
etiq[k]        <- paste0("[",round(lim_inf[k],4)," - ",round(lim_sup[k],4),"]")

bind_rows(
  data.frame(Intervalo=etiq, MC=round(mc,4), 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("*Longitud, 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))
Tabla N°1: Distribución de Frecuencias
Longitud, Kansas (n=47,757)
Intervalo MC ni hi% hi Ni (asc) Hi (asc) Ni (desc) Hi (desc)
[-102.0439 - -101.3013) -101.6726 3926 8.22 0.0822 3926 0.0822 47757 1.0000
[-101.3013 - -100.5587) -100.9300 7021 14.70 0.1470 10947 0.2292 43831 0.9178
[-100.5587 - -99.8161) -100.1874 3872 8.11 0.0811 14819 0.3103 36810 0.7708
[-99.8161 - -99.0735) -99.4448 7460 15.62 0.1562 22279 0.4665 32938 0.6897
[-99.0735 - -98.3309) -98.7022 9191 19.25 0.1925 31470 0.6590 25478 0.5335
[-98.3309 - -97.5883) -97.9596 3517 7.36 0.0736 34987 0.7326 16287 0.3410
[-97.5883 - -96.8457) -97.2170 2406 5.04 0.0504 37393 0.7830 12770 0.2674
[-96.8457 - -96.1031) -96.4744 1630 3.41 0.0341 39023 0.8171 10364 0.2170
[-96.1031 - -95.3605) -95.7318 6222 13.03 0.1303 45245 0.9474 8734 0.1829
[-95.3605 - -94.6179] -94.9892 2512 5.26 0.0526 47757 1.0000 2512 0.0526
TOTAL - 47757 100.00 1.0000 47757 1.0000 47757 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,6,6,2))
plot(h_obj, col=grises, border="black", freq=FALSE,
     main="", xlab="", ylab="", las=1, xaxt="n")
axis(1, at=breaks_vec, labels=round(breaks_vec,3), las=1, cex.axis=0.8)
mtext("Densidad de Probabilidad", side=2, line=4.5, cex=1)
mtext("Longitud (grados)",        side=1, line=3.5, cex=1)
mtext("Histograma General - Longitud, Kansas, EE.UU.",
      side=3, line=3, cex=0.95, font=2)

6. Conjetura

El histograma corresponde a una variable continua. Al no poder conjeturar ningún modelo inferencial se decidió trabajar por partes dividiendo el grafico.

set.seed(42)

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

# Numero de bins adaptado al tamano de la (sub)zona: evita que sub-zonas
# angostas hereden los mismos 10 bins del histograma completo, lo que
# dejaria muy pocos datos por bin y volveria el Pearson ruidoso/bajo.
k_adaptativo <- function(n) {
  max(3, min(k, ceiling(1 + 3.322 * log10(n))))
}

# Calcula el Pearson (correlacion hi_obs vs hi_teo) de un modelo ya ajustado
# en una zona [lim_min, lim_max], usando bins adaptados al tamano de esa zona.
# Se usa tanto para ELEGIR el modelo (evaluar_zona) como para APROBARLO
# (validar), de modo que ambos pasos midan lo mismo.
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)
}

# Funcion central: evalua Normal, Log-Normal y Exponencial con offset correcto
evaluar_zona <- function(datos, n_samp = 350, forzar = NULL, candidatos = NULL,
                          lim_min = NULL, lim_max = NULL) {
  if (length(datos) < 30) return(NULL)
  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

  # Limites de la zona para el calculo de Pearson: si no se pasan
  # explicitamente, se usa el rango real de los datos de la zona.
  z_min <- if (is.null(lim_min)) min(datos) else lim_min
  z_max <- if (is.null(lim_max)) max(datos) else lim_max

  resultados <- list()

  # Normal (datos originales, sin desplazamiento)
  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, z_min, z_max),
      aic = -2*f$loglik + 2*length(f$estimate), offset = 0)
  }, error = function(e){})

  # Log-Normal (datos desplazados: x - min_zona + 0.001)
  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, z_min, z_max),
      aic = -2*f$loglik + 2*length(f$estimate), offset = offset)
  }, error = function(e){})

  # Exponencial (datos desplazados: x - min_zona + 0.001)
  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, z_min, z_max),
      aic = -2*f$loglik + 2*length(f$estimate), offset = offset)
  }, error = function(e){})

  if (length(resultados) == 0) return(NULL)

  if (!is.null(forzar) && !is.null(resultados[[forzar]])) {
    mejor <- forzar
  } else {
    pool <- resultados
    if (!is.null(candidatos)) {
      pool_filtrado <- resultados[names(resultados) %in% candidatos]
      if (length(pool_filtrado) > 0) pool <- pool_filtrado
    }
    # Se elige por Pearson (mismo criterio con el que luego se aprueba/
    # rechaza en validar()), no por el p-valor de KS: son metricas
    # distintas y optimizar una no garantiza la otra.
    mejor <- NULL
    mejor_pearson <- -Inf
    for (nombre in names(pool)) {
      valor_pearson <- pool[[nombre]]$pearson
      if (!is.na(valor_pearson) && valor_pearson > mejor_pearson) {
        mejor_pearson <- valor_pearson
        mejor <- nombre
      }
    }
  }
  c(resultados[[mejor]], list(modelo = mejor))
}

# Ajusta un modelo forzado especifico (sin comparar con otros)
ajustar_forzado <- function(datos, modelo) {
  offset <- min(datos) - 0.001
  d_pos  <- datos - offset
  if (modelo == "normal") {
    f <- fitdistr(datos, "normal")
    list(fit = f, modelo = "normal", offset = 0,
         pval = ks.test(sample(datos, min(400, length(datos))),
                        "pnorm", mean = f$estimate["mean"],
                        sd = f$estimate["sd"])$p.value)
  } else if (modelo == "lognormal") {
    f <- fitdistr(d_pos, "lognormal")
    list(fit = f, modelo = "lognormal", offset = offset,
         pval = ks.test(sample(d_pos, min(400, length(d_pos))),
                        "plnorm", meanlog = f$estimate["meanlog"],
                        sdlog = f$estimate["sdlog"])$p.value)
  } else if (modelo == "exponential") {
    f <- fitdistr(d_pos, "exponential")
    list(fit = f, modelo = "exponential", offset = offset,
         pval = ks.test(sample(d_pos, min(400, length(d_pos))),
                        "pexp", rate = f$estimate["rate"])$p.value)
  }
}

# Valida un modelo (Pearson + KS) contra el histograma real en [lim_min, lim_max]
# Usa calc_pearson_zona() con bins adaptados al tamano de la zona -- el mismo
# criterio y los mismos bins con los que evaluar_zona() eligio el modelo --
# en lugar de los k=10 bins fijos del histograma general.
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)
}

# Cortes fijos para Oeste / Centro / Este
corte1 <- -99.810
corte2 <- -96.846
lim_oeste_min <- -102.044
lim_oeste_max <- corte1
lim_este_min  <- corte2
lim_este_max  <- -94.618

x_z1 <- x[x >= lim_oeste_min & x <  corte1]
x_z2 <- x[x >= corte1        & x <  corte2]
x_este_total <- x[x >= lim_este_min & x <= lim_este_max]

set.seed(42)
r2f <- ajustar_forzado(x_z2, "normal")

# --- Busqueda automatica de particion (1, 2 o 3 sub-zonas) ---
# Funcion generica (antes especifica de Zona Este): prueba distintos
# puntos de corte y se queda con la primera combinacion donde TODAS las
# sub-zonas resultantes aprueban (Pearson > 70%). Si ninguna combinacion
# aprueba del todo, devuelve la mejor encontrada (mayor Pearson minimo
# entre sub-zonas) y la marca como no aprobada. Se usa para Oeste y Este.
buscar_particion_zona <- function(datos, lim_min, lim_max, max_partes = 3, pts_extra = numeric(0)) {
  mejor_particion <- NULL
  mejor_score     <- -1000

  for (n_partes in 1:max_partes) {

    # Armamos la lista de combinaciones de cortes a probar
    lista_cortes <- list()
    if (n_partes == 1) {
      lista_cortes[[1]] <- numeric(0)
    } else {
      pts <- c(pts_extra, quantile(datos, probs = seq(0.2, 0.8, by = 0.1), names = FALSE))
      pts <- unique(pts)
      pts <- sort(pts[pts > lim_min & pts < lim_max])

      if (n_partes == 2) {
        for (i in 1:length(pts)) {
          lista_cortes[[length(lista_cortes) + 1]] <- pts[i]
        }
      } else {
        for (i in 1:(length(pts) - 1)) {
          for (j in (i + 1):length(pts)) {
            if (pts[j] - pts[i] > 0.05) {
              lista_cortes[[length(lista_cortes) + 1]] <- c(pts[i], pts[j])
            }
          }
        }
      }
    }

    # Probamos cada combinacion de cortes
    for (cortes in lista_cortes) {
      limites <- sort(unique(c(lim_min, cortes, lim_max)))
      if (length(limites) - 1 != n_partes) next

      partes <- list()
      combinacion_valida <- TRUE
      for (i in 1:n_partes) {
        d_sub <- datos[datos >= limites[i] & datos <= limites[i + 1]]
        if (length(d_sub) < 30) { combinacion_valida <- FALSE; break }
        r_sub <- evaluar_zona(d_sub, candidatos = c("normal", "lognormal", "exponential"),
                               lim_min = limites[i], lim_max = limites[i + 1])
        if (is.null(r_sub)) { combinacion_valida <- FALSE; break }
        v_sub <- validar(d_sub, r_sub, limites[i], limites[i + 1])
        partes[[i]] <- list(datos = d_sub, res = r_sub, val = v_sub,
                             lim_min = limites[i], lim_max = limites[i + 1])
      }
      if (!combinacion_valida) next

      pearson_minimo <- 1000
      todas_ok <- TRUE
      for (i in 1:n_partes) {
        if (partes[[i]]$val$pearson < pearson_minimo) pearson_minimo <- partes[[i]]$val$pearson
        if (partes[[i]]$val$val != "APROBADO") todas_ok <- FALSE
      }

      if (todas_ok) {
        return(list(n_partes = n_partes, limites = limites, partes = partes, todas_aprobadas = TRUE))
      }
      if (pearson_minimo > mejor_score) {
        mejor_score     <- pearson_minimo
        mejor_particion <- list(n_partes = n_partes, limites = limites, partes = partes, todas_aprobadas = FALSE)
      }
    }
  }
  mejor_particion
}

particion_oeste <- buscar_particion_zona(x_z1, lim_oeste_min, lim_oeste_max, max_partes = 3)
particion_este  <- buscar_particion_zona(x_este_total, lim_este_min, lim_este_max,
                                          max_partes = 3, pts_extra = -96.103)

reportar_particion <- function(particion, etiqueta) {
  estado <- if (particion$todas_aprobadas) "TODAS APROBADAS" else "mejor opcion encontrada (no todas aprobadas)"
  cat("Zona", etiqueta, "dividida en", particion$n_partes, "parte(s) -", estado, ". Cortes:",
      paste(round(particion$limites, 3), collapse = " | "), "\n")
  for (i in seq_along(particion$partes)) {
    p <- particion$partes[[i]]
    cat(" -", etiqueta, i, ":", round(p$lim_min,3), "a", round(p$lim_max,3),
        "| Modelo:", lbl[p$res$modelo], "| Pearson:", p$val$pearson,
        "| KS p:", p$val$pval, "|", p$val$val, "\n")
  }
}
reportar_particion(particion_este, "Este")
## Zona Este dividida en 1 parte(s) - TODAS APROBADAS . Cortes: -96.846 | -94.618 
##  - Este 1 : -96.846 a -94.618 | Modelo: Normal | Pearson: 90.44 | KS p: 0.0466 | APROBADO
# De las sub-zonas de Oeste, solo se reportan formalmente las que APRUEBAN
# (Pearson > 70%); las rechazadas se documentan en la nota desplegable.
partes_oeste_ok         <- list()
partes_oeste_rechazadas <- list()
for (i in 1:length(particion_oeste$partes)) {
  p <- particion_oeste$partes[[i]]
  if (p$val$val == "APROBADO") {
    partes_oeste_ok[[length(partes_oeste_ok) + 1]] <- p
  } else {
    partes_oeste_rechazadas[[length(partes_oeste_rechazadas) + 1]] <- p
  }
}

mejor_pearson_oeste_rechazada <- NA
if (length(partes_oeste_rechazadas) > 0) {
  mejor_pearson_oeste_rechazada <- partes_oeste_rechazadas[[1]]$val$pearson
  for (i in 1:length(partes_oeste_rechazadas)) {
    if (partes_oeste_rechazadas[[i]]$val$pearson > mejor_pearson_oeste_rechazada) {
      mejor_pearson_oeste_rechazada <- partes_oeste_rechazadas[[i]]$val$pearson
    }
  }
}

# Si solo queda una sub-zona aprobada se llama simplemente "Oeste"
# (sin numerar); si quedan varias, se numeran Oeste 1, Oeste 2, ...
etiqueta_oeste <- function(i) {
  if (length(partes_oeste_ok) == 1) {
    return("Oeste")
  }
  paste("Oeste", i)
}

6.1 Histograma con Cortes por Zona

for (i in seq_along(partes_oeste_ok)) {
  p <- partes_oeste_ok[[i]]
  cat("Zona", etiqueta_oeste(i), ":", length(p$datos), "obs | Modelo:", lbl[p$res$modelo],
      "| KS p =", p$val$pval, "\n")
}
## Zona Oeste : 10392 obs | Modelo: Log-Normal | KS p = 0.3559
cat("Zona Centro:", length(x_z2), "obs | Modelo:", lbl[r2f$modelo],
    "| KS p =", round(r2f$pval, 4), "\n")
## Zona Centro: 22545 obs | Modelo: Normal | KS p = 0.0187
for (i in seq_along(particion_este$partes)) {
  p <- particion_este$partes[[i]]
  cat("Zona Este", i, ":", length(p$datos), "obs | Modelo:", lbl[p$res$modelo],
      "| KS p =", p$val$pval, "\n")
}
## Zona Este 1 : 10365 obs | Modelo: Normal | KS p = 0.0466
h_obj         <- hist(x, breaks = breaks_vec, plot = FALSE)
h_obj$density <- hi_dec

par(mar = c(5,6,6,2))
plot(h_obj, col = grises, border = "black", freq = FALSE,
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = breaks_vec, labels = round(breaks_vec,3), las = 1, cex.axis = 0.8)
abline(v = corte1, col = "black", lty = 2, lwd = 2)
abline(v = corte2, col = "black", lty = 2, lwd = 2)

cortes_oeste <- particion_oeste$limites
cortes_este  <- particion_este$limites
if (length(cortes_oeste) > 2) {
  for (cc in cortes_oeste[2:(length(cortes_oeste) - 1)]) {
    abline(v = cc, col = "darkgreen", lty = 3, lwd = 2)
  }
}
if (length(cortes_este) > 2) {
  for (cc in cortes_este[2:(length(cortes_este) - 1)]) {
    abline(v = cc, col = "blue", lty = 3, lwd = 2)
  }
}

yt <- max(hi_dec)
for (i in seq_along(partes_oeste_ok)) {
  p <- partes_oeste_ok[[i]]
  text(mean(c(p$lim_min, p$lim_max)), yt*0.88,
       paste0(etiqueta_oeste(i), "\n(", round(p$lim_min,2), " a ", round(p$lim_max,2), ")\n", lbl[p$res$modelo]),
       cex=0.72, font=2)
}
text(mean(c(corte1, corte2)), yt*0.88,
     paste0("Centro\n(", corte1, " a ", corte2, ")\n", lbl[r2f$modelo]),
     cex=0.82, font=2)
for (i in seq_along(particion_este$partes)) {
  p <- particion_este$partes[[i]]
  text(mean(c(p$lim_min, p$lim_max)), yt*0.60,
       paste0("Este ", i, "\n(", round(p$lim_min,2), " a ", round(p$lim_max,2), ")\n", lbl[p$res$modelo]),
       cex=0.72, font=2)
}
mtext("Densidad de Probabilidad", side=2, line=4.5, cex=1)
mtext("Longitud (grados)",        side=1, line=3.5, cex=1)
mtext(paste0("Cortes: ", corte1, ", ", corte2,
             if (length(cortes_oeste) > 2) paste0(" | sub-cortes Oeste (",
               paste(round(cortes_oeste[2:(length(cortes_oeste)-1)],3), collapse=", "), ")") else "",
             if (length(cortes_este) > 2) paste0(" | sub-cortes Este (",
               paste(round(cortes_este[2:(length(cortes_este)-1)],3), collapse=", "), ")") else "",
             " - Kansas, EE.UU."),
      side=3, line=3, cex=0.88, font=2)

¿Por qué se excluye “Oeste 1” del ajuste formal? Una de las sub-zonas de Oeste no superó 70% de Pearson (mejor resultado: 66.16%), por lo que se excluye; la otra sub-zona sí aprueba y se reporta como Oeste.

6.2 Histogramas Individuales por Zona

plot_zona <- function(datos, res, titulo, pal, lim_min, lim_max) {
  offset <- res$offset
  n_z    <- length(datos)

  # Mismos intervalos (ancho c_amp) que el histograma principal (seccion 5.1)
  n_bins_z   <- ceiling((lim_max - lim_min) / c_amp)
  lim_inf_z  <- lim_min + (0:(n_bins_z - 1)) * c_amp
  lim_sup_z  <- lim_inf_z + c_amp
  lim_sup_z[n_bins_z] <- lim_max
  brks <- c(lim_inf_z, lim_sup_z[n_bins_z])

  h_z    <- hist(datos, breaks = brks, plot = FALSE)
  hi_z   <- h_z$counts / n_z                # hi: frecuencia relativa (diapo 17)
  dens_z <- hi_z / diff(brks)                # densidad real = hi / amplitud (diapo 32, 37)
  h_z$density <- dens_z
  xs     <- seq(lim_min, lim_max, length.out = 500)   # X: valores de Longitud
  nb     <- length(brks) - 1

  par(mar = c(5,6,6,2))
  plot(h_z, col = pal[seq_len(nb)], border = "black", freq = FALSE,
       main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
  axis(1, at = brks, labels = round(brks,3), las = 1, cex.axis = 0.8)

  # Y: f(x) calculada con la funcion de densidad teorica de la diapositiva
  # correspondiente -- Normal (diapo 57), Log-Normal (diapo 60) o Exponencial (diapo 63)
  if (res$modelo == "lognormal") {
    xs_pos <- xs - offset; xs_pos[xs_pos <= 0] <- 1e-9
    ys <- dlnorm(xs_pos,
                 meanlog = res$fit$estimate["meanlog"],
                 sdlog   = res$fit$estimate["sdlog"])
  } else if (res$modelo == "exponential") {
    xs_pos <- xs - offset; xs_pos[xs_pos <= 0] <- 1e-9
    ys <- dexp(xs_pos, rate = res$fit$estimate["rate"])
  } else {
    ys <- dnorm(xs,
                mean = res$fit$estimate["mean"],
                sd   = res$fit$estimate["sd"])
  }
  lines(xs, ys, col = "black", lwd = 2.5)

  mtext("Densidad de Probabilidad", side=2, line=4.5, cex=1)
  mtext("Longitud (grados)",        side=1, line=3.5, cex=1)
  mtext(titulo, side=3, line=3, cex=0.95, font=2)
  legend("topright",
         legend = c("Histograma", paste0("Curva ", lbl[res$modelo])),
         fill   = c("gray55", NA), border = c("black", NA),
         lty    = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.85)
}

par(mfrow = c(length(partes_oeste_ok) + 1 + particion_este$n_partes, 1))
for (i in seq_along(partes_oeste_ok)) {
  p <- partes_oeste_ok[[i]]
  plot_zona(p$datos, p$res,
            paste0("Zona ", etiqueta_oeste(i), " (", round(p$lim_min,3), " a ", round(p$lim_max,3), ") - ", lbl[p$res$modelo]),
            gray(seq(0.20,0.65,length.out=k)), p$lim_min, p$lim_max)
}
plot_zona(x_z2, r2f,
          paste0("Zona Centro (", corte1, " a ", corte2, ") - ", lbl[r2f$modelo]),
          gray(seq(0.30,0.75,length.out=k)), corte1, corte2)
for (i in seq_along(particion_este$partes)) {
  p <- particion_este$partes[[i]]
  plot_zona(p$datos, p$res,
            paste0("Zona Este ", i, " (", round(p$lim_min,3), " a ", round(p$lim_max,3), ") - ", lbl[p$res$modelo]),
            gray(seq(0.40,0.85,length.out=k)), p$lim_min, p$lim_max)
}

par(mfrow = c(1,1))

7. Cálculo de Parámetros y Probabilidades

v2 <- validar(x_z2, r2f, corte1, corte2)

cat("\n--- [", lbl[r2f$modelo], "] Zona Centro - Parametros ---\n")
## 
## --- [ Normal ] Zona Centro - Parametros ---
print(round(r2f$fit$estimate, 6))
##      mean        sd 
## -98.67543   0.73205
cat("Pearson R%:", v2$pearson, "| KS p-valor:", v2$pval, "\n")
## Pearson R%: 79.44 | KS p-valor: 0.4832
filas_zona <- function(particion, etiqueta) {
  filas <- list()
  for (i in 1:length(particion$partes)) {
    p <- particion$partes[[i]]
    cat("\n--- [", lbl[p$res$modelo], "] Zona", etiqueta, i, "- Parametros ---\n")
    print(round(p$res$fit$estimate, 6))
    cat("Pearson R%:", p$val$pearson, "| KS p-valor:", p$val$pval, "\n")
    fila <- data.frame(
      Zona = sprintf("%s %d (%.3f a %.3f)", etiqueta, i, p$lim_min, p$lim_max),
      Modelo = as.character(lbl[p$res$modelo]),
      Pearson = p$val$pearson,
      KS_p = p$val$pval,
      Validacion = p$val$val,
      stringsAsFactors = FALSE
    )
    filas[[i]] <- fila
  }
  filas
}

filas_oeste <- list()
for (i in 1:length(partes_oeste_ok)) {
  p <- partes_oeste_ok[[i]]
  cat("\n--- [", lbl[p$res$modelo], "] Zona", etiqueta_oeste(i), "- Parametros ---\n")
  print(round(p$res$fit$estimate, 6))
  cat("Pearson R%:", p$val$pearson, "| KS p-valor:", p$val$pval, "\n")
  fila <- data.frame(
    Zona = sprintf("%s (%.3f a %.3f)", etiqueta_oeste(i), p$lim_min, p$lim_max),
    Modelo = as.character(lbl[p$res$modelo]),
    Pearson = p$val$pearson,
    KS_p = p$val$pval,
    Validacion = p$val$val,
    stringsAsFactors = FALSE
  )
  filas_oeste[[i]] <- fila
}
## 
## --- [ Log-Normal ] Zona Oeste - Parametros ---
##   meanlog     sdlog 
## -0.901152  0.978237 
## Pearson R%: 87.05 | KS p-valor: 0.3559
filas_este <- filas_zona(particion_este, "Este")
## 
## --- [ Normal ] Zona Este 1 - Parametros ---
##       mean         sd 
## -95.653433   0.455534 
## Pearson R%: 90.44 | KS p-valor: 0.0466
tabla_resumen <- bind_rows(
  bind_rows(filas_oeste),
  data.frame(Zona = "Centro (-99.810 a -96.846)", Modelo = as.character(lbl[r2f$modelo]),
             Pearson = v2$pearson, KS_p = v2$pval, Validacion = v2$val, stringsAsFactors = FALSE),
  bind_rows(filas_este)
)

n_filas <- nrow(tabla_resumen)

tabla_resumen %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N\u00b02: Resumen de Validación por Zona**"),
    subtitle = md("*Pearson (R%) y Kolmogorov-Smirnov (p-valor) por zona*")
  ) %>%
  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_fill(color="#F5F5F5"), locations=cells_body(rows=seq(1, n_filas, by=2))) %>%
  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")
Tabla N°2: Resumen de Validación por Zona
Pearson (R%) y Kolmogorov-Smirnov (p-valor) por zona
Zona Modelo Aplicado Pearson (R %) K-S (p-valor) Validación
Oeste (-101.211 a -99.810) Log-Normal 87.05 0.3559 APROBADO
Centro (-99.810 a -96.846) Normal 79.44 0.4832 APROBADO
Este 1 (-96.846 a -94.618) Normal 90.44 0.0466 APROBADO
Autor: Fernando Almeida

8. Intervalo de Confianza

El Intervalo de Confianza representa el puente fundamental entre los modelos empíricos observados (Normal, Log-Normal y Exponencial) y la estimación poblacional. Aunque la distribución original de la Longitud presenta asimetría y comportamiento diferenciado por zona, el TLC garantiza que la distribución de las medias muestrales tenderá a la normalidad debido al volumen masivo de datos (\(n=\) 47,757).

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      = "Longitud Promedio Kansas (\u00b0)",
  Lim_Inferior   = round(lim_inf_ic, 3),
  Media_Muestral = round(media_muestral, 3),
  Lim_Superior   = round(lim_sup_ic, 3),
  Error_Estandar = paste0("+/- ", round(margen, 4)),
  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 Longitud*")
  ) %>%
  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")
TABLA N° 3: ESTIMACIÓN DE LA MEDIA POBLACIONAL
Inferencia Estadística para la Variable Longitud
Parámetro Lim_Inferior Media_Muestral Lim_Superior Error_Estándar Confianza
Longitud Promedio Kansas (°) -98.737 -98.719 -98.701 +/- 0.0178 95% (Z=1.96)
Autor: Fernando Almeida

9. Conclusiones

Se trabajó con la variable LONGITUDE dividida en 3 zona(s) geográficas con cortes fijos en −99.810° y −96.846°, excluyendo del ajuste la sub-zona de Oeste que no superó Pearson (ver nota en la Sección 6). Se aplicó el modelo Normal forzado en la zona Centro, y en las zonas Oeste y Este el modelo se eligió automáticamente por mayor correlación de Pearson en cada sub-zona. Los parámetros estimados fueron: Zona Oeste (meanlog = -0.9012, sdlog = 0.9782), Zona Centro (media = -98.6754, sd = 0.7321), Zona Este (media = -95.6534, sd = 0.4555). La prueba de Pearson obtuvo 87.05%, 79.44%, 90.44% respectivamente, y los p-valores Kolmogorov-Smirnov fueron 0.3559, 0.4832, 0.0466. Por tanto, 0 zona(s) no superaron la prueba KS.

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