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 SiO2 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[["SiO2 (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 
##               189                23               166
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 SiO2 en basaltos (muestra aleatoria)")
Tabla N° 1: distribución de SiO2 en basaltos (muestra aleatoria)
Li Ls MC ni hi
45 46 45.5 5 3.01
46 47 46.5 5 3.01
47 48 47.5 14 8.43
48 49 48.5 27 16.27
49 50 49.5 37 22.29
50 51 50.5 41 24.70
51 52 51.5 22 13.25
52 53 52.5 8 4.82
53 54 53.5 4 2.41
54 55 54.5 3 1.81

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 SiO2 en basaltos", xlab = "SiO2 (wt%)", ylab = "Cantidad (ni)")

El histograma es una campana casi simétrica alrededor de 49–50 %: la forma de un modelo normal.

4. Conjetura del modelo

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

ajuste <- fitdistr(x, "normal")
ajuste                                                 # parámetros estimados y su error estándar
##       mean           sd     
##   49.79319963    1.77434853 
##  ( 0.13771615) ( 0.09738003)
media <- ajuste$estimate[["mean"]]
desv <- ajuste$estimate[["sd"]]
techo <- max(h$density, dnorm(seq(min(x), max(x), length.out = 400), media, desv)) * 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 normal de SiO2", xlab = "SiO2 (wt%)", ylab = "Densidad de probabilidad")
curve(dnorm(x, media, desv), add = TRUE, col = "#d95926", lwd = 3)   # curva del modelo
legend("topright", c("datos", "modelo 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 <- pnorm(Ls, media, desv) - pnorm(Li, media, desv)                                   # 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
45 46 5 2.13
46 47 5 6.88
47 48 14 16.33
48 49 27 28.44
49 50 37 36.35
50 51 41 34.10
51 52 22 23.47
52 53 8 11.86
53 54 4 4.40
54 55 3 1.20

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 = "SiO2", 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
SiO2 97.7 10.32 16.92 se acepta el modelo

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

6. Estimación de probabilidades

¿Cuál es la probabilidad de que un basalto tenga más de 51 % de SiO2 (cerca del límite de 52 % con las rocas intermedias)?

p1 <- 1 - pnorm(51, media, desv)
p1 * 100                                               # en %
## [1] 24.82088
xs <- seq(min(x), max(x), length.out = 400)
plot(xs, dnorm(xs, media, desv), type = "l", lwd = 2, las = 1, main = "Gráfica N° 4: P(SiO2 > 51 %)",
     xlab = "SiO2 (wt%)", ylab = "Densidad de probabilidad")
zona <- xs[xs >= 51]                                    # parte de la curva que se pinta
polygon(c(zona[1], zona, zona[length(zona)]), c(0, dnorm(zona, media, desv), 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 48 y 50 % de SiO2?

p2 <- pnorm(50, media, desv) - pnorm(48, media, desv)
c(probabilidad = p2 * 100, basaltos_esperados_de_50 = p2 * 50)
##             probabilidad basaltos_esperados_de_50 
##                 39.02935                 19.51468

7. Intervalos de confianza

Intervalo de confianza para la media poblacional de SiO2 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 
##   48.671   49.396   50.121
n_necesario <- ceiling((z * sigma / (E / 2))^2)        # tamaño para un error la mitad de grande
n_necesario
## [1] 756

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[["SiO2 (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 
##             49.25696              1.00000
plot(IC, c(1, 1, 1), type = "b", pch = c(91, 19, 93), cex = 2, yaxt = "n", ylab = "", xlab = "SiO2 (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")

Teorema del Límite Central. Si se sacan muchas muestras de 30 basaltos de la base y se calcula la media de cada una, esas medias forman una campana normal centrada en la media real, con desviación σ/√30, aunque los datos originales no sean normales.

set.seed(2026)
medias <- replicate(1000, mean(sample(poblacion, 30)))  # 1000 muestras de 30
hist(medias, breaks = 30, freq = FALSE, col = "#86b6ef", border = "white", las = 1,
     main = "Gráfica N° 6: Medias de 1000 muestras de 30 basaltos", xlab = "Media de la muestra", ylab = "Densidad")
curve(dnorm(x, mean(poblacion), sd(poblacion) / sqrt(30)), add = TRUE, col = "#d95926", lwd = 3)

c(desviacion_de_las_medias = sd(medias), sigma_sobre_raiz_n = sd(poblacion) / sqrt(30))
## desviacion_de_las_medias       sigma_sobre_raiz_n 
##                0.7832842                0.7829908

8. Conclusión

El comportamiento de SiO2 en los basaltos se representa con un modelo normal de parámetros media = 49.79 y desviación = 1.77 (chi-cuadrado 10.32 < umbral 16.92; Pearson 97.7 %). Los basaltos tienen una sílice muy constante: el SiO2 los define. Con el 95 % de confianza, la media poblacional está entre 48.67 y 50.12 % en peso, y la media real de la base (49.26) cae dentro del intervalo.