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 Mg# 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[["FeOT (wt%)"]]                             # hierro total como FeO
otra <- muestra[["FeO (wt%)"]] + 0.8998 * muestra[["Fe2O3 (wt%)"]]   # si falta, se calcula: FeOT = FeO + 0,8998·Fe2O3
x[is.na(x)] <- otra[is.na(x)]
otra <- 0.8998 * muestra[["Fe2O3T (wt%)"]]
x[is.na(x)] <- otra[is.na(x)]
mg <- muestra[["MgO (wt%)"]] / 40.30                       # moles de MgO (masa molar 40,30)
fe <- x / 71.85                                      # moles de FeO (masa molar 71,85)
x <- mg / (mg + fe)                                  # número de magnesio, entre 0 y 1
x[x >= 1] <- NA                                    # sin hierro medido no hay Mg# (la beta exige < 1)
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 
##               129                 8               121
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 Mg# en basaltos (muestra aleatoria)")
Tabla N° 1: distribución de Mg# en basaltos (muestra aleatoria)
Li Ls MC ni hi
0.25 0.30 0.275 1 0.83
0.30 0.35 0.325 2 1.65
0.35 0.40 0.375 10 8.26
0.40 0.45 0.425 14 11.57
0.45 0.50 0.475 12 9.92
0.50 0.55 0.525 23 19.01
0.55 0.60 0.575 25 20.66
0.60 0.65 0.625 23 19.01
0.65 0.70 0.675 6 4.96
0.70 0.75 0.725 2 1.65
0.75 0.80 0.775 3 2.48

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 Mg# en basaltos", xlab = "Mg# = Mg/(Mg+Fe)", ylab = "Cantidad (ni)")

El número de magnesio (Mg#) es la proporción de magnesio frente a magnesio más hierro: vale entre 0 y 1. Mide qué tan primitivo es un magma (los magmas recién salidos del manto tienen Mg# cercano a 0,7). Una variable acotada entre 0 y 1 se describe con el modelo beta, el único de la plantilla para proporciones.

4. Conjetura del modelo

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

ajuste <- fitdistr(x, "beta", start = list(shape1 = 2, shape2 = 2))
ajuste                                                 # parámetros estimados y su error estándar
##     shape1      shape2  
##   13.083151   11.301773 
##  ( 1.667895) ( 1.436445)
a <- ajuste$estimate[["shape1"]]
b <- ajuste$estimate[["shape2"]]
techo <- max(h$density, dbeta(seq(min(x), max(x), length.out = 400), a, b)) * 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 beta de Mg#", xlab = "Mg# = Mg/(Mg+Fe)", ylab = "Densidad de probabilidad")
curve(dbeta(x, a, b), add = TRUE, col = "#d95926", lwd = 3)   # curva del modelo
legend("topright", c("datos", "modelo beta"), 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 <- pbeta(Ls, a, b) - pbeta(Li, a, b)                                   # 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.25 0.30 1 0.79
0.30 0.35 2 2.76
0.35 0.40 10 6.89
0.40 0.45 14 13.07
0.45 0.50 12 19.49
0.50 0.55 23 23.22
0.55 0.60 25 22.18
0.60 0.65 23 16.81
0.65 0.70 6 9.86
0.70 0.75 2 4.27
0.75 0.80 3 1.26

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 = "Mg#", 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
Mg# 91.56 12.38 18.31 se acepta el modelo

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

6. Estimación de probabilidades

¿Cuál es la probabilidad de que un basalto tenga Mg# mayor que 0,65 (magma primitivo, poco evolucionado)?

p1 <- 1 - pbeta(0.65, a, b)
p1 * 100                                               # en %
## [1] 12.9147
xs <- seq(min(x), max(x), length.out = 400)
plot(xs, dbeta(xs, a, b), type = "l", lwd = 2, las = 1, main = "Gráfica N° 4: P(Mg# > 0,65)",
     xlab = "Mg# = Mg/(Mg+Fe)", ylab = "Densidad de probabilidad")
zona <- xs[xs >= 0.65]                                    # parte de la curva que se pinta
polygon(c(zona[1], zona, zona[length(zona)]), c(0, dbeta(zona, a, b), 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 Mg# entre 0,5 y 0,6?

p2 <- pbeta(0.6, a, b) - pbeta(0.5, a, b)
c(probabilidad = p2 * 100, basaltos_esperados_de_50 = p2 * 50)
##             probabilidad basaltos_esperados_de_50 
##                 37.51996                 18.75998

7. Intervalos de confianza

Intervalo de confianza para la media poblacional de Mg# 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.530    0.551    0.573
n_necesario <- ceiling((z * sigma / (E / 2))^2)        # tamaño para un error la mitad de grande
n_necesario
## [1] 516

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[["FeOT (wt%)"]]                             # hierro total como FeO
otra <- base[["FeO (wt%)"]] + 0.8998 * base[["Fe2O3 (wt%)"]]   # si falta, se calcula: FeOT = FeO + 0,8998·Fe2O3
poblacion[is.na(poblacion)] <- otra[is.na(poblacion)]
otra <- 0.8998 * base[["Fe2O3T (wt%)"]]
poblacion[is.na(poblacion)] <- otra[is.na(poblacion)]
mg <- base[["MgO (wt%)"]] / 40.30                       # moles de MgO (masa molar 40,30)
fe <- poblacion / 71.85                                      # moles de FeO (masa molar 71,85)
poblacion <- mg / (mg + fe)                                  # número de magnesio, entre 0 y 1
poblacion[poblacion >= 1] <- NA                                    # sin hierro medido no hay Mg# (la beta exige < 1)
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.5655675            1.0000000
plot(IC, c(1, 1, 1), type = "b", pch = c(91, 19, 93), cex = 2, yaxt = "n", ylab = "", xlab = "Mg# = Mg/(Mg+Fe)",
     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 Mg# en los basaltos se representa con un modelo beta de parámetros forma a = 13.08 y forma b = 11.3 (media a/(a+b) = 0.537) (chi-cuadrado 12.38 < umbral 18.31; Pearson 91.6 %). La mayoría de los basaltos ya perdió parte de su magnesio al cristalizar olivino: no son magmas primitivos. Con el 95 % de confianza, la media poblacional está entre 0.53 y 0.57 (proporción), y la media real de la base (0.57) cae dentro del intervalo.