1. CARGA DE LIBRERÍAS Y DATOS

# Librerías utilizadas para tablas, transformación y gráficas.
library(gt)
library(tidyr)
library(ggplot2)
library(knitr)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Cargar el archivo de trabajo.
datos <- read.csv("dataset_geologico_limpio_80....csv",
                  header = TRUE, sep = ",", dec = ".",
                  stringsAsFactors = FALSE)

2. SELECCIÓN DE LA VARIABLE

# Extraer NMODES como variable cuantitativa discreta.
nmodos <- round(as.numeric(datos$NMODES))
nmodos <- na.omit(nmodos)
nmodos <- nmodos[nmodos >= 0]

cat("Número de observaciones:", length(nmodos), "\n")
## Número de observaciones: 27784
cat("Valor mínimo:", min(nmodos), "\n")
## Valor mínimo: 0
cat("Valor máximo:", max(nmodos), "\n")
## Valor máximo: 9

NMODES representa el número de modos granulométricos identificados en cada muestra geológica. Es una variable cuantitativa discreta porque corresponde a un conteo y toma valores enteros.

Los valores representan la cantidad de modos identificados: 0 indica que no se registró ningún modo; 1, una muestra unimodal; 2, una muestra bimodal; y 3 o más, una muestra multimodal.

3. TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

# Construir la distribución de frecuencias completa.
valores <- 0:max(nmodos)
categorias <- ifelse(valores == 0, "Sin modos",
              ifelse(valores == 1, "Unimodal",
              ifelse(valores == 2, "Bimodal",
                     paste(valores, "modos"))))
ni <- table(factor(nmodos, levels = valores))

TDF <- data.frame(Valor = valores, Categoria = categorias,
                  ni = as.numeric(ni)) %>%
  mutate(hi = ni / sum(ni) * 100,
         Ni_asc = cumsum(ni),
         Ni_dsc = rev(cumsum(rev(ni))),
         Hi_asc = cumsum(hi),
         Hi_dsc = rev(cumsum(rev(hi))))

# Incorporar la fila de comprobación de totales.
tabla_frecuencias <- TDF %>%
  mutate(Valor = as.character(Valor)) %>%
  rbind(data.frame(Valor = "TOTAL", Categoria = "",
                   ni = sum(TDF$ni), hi = 100,
                   Ni_asc = NA, Ni_dsc = NA, Hi_asc = NA, Hi_dsc = NA))

tabla_frecuencias %>% gt() %>%
  tab_header(title = md("**Tabla N.° 1**"),
             subtitle = "Distribución de frecuencias de NMODES") %>%
  fmt_number(columns = c(hi, Hi_asc, Hi_dsc), decimals = 6) %>%
  sub_missing(columns = everything(), missing_text = "")
Tabla N.° 1
Distribución de frecuencias de NMODES
Valor Categoria ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
0 Sin modos 6790 24.438526 6790 27784 24.438526 100.000000
1 Unimodal 9603 34.563058 16393 20994 59.001584 75.561474
2 Bimodal 6796 24.460121 23189 11391 83.461705 40.998416
3 3 modos 3041 10.945148 26230 4595 94.406853 16.538295
4 4 modos 1163 4.185862 27393 1554 98.592715 5.593147
5 5 modos 301 1.083357 27694 391 99.676073 1.407285
6 6 modos 65 0.233948 27759 90 99.910020 0.323927
7 7 modos 19 0.068385 27778 25 99.978405 0.089980
8 8 modos 3 0.010798 27781 6 99.989202 0.021595
9 9 modos 3 0.010798 27784 3 100.000000 0.010798
TOTAL 27784 100.000000



4. GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIAS

# Gráfica de la distribución observada.
# Ampliar el margen inferior para mostrar todos los nombres completos.
par(mar = c(10, 5, 4, 2))

barplot(TDF$ni, names.arg = TDF$Categoria, col = "gray75", border = "gray30",
        main = "Gráfica N.° 1\nDistribución del número de modos granulométricos",
        xlab = "", ylab = "Frecuencia absoluta",
        las = 2, cex.names = 0.80)

mtext("Clasificación modal", side = 1, line = 7)

# Restaurar los márgenes predeterminados.
par(mar = c(5.1, 4.1, 4.1, 2.1))

5. CONJETURA DEL MODELO

Mediante la gráfica se observa una mayor concentración en los primeros valores y una disminución de la frecuencia conforme aumenta el número de modos. Debido a que NMODES es una variable de conteo discreta, se plantea como conjetura que presenta una tendencia hacia un modelo de Poisson.

6. CÁLCULO DE PARÁMETROS

# Estimar lambda mediante la media muestral.
lambda <- mean(nmodos)
N <- length(nmodos)

cat("Parámetro lambda:", round(lambda, 6), "\n")
## Parámetro lambda: 1.405449
cat("Tamaño de la muestra:", N, "\n")
## Tamaño de la muestra: 27784

El parámetro \(\lambda\) representa el número promedio de modos granulométricos identificados en cada muestra geológica.

7. REALIDAD Y MODELO

# Calcular probabilidades y frecuencias esperadas.
P_esperada <- dpois(valores, lambda = lambda)
P_esperada[length(P_esperada)] <- P_esperada[length(P_esperada)] +
  (1 - ppois(max(valores), lambda = lambda))
Fe <- N * P_esperada

comparativa <- TDF %>%
  mutate(Frecuencia_Esperada = Fe,
         Diferencia = ni - Frecuencia_Esperada,
         Error_Porcentual = abs(Diferencia) / ifelse(ni == 0, 1, ni) * 100)

comparativa %>% gt() %>%
  tab_header(title = md("**Tabla N.° 2**"),
             subtitle = "Frecuencias observadas y esperadas del modelo de Poisson") %>%
  fmt_number(columns = c(Frecuencia_Esperada, Diferencia, Error_Porcentual),
             decimals = 6)
Tabla N.° 2
Frecuencias observadas y esperadas del modelo de Poisson
Valor Categoria ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc Frecuencia_Esperada Diferencia Error_Porcentual
0 Sin modos 6790 24.43852577 6790 27784 24.43853 100.00000000 6,814.216803 −24.216803 0.356654
1 Unimodal 9603 34.56305788 16393 20994 59.00158 75.56147423 9,577.035414 25.964586 0.270380
2 Bimodal 6796 24.46012093 23189 11391 83.46170 40.99841635 6,730.018282 65.981718 0.970890
3 3 modos 3041 10.94514829 26230 4595 94.40685 16.53829542 3,152.899557 −111.899557 3.679696
4 4 modos 1163 4.18586237 27393 1554 98.59272 5.59314714 1,107.810024 55.189976 4.745484
5 5 modos 301 1.08335733 27694 391 99.67607 1.40728477 311.394138 −10.394138 3.453202
6 6 modos 65 0.23394760 27759 90 99.91002 0.32392744 72.941439 −7.941439 12.217599
7 7 modos 19 0.06838468 27778 25 99.97840 0.08997984 14.645069 4.354931 22.920687
8 8 modos 3 0.01079758 27781 6 99.98920 0.02159516 2.572863 0.427137 14.237913
9 9 modos 3 0.01079758 27784 3 100.00000 0.01079758 0.466410 2.533590 84.452986
 
# Gráfica comparativa del modelo.
grafico <- comparativa %>% select(Categoria, ni, Frecuencia_Esperada) %>%
  pivot_longer(cols = c(ni, Frecuencia_Esperada),
               names_to = "Distribucion", values_to = "Frecuencia")

# Conservar el orden natural de las categorías: 0, 1, 2, 3, ... modos.
grafico$Categoria <- factor(
  grafico$Categoria,
  levels = categorias,
  ordered = TRUE
)

# Mostrar primero la frecuencia observada y luego la esperada.
grafico$Distribucion <- factor(
  grafico$Distribucion,
  levels = c("ni", "Frecuencia_Esperada")
)

ggplot(grafico, aes(Categoria, Frecuencia, fill = Distribucion)) +
  geom_col(position = "dodge") +
  scale_fill_manual(values = c("ni" = "darkred",
                               "Frecuencia_Esperada" = "darkblue"),
                    labels = c("Observada", "Esperada")) +
  labs(title = "Gráfica N.° 2\nComparación entre frecuencias observadas y esperadas",
       subtitle = paste("lambda =", round(lambda, 6)),
       x = "Clasificación modal", y = "Frecuencia absoluta",
       fill = "Distribución") +
  theme_classic() +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1, size = 10),
    axis.title.x = element_text(margin = margin(t = 12)),
    legend.position = "bottom",
    plot.margin = margin(10, 20, 20, 10)
  )

8. TESTS DE APROBACIÓN

# Correlación de Pearson entre frecuencias relativas.
Fo_rel <- comparativa$ni / sum(comparativa$ni)
Fe_rel <- comparativa$Frecuencia_Esperada /
  sum(comparativa$Frecuencia_Esperada)
coef_pearson <- cor(Fo_rel, Fe_rel, method = "pearson")

plot(Fo_rel, Fe_rel, pch = 19, col = "darkblue", cex = 1.25,
     xlab = "Frecuencia relativa observada",
     ylab = "Frecuencia relativa esperada",
     main = "Gráfica N.° 3\nCorrelación observada y esperada",
     xlim = c(0, max(Fo_rel) * 1.15),
     ylim = c(0, max(Fe_rel) * 1.15))
abline(0, 1, col = "red", lwd = 2)

cat("Coeficiente de Pearson:", round(coef_pearson, 6), "\n")
## Coeficiente de Pearson: 0.999911
cat("Porcentaje de ajuste:", round(coef_pearson * 100, 2), "%\n")
## Porcentaje de ajuste: 99.99 %
# Agrupar clases consecutivas cuando la frecuencia esperada sea menor que cinco.
Fo_chi <- c(); Fe_chi <- c(); fo_temp <- 0; fe_temp <- 0
for (i in 1:length(Fe)) {
  fo_temp <- fo_temp + comparativa$ni[i]
  fe_temp <- fe_temp + Fe[i]
  if (fe_temp >= 5) {
    Fo_chi <- c(Fo_chi, fo_temp); Fe_chi <- c(Fe_chi, fe_temp)
    fo_temp <- 0; fe_temp <- 0
  }
}
if (fe_temp > 0) {
  ultimo <- length(Fe_chi)
  Fo_chi[ultimo] <- Fo_chi[ultimo] + fo_temp
  Fe_chi[ultimo] <- Fe_chi[ultimo] + fe_temp
}

Chi2 <- sum((Fo_chi - Fe_chi)^2 / Fe_chi)
gl <- max(length(Fo_chi) - 2, 1)
p_valor <- 1 - pchisq(Chi2, gl)

cat("Estadístico chi-cuadrado:", round(Chi2, 6), "\n")
## Estadístico chi-cuadrado: 11.76219
cat("Grados de libertad:", gl, "\n")
## Grados de libertad: 6
cat("P-valor:", round(p_valor, 6), "\n")
## P-valor: 0.067489
if (p_valor > 0.05) {
  cat("DECISION: APRUEBA. No se rechaza H0; el ajuste es adecuado.\n")
} else {
  cat("DECISION: NO APRUEBA. Se rechaza H0; el ajuste no es adecuado.\n")
}
## DECISION: APRUEBA. No se rechaza H0; el ajuste es adecuado.

Resumen de los tests de aprobación

criterio_pearson <- 0.70
alpha <- 0.05

decision_pearson <- ifelse(
  coef_pearson >= criterio_pearson,
  "APRUEBA",
  "NO APRUEBA"
)

decision_chi <- ifelse(
  p_valor > alpha,
  "APRUEBA",
  "NO APRUEBA"
)

tabla_tests <- data.frame(
  Prueba = c("Correlación de Pearson", "Chi-cuadrado"),
  Resultado = c(
    paste0(
      "r = ", round(coef_pearson, 4),
      " (", round(coef_pearson * 100, 2), " %)"
    ),
    paste0(
      "X² = ", round(Chi2, 4),
      "; p = ", format.pval(p_valor, digits = 4)
    )
  ),
  Criterio = c(
    "Aprueba si r >= 0.70",
    "Aprueba si p > 0.05"
  ),
  Decision = c(decision_pearson, decision_chi)
)

tabla_tests %>%
  gt() %>%
  tab_header(
    title = md("**Tabla resumen de los tests de aprobación**"),
    subtitle = "Evaluación del ajuste al modelo de Poisson"
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = Decision)
  ) %>%
  tab_source_note(
    source_note = md(
      "__Pearson evalúa asociación; chi-cuadrado evalúa bondad de ajuste.__"
    )
  )
Tabla resumen de los tests de aprobación
Evaluación del ajuste al modelo de Poisson
Prueba Resultado Criterio Decision
Correlación de Pearson r = 0.9999 (99.99 %) Aprueba si r >= 0.70 APRUEBA
Chi-cuadrado X² = 11.7622; p = 0.06749 Aprueba si p > 0.05 APRUEBA
Pearson evalúa asociación; chi-cuadrado evalúa bondad de ajuste.

9. CÁLCULO DE PROBABILIDADES

9.1 Probabilidad puntual

# Seleccionar como referencia el número esperado de modos.
x <- round(lambda)
categoria_x <- ifelse(x == 0, "Sin modos",
               ifelse(x == 1, "Unimodal",
               ifelse(x == 2, "Bimodal", paste(x, "modos"))))

# Calcular la probabilidad de observar exactamente esa categoría.
prob_puntual <- dpois(x, lambda)
cat("P(X =", x, ") =", round(prob_puntual, 6), "\n")
## P(X = 1 ) = 0.344696
cat("Interpretación: probabilidad de obtener una muestra", categoria_x, "\n")
## Interpretación: probabilidad de obtener una muestra Unimodal

9.2 Probabilidad acumulada

# Calcular la probabilidad de observar como máximo x modos.
prob_acumulada <- ppois(x, lambda)
cat("P(X <=", x, ") =", round(prob_acumulada, 6), "\n")
## P(X <= 1 ) = 0.589953
cat("Interpretación: probabilidad de identificar", x,
    "modos o menos en una muestra.\n")
## Interpretación: probabilidad de identificar 1 modos o menos en una muestra.
tabla_probabilidades <- data.frame(
  Tipo = c("Puntual", "Acumulada"),
  Evento = c(paste("Muestra", categoria_x),
             paste(x, "modos o menos")),
  Probabilidad = c(prob_puntual, prob_acumulada)
)

tabla_probabilidades %>% gt() %>%
  fmt_number(columns = Probabilidad, decimals = 6)
Tipo Evento Probabilidad
Puntual Muestra Unimodal 0.344696
Acumulada 1 modos o menos 0.589953

10. INTERVALOS DE CONFIANZA

# Intervalos normales para la media del conteo.
error <- sd(nmodos) / sqrt(length(nmodos))

# Conservar primero los límites originales con decimales.
limites_inferiores <- mean(nmodos) -
  c(1, 1.96, 2.576) * error
limites_superiores <- mean(nmodos) +
  c(1, 1.96, 2.576) * error

tabla_ic_original <- data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = limites_inferiores,
  Limite_Superior = limites_superiores
)

tabla_ic_original %>%
  gt() %>%
  fmt_number(columns = 2:3, decimals = 6) %>%
  tab_header(
    title = md("**Tabla de intervalos de confianza originales**"),
    subtitle = "Estimación decimal de la media del número de modos"
  )
Tabla de intervalos de confianza originales
Estimación decimal de la media del número de modos
Nivel Limite_Inferior Limite_Superior
68% 1.398331 1.412567
95% 1.391497 1.419401
99% 1.387112 1.423786

Debido a que NMODES es una variable cuantitativa discreta, los límites decimales se interpretan mediante los valores puntuales posibles de la variable. Por esta razón, cada límite se aproxima al número puntual de modos más cercano. En este caso, las estimaciones de los tres niveles de confianza se concentran en el valor puntual de 1 modo.

# Aproximar los límites al valor puntual más cercano de la variable.
tabla_ic_puntual <- data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = round(limites_inferiores),
  Limite_Superior = round(limites_superiores)
)

tabla_ic_puntual %>%
  gt() %>%
  fmt_number(columns = 2:3, decimals = 0) %>%
  tab_header(
    title = md("**Tabla de intervalos de confianza puntuales**"),
    subtitle = "Aproximación al número puntual de modos más cercano"
  )
Tabla de intervalos de confianza puntuales
Aproximación al número puntual de modos más cercano
Nivel Limite_Inferior Limite_Superior
68% 1 1
95% 1 1
99% 1 1

11. CONCLUSIÓN

El comportamiento del número de modos granulométricos (NMODES) se explica mediante un modelo de Poisson con parámetro \(\lambda=1.4054\) modos por muestra. Podemos afirmar con un 95 % de confianza que la media aritmética real del número de modos granulométricos se encuentra entre 1.3915 y 1.4194 modos, con una desviación estándar teórica de \(\sigma=\sqrt{\lambda}=1.1855\) modos.