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)

# 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

2. TABLA Y GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIAS

NMODES representa el número de modos granulométricos identificados en cada muestra geológica. Al ser un conteo de ocurrencias por muestra, puede evaluarse mediante un modelo de Poisson.

Los valores del eje horizontal 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.

# 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 



# 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))

3. CONJETURA DEL MODELO

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

# 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)
  )

4. 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("No se rechaza H0: el ajuste es adecuado.\n")
## No se rechaza H0: el ajuste es adecuado.

5. CÁLCULO DE PROBABILIDADES

5.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

5.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.

5.3 PROBABILIDAD COMPLEMENTARIA

# Calcular la probabilidad de observar más de x modos.
prob_complementaria <- 1 - ppois(x, lambda)
cat("P(X >", x, ") =", round(prob_complementaria, 6), "\n")
## P(X > 1 ) = 0.410047
cat("Interpretación: probabilidad de identificar más de", x,
    "modos en una muestra.\n")
## Interpretación: probabilidad de identificar más de 1 modos en una muestra.

5.4 TABLA RESUMEN DE PROBABILIDADES

tabla_probabilidades <- data.frame(
  Tipo = c("Puntual", "Acumulada", "Complementaria"),
  Evento = c(paste("Muestra", categoria_x),
             paste(x, "modos o menos"),
             paste("Más de", x, "modos")),
  Probabilidad = c(prob_puntual, prob_acumulada, prob_complementaria)
)
tabla_probabilidades %>% gt() %>%
  fmt_number(columns = Probabilidad, decimals = 6)
Tipo Evento Probabilidad
Puntual Muestra Unimodal 0.344696
Acumulada 1 modos o menos 0.589953
Complementaria Más de 1 modos 0.410047

6. INTERVALOS DE CONFIANZA

# Intervalos normales para la media del conteo.
error <- sd(nmodos) / sqrt(length(nmodos))
tabla_ic <- data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = mean(nmodos) - c(1, 1.96, 2.576) * error,
  Limite_Superior = mean(nmodos) + c(1, 1.96, 2.576) * error)
tabla_ic %>% gt() %>% fmt_number(columns = 2:3, decimals = 6)
Nivel Limite_Inferior Limite_Superior
68% 1.398331 1.412567
95% 1.391497 1.419401
99% 1.387112 1.423786

7. CONCLUSIÓN

La variable número de modos granulométricos (NMODES) se explica mediante un modelo de Poisson, con una media aritmética de 1.41 modos por muestra y una desviación estándar de 1.19. El parámetro lambda representa el número promedio de modos identificados en cada muestra geológica.

De esta manera se calcularon probabilidades como, por ejemplo, que al seleccionar aleatoriamente una muestra geológica, la probabilidad de identificar exactamente 1 modo es de 34.47 %, mientras que la probabilidad de identificar 1 modo o menos es de 59 %.

Mediante el teorema del límite central, se estima que la media poblacional del número de modos granulométricos se encuentra entre 1.39 y 1.42 modos, con un 95 % de confianza. La prueba chi-cuadrado obtuvo un p-valor de 0.0675; al ser superior a 0.05, no se rechaza la hipótesis nula y el modelo de Poisson se considera adecuado.