1. CARGA DE LIBRERÍAS Y DATOS

library(gt)
library(tidyr)
library(ggplot2)
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
ruta_dataset <- "dataset_geologico_limpio_80.....csv"

if (!file.exists(ruta_dataset)) {
  message("Seleccione manualmente el archivo CSV del análisis.")
  ruta_dataset <- file.choose()
}

datos <- read.csv(
  ruta_dataset,
  header = TRUE, sep = ",", dec = ".",
  stringsAsFactors = FALSE
)

2. SELECCIÓN DE LA VARIABLE

clasificacion <- trimws(toupper(datos$CLASSIFICATION))
clasificacion <- clasificacion[
  !is.na(clasificacion) & clasificacion != ""
]

# Éxito: presencia de SAND.
# No éxito: cualquier otra clasificación sedimentaria.
# Esta codificación no ordena categorías según sus frecuencias.
exito <- ifelse(clasificacion == "SAND", 1, 0)

n <- 10
numero_bloques <- floor(length(exito) / n)
exito_completo <- exito[seq_len(numero_bloques * n)]
conteo_exitos <- rowSums(matrix(exito_completo, ncol = n, byrow = TRUE))

cat("Registros válidos:", length(clasificacion), "\n")
## Registros válidos: 27784
cat("Muestras clasificadas como SAND:", sum(exito), "\n")
## Muestras clasificadas como SAND: 15999
cat("Muestras con otras clasificaciones:", sum(exito == 0), "\n")
## Muestras con otras clasificaciones: 11785
cat("Bloques de 10 muestras:", numero_bloques, "\n")
## Bloques de 10 muestras: 2778

La variable CLASSIFICATION identifica la clasificación sedimentaria de cada muestra. Para aplicar el modelo binomial se plantea la siguiente pregunta: ¿cuántas muestras clasificadas como SAND aparecen en cada bloque de 10 observaciones?

Cada muestra se analiza con respecto a ese evento y puede producir dos resultados:

  • Éxito estadístico: la muestra está clasificada como SAND.
  • No éxito estadístico: la muestra pertenece a cualquier otra clasificación, como CLAY, SILT, GRAVEL, SILTY CLAY, entre otras.

En este contexto, las palabras éxito y no éxito solo permiten distinguir entre la presencia y la ausencia de SAND.

3. TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

valores <- 0:n
ni <- as.numeric(table(factor(conteo_exitos, levels = valores)))

TDF <- data.frame(Exitos = valores, ni = 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)))
  )

tabla_frecuencias <- bind_rows(
  mutate(TDF, Exitos = as.character(Exitos)),
  data.frame(
    Exitos = "TOTAL", 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 del número de clasificaciones SAND por bloque"
  ) %>%
  fmt_number(columns = c(hi, Hi_asc, Hi_dsc), decimals = 2) %>%
  sub_missing(columns = everything(), missing_text = "")
Tabla N.° 1
Distribución del número de clasificaciones SAND por bloque
Exitos ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
0 0 0.00 0 2778 0.00 100.00
1 12 0.43 12 2778 0.43 100.00
2 49 1.76 61 2766 2.20 99.57
3 161 5.80 222 2717 7.99 97.80
4 366 13.17 588 2556 21.17 92.01
5 598 21.53 1186 2190 42.69 78.83
6 680 24.48 1866 1592 67.17 57.31
7 535 19.26 2401 912 86.43 32.83
8 283 10.19 2684 377 96.62 13.57
9 79 2.84 2763 94 99.46 3.38
10 15 0.54 2778 15 100.00 0.54
TOTAL 2778 100.00



4. GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIAS

barplot(
  TDF$ni,
  names.arg = TDF$Exitos,
  col = "gray75", border = "gray30", space = 0.15,
  main = "Gráfica N.° 1\nClasificaciones SAND por bloque de 10 muestras",
  xlab = "Número de éxitos",
  ylab = "Frecuencia absoluta"
)

5. CONJETURA DEL MODELO

Mediante la gráfica se observa cuántas muestras SAND aparecen en cada bloque de 10 observaciones. El modelo no intenta ordenar todas las clasificaciones; únicamente distingue si SAND está presente o ausente en cada muestra. Debido a que existe un número fijo de ensayos, dos resultados posibles respecto al evento y una probabilidad constante de presencia, se plantea como conjetura un modelo binomial.

6. CÁLCULO DE PARÁMETROS

p <- mean(exito_completo)
N <- length(conteo_exitos)

cat("Número de ensayos por bloque (n):", n, "\n")
## Número de ensayos por bloque (n): 10
cat("Probabilidad estimada de éxito (p):", round(p, 6), "\n")
## Probabilidad estimada de éxito (p): 0.575846
cat("Número de bloques:", N, "\n")
## Número de bloques: 2778

El parámetro \(n=10\) representa el número fijo de muestras revisadas en cada bloque. El parámetro \(p\) representa la probabilidad estimada de encontrar una muestra SAND; por su parte, \(1-p\) representa la probabilidad de encontrar cualquiera de las demás clasificaciones. Por ello, la variable analizada es el conteo de muestras SAND dentro de cada bloque y puede tomar valores puntuales entre 0 y 10.

7. REALIDAD Y MODELO

Fe <- N * dbinom(valores, size = n, prob = p)

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 = "Realidad observada y modelo binomial"
  ) %>%
  fmt_number(
    columns = c(Frecuencia_Esperada, Diferencia, Error_Porcentual),
    decimals = 4
  )
Tabla N.° 2
Realidad observada y modelo binomial
Exitos ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc Frecuencia_Esperada Diferencia Error_Porcentual
0 0 0.0000000 0 2778 0.0000000 100.0000000 0.5236 −0.5236 52.3563
1 12 0.4319654 12 2778 0.4319654 100.0000000 7.1081 4.8919 40.7661
2 49 1.7638589 61 2766 2.1958243 99.5680346 43.4257 5.5743 11.3762
3 161 5.7955364 222 2717 7.9913607 97.8041757 157.2164 3.7836 2.3501
4 366 13.1749460 588 2556 21.1663067 92.0086393 373.5240 −7.5240 2.0558
5 598 21.5262779 1186 2190 42.6925846 78.8336933 608.5307 −10.5307 1.7610
6 680 24.4780418 1866 1592 67.1706263 57.3074154 688.4682 −8.4682 1.2453
7 535 19.2584593 2401 912 86.4290857 32.8293737 534.1073 0.8927 0.1669
8 283 10.1871850 2684 377 96.6162707 13.5709143 271.9208 11.0792 3.9149
9 79 2.8437725 2763 94 99.4600432 3.3837293 82.0375 −3.0375 3.8450
10 15 0.5399568 2778 15 100.0000000 0.5399568 11.1377 3.8623 25.7487
 
grafico <- comparativa %>%
  transmute(
    Exitos,
    Realidad = ni,
    Modelo = Frecuencia_Esperada
  ) %>%
  pivot_longer(
    cols = c(Realidad, Modelo),
    names_to = "Distribucion",
    values_to = "Frecuencia"
  )

ggplot(grafico, aes(factor(Exitos), Frecuencia, fill = Distribucion)) +
  geom_col(position = "dodge") +
  scale_fill_manual(
    values = c("Realidad" = "darkred", "Modelo" = "darkblue"),
    breaks = c("Realidad", "Modelo"),
    labels = c("Realidad observada", "Modelo binomial")
  ) +
  labs(
    title = "Gráfica N.° 2\nComparación entre realidad y modelo",
    subtitle = paste("n =", n, "; p =", round(p, 6)),
    x = "Número de éxitos",
    y = "Frecuencia absoluta",
    fill = "Distribución"
  ) +
  theme_classic() +
  theme(legend.position = "bottom")

8. TESTS DE APROBACIÓN

Fo_rel <- comparativa$ni / sum(comparativa$ni)
Fe_rel <- comparativa$Frecuencia_Esperada /
  sum(comparativa$Frecuencia_Esperada)
coef_pearson <- cor(Fo_rel, Fe_rel)

Fo_chi <- c(); Fe_chi <- c(); fo_temp <- 0; fe_temp <- 0
for (i in seq_along(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 <- pchisq(Chi2, gl, lower.tail = FALSE)

decision_pearson <- ifelse(coef_pearson >= 0.70, "APRUEBA", "NO APRUEBA")
decision_chi <- ifelse(p_valor > 0.05, "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 binomial"
  ) %>%
  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 binomial
Prueba Resultado Criterio Decision
Correlación de Pearson r = 0.9998 (99.98 %) Aprueba si r >= 0.70 APRUEBA
Chi-cuadrado X² = 5.6498; p = 0.6864 Aprueba si p > 0.05 APRUEBA
Pearson evalúa asociación; chi-cuadrado evalúa bondad de ajuste.

9. CÁLCULO DE PROBABILIDADES

x <- round(n * p)
prob_puntual <- dbinom(x, size = n, prob = p)
prob_acumulada <- pbinom(x, size = n, prob = p)

data.frame(
  Tipo = c("Puntual", "Acumulada"),
  Evento = c(
    paste("Exactamente", x, "éxitos"),
    paste(x, "éxitos o menos")
  ),
  Probabilidad = c(prob_puntual, prob_acumulada)
) %>%
  gt() %>%
  fmt_number(columns = Probabilidad, decimals = 6)
Tipo Evento Probabilidad
Puntual Exactamente 6 éxitos 0.247829
Acumulada 6 éxitos o menos 0.676313

10. INTERVALOS DE CONFIANZA

media_binomial <- n * p
desv_binomial <- sqrt(n * p * (1 - p))
z <- c(1, 1.96, 2.576)
limites_inferiores <- pmax(0, media_binomial - z * desv_binomial)
limites_superiores <- pmin(n, media_binomial + z * desv_binomial)

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 = "Límites decimales del número de éxitos"
  )
Tabla de intervalos de confianza originales
Límites decimales del número de éxitos
Nivel Limite_Inferior Limite_Superior
68% 4.195618 7.321301
95% 2.695290 8.821629
99% 1.732579 9.784339

Debido a que el número de éxitos es una variable discreta, sus intervalos se interpretan mediante valores puntuales. Por ello, los límites decimales se aproximan a los valores puntuales más cercanos.

data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = round(limites_inferiores),
  Limite_Superior = round(limites_superiores)
) %>%
  gt() %>%
  fmt_number(columns = 2:3, decimals = 0) %>%
  tab_header(
    title = md("**Tabla de intervalos de confianza puntuales**"),
    subtitle = "Aproximación a valores puntuales"
  )
Tabla de intervalos de confianza puntuales
Aproximación a valores puntuales
Nivel Limite_Inferior Limite_Superior
68% 4 7
95% 3 9
99% 2 10

11. CONCLUSIÓN

El comportamiento del número de muestras clasificadas como SAND en bloques de 10 observaciones se explica mediante un modelo binomial con parámetros \(n=10\) y \(p=0.5758\). Podemos afirmar con un 95 % de confianza que la media aritmética real del número de muestras SAND por bloque se encuentra entre 5.6992 y 5.8177, con una desviación estándar teórica de \(\sigma=\sqrt{np(1-p)}=1.5628\) muestras.