Inferencia de Zr en basaltos: modelo log-normal
Unidad 2 · Estadística inferencial · Geoquímica de rocas con datos de EarthChem PetDB
Cristoffer Rumiguano, Ana Gabriela Fernández, Evelyn Bustamante, Jessica Haro, Jhonattan Ligña
Noviembre de 2026
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 762. Tabla de distribución de cantidad
Se toma Zr 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[["Zr (ppm)"]]
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
## 185 8 177h <- 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 Zr en basaltos (muestra aleatoria)")| Li | Ls | MC | ni | hi |
|---|---|---|---|---|
| 0 | 50 | 25 | 17 | 9.60 |
| 50 | 100 | 75 | 59 | 33.33 |
| 100 | 150 | 125 | 47 | 26.55 |
| 150 | 200 | 175 | 32 | 18.08 |
| 200 | 250 | 225 | 13 | 7.34 |
| 250 | 300 | 275 | 9 | 5.08 |
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 Zr en basaltos", xlab = "Zr (ppm)", ylab = "Cantidad (ni)")
El histograma tiene una cola larga hacia la derecha y no baja de cero: la forma típica de los elementos traza, que se describe con un modelo log-normal (el logaritmo de Zr es normal).
4. Conjetura del modelo
Por la forma del histograma se conjetura un modelo
log-normal. fitdistr() estima sus parámetros con
los datos (máxima verosimilitud).
## meanlog sdlog
## 4.65715793 0.55621338
## (0.04180756) (0.02956241)media_log <- ajuste$estimate[["meanlog"]]
desv_log <- ajuste$estimate[["sdlog"]]
techo <- max(h$density, dlnorm(seq(min(x), max(x), length.out = 400), media_log, desv_log)) * 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 log-normal de Zr", xlab = "Zr (ppm)", ylab = "Densidad de probabilidad")
curve(dlnorm(x, media_log, desv_log), add = TRUE, col = "#d95926", lwd = 3) # curva del modelo
legend("topright", c("datos", "modelo log-normal"), 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 <- plnorm(Ls, media_log, desv_log) - plnorm(Li, media_log, desv_log) # 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")| Li | Ls | Fo | Fe |
|---|---|---|---|
| 0 | 50 | 17 | 15.96 |
| 50 | 100 | 59 | 65.95 |
| 100 | 150 | 47 | 48.62 |
| 150 | 200 | 32 | 24.43 |
| 200 | 250 | 13 | 11.40 |
| 250 | 300 | 9 | 5.34 |
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 = "Zr", Pearson = round(pearson, 2), Chi_cuadrado = round(x2, 2),
Umbral = round(umbral, 2), Resultado = resultado),
caption = "Tabla N° 3: pruebas de bondad de ajuste")| Variable | Pearson | Chi_cuadrado | Umbral | Resultado |
|---|---|---|---|---|
| Zr | 98.79 | 5.93 | 11.07 | se acepta el modelo |
El chi-cuadrado (5.93) es menor que el umbral (11.07): se acepta el modelo. La correlación entre lo observado y lo esperado es de 98.8 %.
6. Estimación de probabilidades
¿Cuál es la probabilidad de que un basalto tenga más de 150 ppm de Zr (basaltos enriquecidos, como los de islas oceánicas)?
## [1] 26.2549xs <- seq(min(x), max(x), length.out = 400)
plot(xs, dlnorm(xs, media_log, desv_log), type = "l", lwd = 2, las = 1, main = "Gráfica N° 4: P(Zr > 150 ppm)",
xlab = "Zr (ppm)", ylab = "Densidad de probabilidad")
zona <- xs[xs >= 150] # parte de la curva que se pinta
polygon(c(zona[1], zona, zona[length(zona)]), c(0, dlnorm(zona, media_log, desv_log), 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 50 y 100 ppm de Zr (típico de dorsales)?
p2 <- plnorm(100, media_log, desv_log) - plnorm(50, media_log, desv_log)
c(probabilidad = p2 * 100, basaltos_esperados_de_50 = p2 * 50)## probabilidad basaltos_esperados_de_50
## 37.25877 18.629387. Intervalos de confianza
Intervalo de confianza para la media poblacional de Zr 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
## 122.027 138.015 154.004n_necesario <- ceiling((z * sigma / (E / 2))^2) # tamaño para un error la mitad de grande
n_necesario## [1] 741Comprobació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[["Zr (ppm)"]]
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
## 130.8504 1.0000plot(IC, c(1, 1, 1), type = "b", pch = c(91, 19, 93), cex = 2, yaxt = "n", ylab = "", xlab = "Zr (ppm)",
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 Zr en los basaltos se representa con un modelo log-normal de parámetros media logarítmica = 4.657 y desviación logarítmica = 0.556 (chi-cuadrado 5.93 < umbral 11.07; Pearson 98.8 %). El circonio se concentra en pocos basaltos enriquecidos y la mayoría tiene valores bajos. Con el 95 % de confianza, la media poblacional está entre 122.03 y 154 ppm, y la media real de la base (130.85) cae dentro del intervalo.