Proyecto estrictamente académico (Estadística y Probabilidades, FIGEMPA – UCE). Datos: EarthChem PetDB (CC BY 4.0), descargados el 13/10/2026. Web del proyecto: https://crisrumiguano-art.github.io/geoquimica-earthchem-ml/

1. Carga de datos

La inferencia se hace con la muestra aleatoria simple de la Lección 3: 1 061 análisis sorteados de los 159 760 de la base (semilla 2026). Con la base entera, el test de chi-cuadrado rechaza cualquier modelo por diferencias mínimas; con una muestra aleatoria mide lo que importa. La base completa se usa al final para comprobar el intervalo de confianza.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center")
library(MASS)          # fitdistr(): estima los parámetros del modelo (viene con R)
dir.create("data/muestras", recursive = TRUE, showWarnings = FALSE)
if (!file.exists("data/muestras/muestra_aleatoria.csv")) {
  download.file("https://github.com/crisrumiguano-art/geoquimica-earthchem-ml/raw/main/data/muestras/muestra_aleatoria.csv", "data/muestras/muestra_aleatoria.csv", mode = "wb")
}
muestra <- read.csv("data/muestras/muestra_aleatoria.csv", check.names = FALSE)
dim(muestra)           # filas y columnas
## [1] 1061   76

2. Tabla de distribución de cantidad

Variable discreta: el número de rocas volcánicas intermedias (andesitas) en un lote de 30 análisis. Es un conteo de un evento poco frecuente (cerca del 6 % de los análisis). Antes se probó con N Analytes (número de analitos) y no sirvió: su varianza (122) es once veces su media (10,8) y el modelo de Poisson exige que sean iguales; el chi-cuadrado salió miles de veces mayor que el umbral.

x <- as.character(muestra[["Sample Type"]])
x <- x[!is.na(x) & x != ""]                            # quita vacíos (en el orden del sorteo)
exito <- x == "igneous:volcanic:intermediate"                            # TRUE = rocas volcánicas intermedias
c(casos = length(x), exitos = sum(exito), proporcion = mean(exito))
##        casos       exitos   proporcion 
## 1.061000e+03 6.000000e+01 5.655042e-02
lote <- 30
B <- floor(length(exito) / lote)
grupo <- rep(1:B, each = lote)
y <- as.vector(tapply(exito[1:(B * lote)], grupo, sum))  # andesitas en cada lote de 30
c(media = mean(y), varianza = var(y))                  # en Poisson deben ser parecidas
##    media varianza 
## 1.714286 1.680672
Fo <- as.vector(table(factor(y, levels = 0:max(y))))   # cuántas veces salió 0, 1, 2, ...
hi <- round(Fo / sum(Fo) * 100, 2)
TDF <- data.frame(valor = 0:max(y), ni = Fo, hi)
knitr::kable(TDF, caption = "Tabla N° 1: distribución de X = rocas intermedias en un lote de 30")
Tabla N° 1: distribución de X = rocas intermedias en un lote de 30
valor ni hi
0 5 14.29
1 13 37.14
2 9 25.71
3 5 14.29
4 1 2.86
5 2 5.71

3. Gráfica de distribución de cantidad

barplot(hi, names.arg = 0:max(y), col = "#86b6ef", border = "white", las = 1,
        main = "Gráfica N° 1: Rocas volcánicas intermedias en lotes de 30 análisis", xlab = "Rocas intermedias en el lote", ylab = "Porcentaje (%)")

4. Conjetura del modelo

X cuenta un evento raro en un lote grande y su media y varianza son parecidas: es el modelo de Poisson. λ se estima con la media.

lambda <- mean(y)                                     # promedio de andesitas por lote
lambda
## [1] 1.714286
valores <- 0:max(y)
P <- dpois(valores, lambda)                                        # probabilidad de cada valor según el modelo
P[length(P)] <- 1 - ppois(max(y) - 1, lambda)                                # el último valor acumula la cola ("o más")
Fe <- P * length(y)                                    # frecuencia esperada
# Para el chi-cuadrado cada clase debe esperar al menos 5 casos: los valores altos se juntan en una clase "o más"
corte <- max(which(Fe >= 5))
Fo2 <- c(Fo[1:(corte - 1)], sum(Fo[corte:length(Fo)]))
Fe2 <- c(Fe[1:(corte - 1)], sum(Fe[corte:length(Fe)]))
clase <- c(0:(corte - 2), paste0(corte - 1, " o más"))
knitr::kable(data.frame(clase, Fo = Fo2, Fe = round(Fe2, 2)), caption = "Tabla N° 2: frecuencias observadas y esperadas")
Tabla N° 2: frecuencias observadas y esperadas
clase Fo Fe
0 5 6.30
1 13 10.81
2 9 9.26
3 o más 8 8.63 
barplot(rbind(Fo2, Fe2), beside = TRUE, names.arg = clase, col = c("#86b6ef", "#d95926"), border = "white", las = 1,
        ylim = c(0, max(Fo2, Fe2) * 1.2),
        main = "Gráfica N° 2: Datos vs modelo Poisson", xlab = "Rocas intermedias en el lote", ylab = "Cantidad")
legend("topright", c("observado (Fo)", "modelo (Fe)"), fill = c("#86b6ef", "#d95926"), bty = "n")

5. Test de aprobación

lim <- range(0, Fo2 / sum(Fo2) * 100, Fe2 / sum(Fo2) * 100)   # mismos ejes en x e y
plot(Fo2 / sum(Fo2) * 100, Fe2 / sum(Fo2) * 100, pch = 19, col = "#3987e5", las = 1, xlim = lim, ylim = lim,
     main = "Gráfica N° 3: Frecuencia observada vs esperada", xlab = "Observada (%)", ylab = "Esperada (%)")
abline(a = 0, b = 1, col = "#d95926", lwd = 2)

pearson <- cor(Fo2, Fe2) * 100
gl <- length(Fo2) - 1
x2 <- sum((Fo2 - Fe2)^2 / Fe2)
umbral <- qchisq(0.95, gl)
resultado <- ifelse(x2 < umbral, "se acepta el modelo", "se rechaza el modelo")
knitr::kable(data.frame(Variable = "X = rocas intermedias en un lote de 30", Pearson = round(pearson, 2), Chi_cuadrado = round(x2, 2),
                        Umbral = round(umbral, 2), Resultado = resultado),
             caption = "Tabla N° 3: pruebas de bondad de ajuste")
Tabla N° 3: pruebas de bondad de ajuste
Variable Pearson Chi_cuadrado Umbral Resultado
X = rocas intermedias en un lote de 30 97.82 0.77 7.81 se acepta el modelo

El chi-cuadrado (0.77) es menor que el umbral (7.81): se acepta el modelo. La correlación entre lo observado y lo esperado es de 97.8 %.

6. Estimación de probabilidades

En un lote de 30 análisis, ¿cuál es la probabilidad de que no haya ninguna roca intermedia?

p1 <- dpois(0, lambda)
p1 * 100
## [1] 18.00923
vp <- 0:8                                       # valores posibles que se dibujan
Pp <- dpois(vp, lambda)
barplot(Pp * 100, names.arg = vp, las = 1, border = "white", col = ifelse(vp == 0, "#d95926", "gray75"),
        main = "Gráfica N° 4: P(X = 0) en un lote de 30", xlab = "Rocas intermedias en el lote", ylab = "Probabilidad (%)")
legend("topright", paste0("P = ", round(p1 * 100, 2), " %"), bty = "n", cex = 1.1)

En 50 lotes de 30 análisis, ¿en cuántos habría 3 o más rocas intermedias?

p3 <- 1 - ppois(2, lambda)
c(probabilidad = p3 * 100, lotes_esperados_de_50 = p3 * 50)
##          probabilidad lotes_esperados_de_50 
##              24.65526              12.32763

7. Intervalos de confianza

Intervalo de confianza del 95 % para la proporción poblacional de rocas intermedias (λ / 30), con la aproximación normal.

p_m <- mean(exito)                                     # proporción en la muestra
z <- qnorm(0.975)
E <- z * sqrt(p_m * (1 - p_m) / length(exito))         # error máximo de estimación
IC <- c(inferior = p_m - E, proporcion = p_m, superior = p_m + E)
round(IC * 100, 2)                                     # en %
##   inferior proporcion   superior 
##       4.27       5.66       7.04

Comprobación con la base completa. La proporción real se conoce: ¿el intervalo la contiene?

if (!file.exists("data/earthchem_limpio.csv.gz")) {
  download.file("https://github.com/crisrumiguano-art/geoquimica-earthchem-ml/raw/main/data/earthchem_limpio.csv.gz", "data/earthchem_limpio.csv.gz", mode = "wb")
}
base <- read.csv(gzfile("data/earthchem_limpio.csv.gz", encoding = "UTF-8"), check.names = FALSE)
todos <- as.character(base[["Sample Type"]])
todos <- todos[!is.na(todos) & todos != ""]
p_real <- mean(todos == "igneous:volcanic:intermediate")
c(proporcion_real = p_real, dentro_del_intervalo = p_real >= IC[["inferior"]] & p_real <= IC[["superior"]])
##      proporcion_real dentro_del_intervalo 
##           0.05831643           1.00000000

8. Conclusión

El número de rocas volcánicas intermedias en un lote de 30 análisis sigue un modelo de Poisson con λ = 1.71. El modelo pasa el test (chi-cuadrado 0.77 < umbral 7.81; Pearson 97.8 %). Con el 95 % de confianza, la proporción poblacional de rocas volcánicas intermedias está entre 4.27 % y 7.04 %, y la proporción real de la base (5.83 %) cae dentro del intervalo.