1. Carga de Librerías

library(readr)
library(dplyr)
library(MASS)
if (!requireNamespace("gt", quietly = TRUE)) install.packages("gt", repos = "https://cloud.r-project.org")
library(gt)
cat("Librerías cargadas correctamente.\n")
## Librerías cargadas correctamente.

2. Carga de Dataset

datos <- read_csv(file.choose(), show_col_types = FALSE)
cat("Archivo cargado correctamente.\n")
## Archivo cargado correctamente.
cat("Total de registros (filas):", nrow(datos), "\n")
## Total de registros (filas): 47757

3. Extracción y Depuración de la Variable

Se extrae la variable cuantitativa discreta RANGE, que representa la división este-oeste del sistema de agrimensura rectangular (PLSS) de Kansas. Sus valores son números enteros en el rango 1–43. Se filtran únicamente los registros válidos dentro de ese rango.

x_raw <- datos %>%
  mutate(RNG = suppressWarnings(as.integer(RANGE))) %>%
  filter(!is.na(RNG), RNG >= 1, RNG <= 43) %>%
  pull(RNG)

n_unique <- length(unique(x_raw))
cat("Observaciones válidas (n):", length(x_raw), "\n")
## Observaciones válidas (n): 47757
cat("Valores únicos:", n_unique, "\n")
## Valores únicos: 43
cat("Como hay", n_unique, "> 10 valores únicos, se agrupa en intervalos de clase (máx. 10).\n")
## Como hay 43 > 10 valores únicos, se agrupa en intervalos de clase (máx. 10).

4. Tabla de Distribución de Frecuencias

Se aplica el criterio de máximo 10 intervalos para variables discretas agrupadas, con notación \([L_i, L_s)\) —cerrado a la izquierda, abierto a la derecha—, excepto el último intervalo que se cierra en ambos extremos.

\[k = 10 \quad \text{(máximo permitido)} \qquad c = \left\lceil \frac{\max - \min}{k} \right\rceil\]

intervalos_cut <- cut(x, breaks = breaks_vec, right = FALSE, include.lowest = TRUE)
freq_abs       <- as.integer(table(intervalos_cut))

li <- lim_inf
ls <- lim_sup

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_intervalo    <- paste0("[", li, " – ", ls, ")")
etiq_intervalo[k] <- paste0("[", li[k], " – ", ls[k] - 1, "]")

tabla_df <- data.frame(
  Intervalo = etiq_intervalo,
  MC        = as.integer(mc),
  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
)

total_row <- data.frame(
  Intervalo = "TOTAL", MC = NA_integer_,
  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
)

bind_rows(tabla_df, total_row) %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N°1: Distribución de Frecuencias**"),
    subtitle = md(paste0("*Variable Cuantitativa Discreta Agrupada: Range, ",
                         "arrendamientos de hidrocarburos, Kansas, EE.UU. (n = ",
                         format(n, big.mark = ","), " registros válidos)*"))
  ) %>%
  cols_label(
    Intervalo = md("**Intervalo [Li – Ls)**"), MC = md("**Marca de Clase**"),
    ni = md("**ni (FA)**"), hi_pct = md("**hi %**"), hi_real = md("**hi (decimal)**"),
    Ni_a = md("**Ni ↑ (FAAa)**"), Hi_a = md("**Hi ↑ (FRAa)**"),
    Ni_d = md("**Ni ↓ (FAAd)**"), Hi_d = md("**Hi ↓ (FRAd)**")
  ) %>%
  tab_spanner(label = md("**Frecuencia Relativa**"), columns = c(hi_pct, hi_real)) %>%
  tab_spanner(label = md("**Acumulada ↑**"),         columns = c(Ni_a, Hi_a)) %>%
  tab_spanner(label = md("**Acumulada ↓**"),         columns = c(Ni_d, Hi_d)) %>%
  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_column_spanners()) %>%
  tab_style(style = cell_fill(color = "#F5F5F5"),
            locations = cells_body(rows = seq(1, nrow(bind_rows(tabla_df, total_row)), by = 2))) %>%
  tab_style(style = list(cell_fill(color = "#D6D6D6"), cell_text(weight = "bold")),
            locations = cells_body(rows = Intervalo == "TOTAL", columns = everything())) %>%
  fmt_missing(columns = everything(), missing_text = "—") %>%
  tab_source_note(source_note = md("*Autor: Leslye Quinchiguango*")) %>%
  tab_options(table.width = pct(100), heading.title.font.size = px(16),
              heading.subtitle.font.size = px(12), table.font.size = px(13),
              data_row.padding = px(6))
Tabla N°1: Distribución de Frecuencias
Variable Cuantitativa Discreta Agrupada: Range, arrendamientos de hidrocarburos, Kansas, EE.UU. (n = 47,757 registros válidos)
Intervalo [Li – Ls) Marca de Clase ni (FA)
Frecuencia Relativa
Acumulada ↑
Acumulada ↓
hi % hi (decimal) Ni ↑ (FAAa) Hi ↑ (FRAa) Ni ↓ (FAAd) Hi ↓ (FRAd)
[1 – 6) 3 3318 6.95 0.0695 3318 0.0695 47757 1.0000
[6 – 11) 8 5137 10.76 0.1076 8455 0.1770 44439 0.9305
[11 – 16) 13 10092 21.13 0.2113 18547 0.3884 39302 0.8230
[16 – 21) 18 10946 22.92 0.2292 29493 0.6176 29210 0.6116
[21 – 26) 23 5443 11.40 0.1140 34936 0.7315 18264 0.3824
[26 – 31) 28 2478 5.19 0.0519 37414 0.7834 12821 0.2685
[31 – 36) 33 5874 12.30 0.1230 43288 0.9064 10343 0.2166
[36 – 41) 38 3279 6.87 0.0687 46567 0.9751 4469 0.0936
[41 – 43] 42 1190 2.49 0.0249 47757 1.0000 1190 0.0249
TOTAL — 47757 100.00 1.0000 47757 1.0000 47757 1.0000
Autor: Leslye Quinchiguango

5. Representación Gráfica Original

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

par(mar = c(5, 6, 6, 2))
plot(h_dens, col = grises, border = "black", freq = FALSE,
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = breaks_vec, labels = breaks_vec, las = 1, cex.axis = 0.9)
mtext("Densidad (%)",  side = 2, line = 4.5, cex = 1)
mtext("Range",         side = 1, line = 3.5, cex = 1)
mtext("Gráfica N°1: Histograma de Densidad (hi %) de la Variable Range,\narrendamientos de hidrocarburos, Kansas, EE.UU.",
      side = 3, line = 3, cex = 0.9, font = 2)

6. Conjetura del Modelo

El histograma de densidad de la variable Range revela un comportamiento diferenciado según la zona del territorio:

  • Zona 1 — Range [1, 31): las barras describen una forma ascendente con pico central, lo que sugiere una distribución Normal. Los arrendamientos se concentran progresivamente hacia los rangos medios y luego decrecen de forma simétrica.

  • Zona 2 — Range [31, 44]: a partir del intervalo 31 se observa una caída pronunciada y continua de las frecuencias, característica de una distribución Exponencial, donde la probabilidad decrece de manera acelerada hacia los valores más altos.

Por lo anterior, se trabaja con dos zonas con modelos distintos. El tramo intermedio entre ambas zonas no se excluye: el corte en 31 es el punto donde cambia el comportamiento distribucional.

corte_z <- 31

x_z1 <- x[x >= x_min  & x <  corte_z]
x_z2 <- x[x >= corte_z & x <= x_max]

cat("Zona 1 [", x_min, ",", corte_z, ") — n:", length(x_z1), "\n")
## Zona 1 [ 1 , 31 ) — n: 37414
cat("Zona 2 [", corte_z, ",", x_max,  "] — n:", length(x_z2), "\n")
## Zona 2 [ 31 , 43 ] — n: 10343
# ── Zona 1: Normal ────────────────────────────────────────────────────────────
fit_z1 <- fitdistr(x_z1, "normal")
mu_z1  <- fit_z1$estimate["mean"]
sd_z1  <- fit_z1$estimate["sd"]

# ── Zona 2: Exponencial (desplazada al mínimo de la zona) ────────────────────
offset_z2 <- min(x_z2) - 0.001
fit_z2    <- fitdistr(x_z2 - offset_z2, "exponential")
rate_z2   <- fit_z2$estimate["rate"]

cat("\n--- Zona 1: Normal ---\n")
## 
## --- Zona 1: Normal ---
cat("μ̂:", round(mu_z1, 4), "| σ̂:", round(sd_z1, 4), "\n")
## μ̂: 15.2986 | σ̂: 6.5852
cat("\n--- Zona 2: Exponencial ---\n")
## 
## --- Zona 2: Exponencial ---
cat("rate (λ̂):", round(rate_z2, 6), "| offset:", round(offset_z2, 4), "\n")
## rate (λ̂): 0.214281 | offset: 30.999

6.1 Gráfica N°2 — Histograma General con Corte por Zona

h_plot        <- hist(x, breaks = breaks_vec, plot = FALSE)
h_plot$density <- hi_dec * 100

# Colorear barras según zona
colores_zona <- ifelse(breaks_vec[-length(breaks_vec)] < corte_z, "gray35", "gray70")

par(mar = c(5, 6, 6, 2))
plot(h_plot, col = colores_zona, border = "black", freq = FALSE,
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = breaks_vec, labels = breaks_vec, las = 1, cex.axis = 0.9)
abline(v = corte_z, col = "black", lty = 2, lwd = 2)

yt <- max(hi_dec * 100)
text(mean(c(x_min, corte_z)), yt * 0.88,
     paste0("Zona 1\n[", x_min, " – ", corte_z, ")\nNormal"),
     cex = 0.85, font = 2)
text(mean(c(corte_z, x_max)), yt * 0.88,
     paste0("Zona 2\n[", corte_z, " – ", x_max, "]\nExponencial"),
     cex = 0.85, font = 2)

mtext("Densidad (%)", side = 2, line = 4.5, cex = 1)
mtext("Range",        side = 1, line = 3.5, cex = 1)
mtext(paste0("Gráfica N°2: Corte por Zona en Range = ", corte_z,
             " — Zona 1: Normal | Zona 2: Exponencial"),
      side = 3, line = 3, cex = 0.9, font = 2)

6.2 Gráficas Individuales por Zona

par(mfrow = c(2, 1))

# ── Zona 1: Normal ────────────────────────────────────────────────────────────
breaks_z1      <- seq(x_min, corte_z, by = c_amp)
if (tail(breaks_z1, 1) < corte_z) breaks_z1 <- c(breaks_z1, corte_z)
h_z1           <- hist(x_z1, breaks = breaks_z1, plot = FALSE)
hi_z1          <- h_z1$counts / length(x_z1) * 100
h_z1$density   <- hi_z1
k_z1           <- length(hi_z1)
grises_z1      <- gray(seq(0.20, 0.65, length.out = k_z1))

xs_z1 <- seq(x_min - 1, corte_z + 1, length.out = 500)
ys_z1 <- dnorm(xs_z1, mean = mu_z1, sd = sd_z1) * c_amp * 100

par(mar = c(5, 6, 6, 2))
plot(h_z1, col = grises_z1, border = "black", freq = FALSE,
     ylim = c(0, max(c(hi_z1, ys_z1)) * 1.25),
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = breaks_z1, labels = breaks_z1, las = 1, cex.axis = 0.9)
lines(xs_z1, ys_z1, col = "black", lwd = 2.5)
mtext("Densidad (%)", side = 2, line = 4.5, cex = 1)
mtext("Range",        side = 1, line = 3.5, cex = 1)
mtext(paste0("Gráfica N°3: Zona 1 [", x_min, " – ", corte_z,
             ") — Normal(μ̂ = ", round(mu_z1, 2), ", σ̂ = ", round(sd_z1, 2), ")"),
      side = 3, line = 3, cex = 0.9, font = 2)
legend("topright",
       legend = c("Histograma (hi %)", paste0("Normal(μ̂=", round(mu_z1,2), ", σ̂=", round(sd_z1,2), ")")),
       fill = c("gray45", NA), border = c("black", NA),
       lty = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.85)

# ── Zona 2: Exponencial ───────────────────────────────────────────────────────
breaks_z2      <- seq(corte_z, x_max + 1, by = c_amp)
if (tail(breaks_z2, 1) <= x_max) breaks_z2 <- c(breaks_z2, x_max + 1)
h_z2           <- hist(x_z2, breaks = breaks_z2, plot = FALSE)
hi_z2          <- h_z2$counts / length(x_z2) * 100
h_z2$density   <- hi_z2
k_z2           <- length(hi_z2)
grises_z2      <- gray(seq(0.40, 0.85, length.out = k_z2))

xs_z2     <- seq(corte_z - 0.5, x_max + 1, length.out = 500)
xs_z2_pos <- xs_z2 - offset_z2
xs_z2_pos[xs_z2_pos <= 0] <- 1e-9
ys_z2 <- dexp(xs_z2_pos, rate = rate_z2) * c_amp * 100

par(mar = c(5, 6, 6, 2))
plot(h_z2, col = grises_z2, border = "black", freq = FALSE,
     ylim = c(0, max(c(hi_z2, ys_z2)) * 1.25),
     main = "", xlab = "", ylab = "", las = 1, xaxt = "n")
axis(1, at = breaks_z2, labels = breaks_z2, las = 1, cex.axis = 0.9)
lines(xs_z2, ys_z2, col = "black", lwd = 2.5)
mtext("Densidad (%)", side = 2, line = 4.5, cex = 1)
mtext("Range",        side = 1, line = 3.5, cex = 1)
mtext(paste0("Gráfica N°4: Zona 2 [", corte_z, " – ", x_max,
             "] — Exponencial(λ̂ = ", round(rate_z2, 4), ")"),
      side = 3, line = 3, cex = 0.9, font = 2)
legend("topright",
       legend = c("Histograma (hi %)", paste0("Exponencial(λ̂=", round(rate_z2, 4), ")")),
       fill = c("gray65", NA), border = c("black", NA),
       lty = c(NA, 1), lwd = c(NA, 2.5), bty = "n", cex = 0.85)

par(mfrow = c(1, 1))

— Hasta aquí el punto 6. Pendiente revisión y aprobación antes de continuar. —

Autor: Leslye Quinchiguango