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

Se toma K2O de las rocas volcánicas máficas (basaltos y afines) de la muestra. Mezclar todos los tipos de roca da una distribución con varios picos que ningún modelo simple describe. Como en la plantilla del curso, los valores atípicos (regla de la caja) quedan fuera del ajuste: el modelo describe a los basaltos típicos. No se borran de la base y se informan aquí.

x <- muestra[["K2O (wt%)"]]
tipo <- muestra[["Sample Type"]]
x <- x[tipo == "igneous:volcanic:mafic" & !is.na(x) & x > 0]   # basaltos con dato
caja <- boxplot(x, plot = FALSE)                       # límites de la caja (bigotes)
atipicos <- x[x < caja$stats[1] | x > caja$stats[5]]
x <- x[x >= caja$stats[1] & x <= caja$stats[5]]        # datos para el modelo
c(basaltos_con_dato = length(x) + length(atipicos), atipicos = length(atipicos), para_el_modelo = length(x))
## basaltos_con_dato          atipicos    para_el_modelo 
##               197                17               180
h <- hist(x, plot = FALSE)                             # intervalos de Sturges (como hist() en la plantilla)
Li <- h$breaks[-length(h$breaks)]                      # límites inferiores
Ls <- h$breaks[-1]                                     # límites superiores
ni <- h$counts                                         # frecuencia absoluta
hi <- round(ni / sum(ni) * 100, 2)                     # frecuencia relativa (%)
TDF <- data.frame(Li, Ls, MC = h$mids, ni, hi)
knitr::kable(TDF, caption = "Tabla N° 1: distribución de K2O en basaltos (muestra aleatoria)")
Tabla N° 1: distribución de K2O en basaltos (muestra aleatoria)
Li Ls MC ni hi
0.0 0.2 0.1 54 30.00
0.2 0.4 0.3 42 23.33
0.4 0.6 0.5 20 11.11
0.6 0.8 0.7 25 13.89
0.8 1.0 0.9 14 7.78
1.0 1.2 1.1 12 6.67
1.2 1.4 1.3 6 3.33
1.4 1.6 1.5 3 1.67
1.6 1.8 1.7 1 0.56
1.8 2.0 1.9 2 1.11
2.0 2.2 2.1 1 0.56

3. Gráfica de distribución de cantidad

hist(x, breaks = h$breaks, col = "#86b6ef", border = "white", las = 1,
     main = "Gráfica N° 1: Distribución de K2O en basaltos", xlab = "K2O (wt%)", ylab = "Cantidad (ni)")

El histograma es máximo cerca de cero y decrece: la mayoría de los basaltos tiene muy poco potasio. Es la forma del modelo exponencial.

4. Conjetura del modelo

Por la forma del histograma se conjetura un modelo exponencial. fitdistr() estima sus parámetros con los datos (máxima verosimilitud).

ajuste <- fitdistr(x, "exponential")
ajuste                                                 # parámetros estimados y su error estándar
##      rate   
##   1.9815646 
##  (0.1476971)
lambda <- ajuste$estimate[["rate"]]   # tasa
techo <- max(h$density, dexp(seq(min(x), max(x), length.out = 400), lambda)) * 1.1                # alto del gráfico: que entren barras y curva
hist(x, breaks = h$breaks, freq = FALSE, col = "#86b6ef", border = "white", las = 1, ylim = c(0, techo),
     main = "Gráfica N° 2: Histograma y modelo exponencial de K2O", xlab = "K2O (wt%)", ylab = "Densidad de probabilidad")
curve(dexp(x, lambda), add = TRUE, col = "#d95926", lwd = 3)   # curva del modelo
legend("topright", c("datos", "modelo exponencial"), col = c("#86b6ef", "#d95926"), pch = c(15, NA),
       pt.cex = 2, lwd = c(NA, 3), bty = "n")

# Frecuencias observadas (Fo, los datos) y esperadas (Fe, lo que diría el modelo) en cada intervalo
Fo <- ni
P <- pexp(Ls, lambda) - pexp(Li, lambda)                                   # probabilidad de cada intervalo según el modelo
Fe <- P * sum(Fo)
knitr::kable(data.frame(Li, Ls, Fo, Fe = round(Fe, 2)), caption = "Tabla N° 2: frecuencias observadas y esperadas")
Tabla N° 2: frecuencias observadas y esperadas
Li Ls Fo Fe
0.0 0.2 54 58.90
0.2 0.4 42 39.63
0.4 0.6 20 26.66
0.6 0.8 25 17.94
0.8 1.0 14 12.07
1.0 1.2 12 8.12
1.2 1.4 6 5.46
1.4 1.6 3 3.68
1.6 1.8 1 2.47
1.8 2.0 2 1.66
2.0 2.2 1 1.12

5. Test de aprobación

Dos pruebas, como en la plantilla: la correlación de Pearson entre Fo y Fe (si el modelo es bueno, los puntos caen sobre la recta y = x) y el chi-cuadrado, que se compara con el umbral de la tabla al 95 % de confianza.

lim <- range(0, Fo / sum(Fo) * 100, Fe / sum(Fo) * 100)   # mismos ejes en x e y
plot(Fo / sum(Fo) * 100, Fe / sum(Fo) * 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)         # recta y = x: ajuste perfecto

pearson <- cor(Fo, Fe) * 100                           # en %
gl <- length(Fo) - 1                                   # grados de libertad
x2 <- sum((Fo - Fe)^2 / Fe)                            # chi-cuadrado calculado
umbral <- qchisq(0.95, gl)                             # valor crítico (95 %)
resultado <- ifelse(x2 < umbral, "se acepta el modelo", "se rechaza el modelo")
knitr::kable(data.frame(Variable = "K2O", 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
K2O 97.92 8.29 18.31 se acepta el modelo

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

6. Estimación de probabilidades

¿Cuál es la probabilidad de que un basalto tenga menos de 0,2 % de K2O (basalto de dorsal empobrecido)?

p1 <- pexp(0.2, lambda)
p1 * 100                                               # en %
## [1] 32.72039
xs <- seq(min(x), max(x), length.out = 400)
plot(xs, dexp(xs, lambda), type = "l", lwd = 2, las = 1, main = "Gráfica N° 4: P(K2O < 0,2 %)",
     xlab = "K2O (wt%)", ylab = "Densidad de probabilidad")
zona <- xs[xs <= 0.2]                                    # parte de la curva que se pinta
polygon(c(zona[1], zona, zona[length(zona)]), c(0, dexp(zona, lambda), 0), col = rgb(0.85, 0.35, 0.15, 0.5), border = NA)
legend("topright", paste0("P = ", round(p1 * 100, 2), " %"), bty = "n", cex = 1.1)

Si se analizan 50 basaltos nuevos, ¿cuántos tendrían entre 0,5 y 1 % de K2O?

p2 <- pexp(1, lambda) - pexp(0.5, lambda)
c(probabilidad = p2 * 100, basaltos_esperados_de_50 = p2 * 50)
##             probabilidad basaltos_esperados_de_50 
##                 23.34327                 11.67164

7. Intervalos de confianza

Intervalo de confianza para la media poblacional de K2O en basaltos, con todos los basaltos de la muestra (también los atípicos: aquí no se ajusta un modelo, se estima la media). Por el Teorema del Límite Central, la media de una muestra grande se distribuye normal con error estándar σ/√n; se toma σ = la desviación de la muestra (n grande).

todos <- c(x, atipicos)                                # todos los basaltos con dato
n <- length(todos)
media_m <- mean(todos)
sigma <- sd(todos)
z <- qnorm(0.975)                                      # 1,96 para el 95 %
E <- z * sigma / sqrt(n)                               # error máximo de estimación
IC <- c(inferior = media_m - E, media = media_m, superior = media_m + E)
round(IC, 3)
## inferior    media superior 
##    0.614    0.745    0.876
n_necesario <- ceiling((z * sigma / (E / 2))^2)        # tamaño para un error la mitad de grande
n_necesario
## [1] 788

Comprobación con la base completa. La población real se conoce: ¿el intervalo contiene la media de todos los basaltos de la base?

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)
poblacion <- base[["K2O (wt%)"]]
poblacion <- poblacion[base[["Sample Type"]] == "igneous:volcanic:mafic"]
poblacion <- poblacion[!is.na(poblacion) & poblacion > 0]
media_real <- mean(poblacion)
c(media_real = media_real, dentro_del_intervalo = media_real >= IC[["inferior"]] & media_real <= IC[["superior"]])
##           media_real dentro_del_intervalo 
##            0.6236786            1.0000000
plot(IC, c(1, 1, 1), type = "b", pch = c(91, 19, 93), cex = 2, yaxt = "n", ylab = "", xlab = "K2O (wt%)",
     xlim = range(c(IC, media_real)) + c(-E, E), main = "Gráfica N° 5: Intervalo del 95 % y media real")
abline(v = media_real, col = "#d95926", lwd = 2, lty = 2)
legend("topright", c("intervalo (muestra)", "media real (base)"), col = c("black", "#d95926"),
       lty = c(1, 2), lwd = 2, bty = "n")

8. Conclusión

El comportamiento de K2O en los basaltos se representa con un modelo exponencial de parámetros λ = 1.982 (media 1/λ = 0.5) (chi-cuadrado 8.29 < umbral 18.31; Pearson 97.9 %). El potasio separa a los basaltos de dorsal (pobres en K2O) de los de arco e intraplaca. Con el 95 % de confianza, la media poblacional está entre 0.61 y 0.88 % en peso, y la media real de la base (0.62) cae dentro del intervalo.