Proyecto estrictamente académico. 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 y librerías

Se lee la base limpia (159 760 análisis de roca total × 75 variables). La primera vez se descarga del repositorio del proyecto (26 MB).

knitr::opts_chunk$set(
  echo = TRUE,          # muestra el código en el documento
  message = FALSE,      # sin mensajes de carga de paquetes
  warning = FALSE,      # sin advertencias
  fig.align = "center" # gráficos centrados
)
library(e1071)         # asimetría y curtosis (install.packages("e1071"))
if (!file.exists("data/earthchem_limpio.csv.gz")) {
  dir.create("data", showWarnings = FALSE)
  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)
dim(base)              # filas y columnas
## [1] 159760     75

2. Datos

Se toma la columna de la variable y se quitan las filas vacías (vacío = no medido, nunca cero).

x <- base[["Cr (ppm)"]]
tipo <- base[["Sample Type"]]                        # tipo de roca de cada fila
tiene <- !is.na(x)                                     # filas con dato
x <- x[tiene]
tipo <- tipo[tiene]
x <- log10(x)                                          # escala log10: 10, 100, 1000 ppm -> 1, 2, 3
n <- length(x)

3. Límites de outliers (regla intercuartílica)

Con el rango intercuartílico (IR = Q3 − Q1) se fijan los límites Q1 − 1,5·IR y Q3 + 1,5·IR. Los datos que quedan fuera no se eliminan: se cuentan en clases abiertas.

Q1 <- quantile(x, 0.25, names = FALSE)
Q3 <- quantile(x, 0.75, names = FALSE)
IR <- Q3 - Q1                                          # rango intercuartílico
lim_inf <- Q1 - 1.5 * IR                               # menos que esto = outlier bajo
lim_sup <- Q3 + 1.5 * IR                               # más que esto = outlier alto
tipico <- x[x >= lim_inf & x <= lim_sup]               # los intervalos se arman con estos datos
# Los outliers NO se borran: se cuentan en las clases abiertas (son otras rocas, no errores).

4. Intervalos (regla de Sturges)

El número de intervalos sale de la regla de Sturges, k = 1 + 3,3·log10(n); la amplitud se redondea a un número cómodo para que los límites sean fáciles de leer.

k <- round(1 + 3.3 * log10(n))                         # número de intervalos
a <- (max(tipico) - min(tipico)) / k                   # amplitud calculada
a <- 0.2                                               # se redondea hacia arriba a un número cómodo
desde <- floor(min(tipico) / a) * a                    # primer límite: múltiplo de a
k <- floor((max(tipico) - desde) / a) + 1              # intervalos que hacen falta para llegar al máximo
Li <- desde + (0:(k - 1)) * a                          # límites inferiores
Ls <- Li + a                                           # límites superiores
MC <- (Li + Ls) / 2                                    # marca de clase

5. Frecuencias

Tabla de distribución de frecuencias: límites, marca de clase, frecuencia absoluta y relativa, y acumuladas ascendentes y descendentes.

dentro <- x[x >= Li[1] & x <= Ls[k]]
ni <- as.vector(table(cut(dentro, breaks = c(Li, Ls[k]), right = FALSE, include.lowest = TRUE)))  # cuenta por intervalo [Li, Ls[
abajo <- sum(x < Li[1])                                # clase abierta "menor que" (outliers bajos)
arriba <- sum(x > Ls[k])                               # clase abierta "mayor o igual" (outliers altos)
hi <- round(ni / n * 100, 2)                           # frecuencia relativa (%)
todos <- c(abajo, ni, arriba)                          # con las clases abiertas, para acumular
cerradas <- 2:(k + 1)                                  # posiciones de los intervalos normales
Ni_asc <- cumsum(todos)[cerradas]                      # acumulada ascendente
Ni_dsc <- rev(cumsum(rev(todos)))[cerradas]            # acumulada descendente
Hi_asc <- round(cumsum(todos) / n * 100, 2)[cerradas]
Hi_dsc <- round(rev(cumsum(rev(todos))) / n * 100, 2)[cerradas]
TDF <- data.frame(Li, Ls, MC, ni, hi, Ni_asc, Ni_dsc, Hi_asc, Hi_dsc)
knitr::kable(TDF, caption = "Tabla N° 1: distribución de frecuencias")
Tabla N° 1: distribución de frecuencias
Li Ls MC ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
0.8 1.0 0.9 649 1.71 2629 35931 6.93 94.78
1.0 1.2 1.1 968 2.55 3597 35282 9.49 93.07
1.2 1.4 1.3 1149 3.03 4746 34314 12.52 90.51
1.4 1.6 1.5 1196 3.15 5942 33165 15.67 87.48
1.6 1.8 1.7 1632 4.30 7574 31969 19.98 84.33
1.8 2.0 1.9 2203 5.81 9777 30337 25.79 80.02
2.0 2.2 2.1 3493 9.21 13270 28134 35.00 74.21
2.2 2.4 2.3 6832 18.02 20102 24641 53.02 65.00
2.4 2.6 2.5 7523 19.84 27625 17809 72.87 46.98
2.6 2.8 2.7 2750 7.25 30375 10286 80.12 27.13
2.8 3.0 2.9 1827 4.82 32202 7536 84.94 19.88
3.0 3.2 3.1 978 2.58 33180 5709 87.52 15.06
3.2 3.4 3.3 2048 5.40 35228 4731 92.92 12.48
3.4 3.6 3.5 2191 5.78 37419 2683 98.70 7.08
3.6 3.8 3.7 334 0.88 37753 492 99.58 1.30 
c(abajo = abajo, arriba = arriba)                      # datos en las clases abiertas (no se borra ninguno)
##  abajo arriba 
##   1980    158

6. Gráficos: histograma + polígono, ojiva y caja

Histograma con polígono de frecuencias, ojivas ascendente y descendente (se cruzan en la mediana) y diagrama de caja.

etiq <- paste0("[", Li, "-", Ls, "[")                  # texto de cada intervalo
nota <- c(if (abajo > 0) paste0("< ", Li[1], ": ", abajo, " (", round(abajo / n * 100, 2), " %)"),
          if (arriba > 0) paste0("≥ ", Ls[k], ": ", arriba, " (", round(arriba / n * 100, 2), " %)"))
layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE), heights = c(1.15, 1))   # arriba 1 gráfico, abajo 2
par(mar = c(8, 5, 4, 1))
h <- barplot(hi, names.arg = etiq, space = 0, las = 2, col = "#3987e5", border = "white",   # histograma
             ylab = "Porcentaje (%)", cex.names = .75, ylim = c(0, max(hi) * 1.15),
             main = "Gráfica No. 13: Histograma y polígono de frecuencias de log10 Cr (ppm)")
lines(c(h[1] - 1, h, h[length(h)] + 1), c(0, hi, 0), col = "#d95926", lwd = 2, type = "b", pch = 16, cex = .6)  # polígono
mtext(if (length(nota) > 0) paste0("Clases abiertas (fuera del gráfico): ", paste(nota, collapse = " · ")) else "",
      side = 1, line = 6.5, cex = .7, col = "gray30")  # nota bajo el histograma
par(mar = c(5, 5, 4, 1))
plot(Ls, Hi_asc, type = "b", pch = 16, col = "#3987e5", ylim = c(0, 100),            # ojiva ascendente
     xlab = "log10 Cr (ppm)", ylab = "Porcentaje acumulado (%)",
     main = "Gráfica No. 14: Ojiva ascendente y descendente")
lines(Li, Hi_dsc, type = "b", pch = 17, col = "#d95926")                              # ojiva descendente
abline(h = 50, lty = 3, col = "gray50")                # línea del 50 %: donde se cruzan está la mediana
legend("right", c("ascendente (≤)", "descendente (≥)"), col = c("#3987e5", "#d95926"), pch = c(16, 17), bty = "n", cex = .85)
boxplot(x, horizontal = TRUE, col = "#86b6ef", xlab = "log10 Cr (ppm)",       # diagrama de caja
        main = "Gráfica No. 15: Diagrama de caja (todos los datos)")

7. Error de agrupación: media y desviación con la TDF vs con los datos

Se comparan la media y la desviación calculadas con la TDF (cada dato = la marca de clase de su intervalo) con las calculadas con los datos; menos de 5 % de error indica que la agrupación representa bien los datos.

xc <- x[x >= Li[1] & x <= Ls[k]]                       # los datos de los intervalos normales
media_agrupada <- sum(MC * ni) / sum(ni)               # cada dato = la marca de clase de su intervalo
desv_agrupada <- sqrt(sum((MC - media_agrupada)^2 * ni) / sum(ni))
c(error_media = abs(media_agrupada - mean(xc)) / mean(xc) * 100,   # en %; menos de 5 % = bien
  error_desviacion = abs(desv_agrupada - sd(xc)) / sd(xc) * 100)
##      error_media error_desviacion 
##       0.02194378       0.27576692

8. Indicadores (con los datos sin agrupar, en log10)

Medidas de tendencia central, dispersión y forma con los datos sin agrupar (asimetría y curtosis con el paquete e1071).

ind <- c(n = n, media = mean(x), mediana = median(x), desviacion = sd(x), varianza = var(x),
  CV = sd(x) / mean(x) * 100, asimetria = skewness(x), curtosis = kurtosis(x))

knitr::kable(data.frame(indicador = names(ind), valor = round(ind, 4)), row.names = FALSE,
             caption = "Tabla N° 2: indicadores estadísticos")
Tabla N° 2: indicadores estadísticos
indicador valor
n 37911.0000
media 2.2816
mediana 2.3711
desviacion 0.7647
varianza 0.5848
CV 33.5166
asimetria -0.8409
curtosis 2.1577

9. Corrección de la agrupación

Con Sturges sobre todo el rango, los valores extremos estiran los intervalos y casi todos los datos caen en uno o dos. Por eso la TDF de arriba usa el rango típico (sin outliers) y clases abiertas.

9.1. TDF “por defecto” (como la plantilla): Sturges sobre TODO el rango, en la escala original

Para justificar la corrección, se arma la TDF como en la plantilla: Sturges sobre todo el rango, sin tratar los extremos.

x0 <- base[["Cr (ppm)"]]
x0 <- x0[!is.na(x0)]
k0 <- round(1 + 3.3 * log10(length(x0)))
a0 <- (max(x0) - min(x0)) / k0                         # amplitud sin redondear
Li0 <- min(x0) + (0:(k0 - 1)) * a0
Ls0 <- Li0 + a0
ni0 <- as.vector(table(cut(x0[x0 <= Ls0[k0]], breaks = c(Li0, Ls0[k0]), right = FALSE, include.lowest = TRUE)))
arriba0 <- sum(x0 > Ls0[k0])                           # por redondeo, el máximo puede quedar afuera
hi0 <- round(c(ni0, arriba0[arriba0 > 0]) / length(x0) * 100, 2)
num <- function(v) trimws(formatC(round(v, 1), format = "fg", digits = 6))   # número a texto: 1 decimal, 6 cifras
etiq0 <- c(paste0("[", num(Li0), "-", num(Ls0), "["), if (arriba0 > 0) paste0("≥ ", num(Ls0[k0])))

9.2. Gráfico: por defecto (izquierda) vs corregida (derecha)

A la izquierda, la agrupación por defecto (casi todo cae en uno o dos intervalos); a la derecha, la corregida.

etiq <- paste0("[", Li, "-", Ls, "[")
nota <- c(if (abajo > 0) paste0("< ", Li[1], ": ", abajo, " (", round(abajo / n * 100, 2), " %)"),
          if (arriba > 0) paste0("≥ ", Ls[k], ": ", arriba, " (", round(arriba / n * 100, 2), " %)"))
par(mfrow = c(1, 2), mar = c(8, 5, 4, 1))              # 2 histogramas lado a lado
barplot(hi0, names.arg = etiq0, space = 0, las = 2, cex.names = .65, col = "gray60", border = "white",
        ylab = "Porcentaje (%)", main = paste0("Por defecto: Sturges sobre el rango completo\nCr, k = ", length(hi0)))
barplot(hi, names.arg = etiq, space = 0, las = 2, cex.names = .65, col = "#3987e5", border = "white",
        ylab = "Porcentaje (%)", main = paste0("Corregida: rango típico (IQR) + clases abiertas\nCr (log10 ppm), k = ", k, " cerradas"))
mtext(if (length(nota) > 0) paste0("Clases abiertas (fuera del gráfico): ", paste(nota, collapse = " · ")) else "",
      side = 1, line = 6.5, cex = .65, col = "gray30")

10. Estrato: rocas volcánicas máficas

El mismo análisis para un solo tipo de roca (basaltos y afines). Al separar por estrato, la distribución deja de mezclar rocas distintas.

10.1. Datos

Se toma la columna de la variable y se quitan las filas vacías (vacío = no medido, nunca cero).

x <- base[["Cr (ppm)"]]
tipo <- base[["Sample Type"]]                        # tipo de roca de cada fila
tiene <- !is.na(x)                                     # filas con dato
x <- x[tiene]
tipo <- tipo[tiene]
x <- log10(x)                                          # escala log10: 10, 100, 1000 ppm -> 1, 2, 3
x <- x[tipo == "igneous:volcanic:mafic"]          # solo rocas volcánicas máficas
n <- length(x)

10.2. Límites de outliers (regla intercuartílica)

Con el rango intercuartílico (IR = Q3 − Q1) se fijan los límites Q1 − 1,5·IR y Q3 + 1,5·IR. Los datos que quedan fuera no se eliminan: se cuentan en clases abiertas.

Q1 <- quantile(x, 0.25, names = FALSE)
Q3 <- quantile(x, 0.75, names = FALSE)
IR <- Q3 - Q1                                          # rango intercuartílico
lim_inf <- Q1 - 1.5 * IR                               # menos que esto = outlier bajo
lim_sup <- Q3 + 1.5 * IR                               # más que esto = outlier alto
tipico <- x[x >= lim_inf & x <= lim_sup]               # los intervalos se arman con estos datos
# Los outliers NO se borran: se cuentan en las clases abiertas (son otras rocas, no errores).

10.3. Intervalos (regla de Sturges)

El número de intervalos sale de la regla de Sturges, k = 1 + 3,3·log10(n); la amplitud se redondea a un número cómodo para que los límites sean fáciles de leer.

k <- round(1 + 3.3 * log10(n))                         # número de intervalos
a <- (max(tipico) - min(tipico)) / k                   # amplitud calculada
a <- 0.2                                               # se redondea hacia arriba a un número cómodo
desde <- floor(min(tipico) / a) * a                    # primer límite: múltiplo de a
k <- floor((max(tipico) - desde) / a) + 1              # intervalos que hacen falta para llegar al máximo
Li <- desde + (0:(k - 1)) * a                          # límites inferiores
Ls <- Li + a                                           # límites superiores
MC <- (Li + Ls) / 2                                    # marca de clase

10.4. Frecuencias

Tabla de distribución de frecuencias: límites, marca de clase, frecuencia absoluta y relativa, y acumuladas ascendentes y descendentes.

dentro <- x[x >= Li[1] & x <= Ls[k]]
ni <- as.vector(table(cut(dentro, breaks = c(Li, Ls[k]), right = FALSE, include.lowest = TRUE)))  # cuenta por intervalo [Li, Ls[
abajo <- sum(x < Li[1])                                # clase abierta "menor que" (outliers bajos)
arriba <- sum(x > Ls[k])                               # clase abierta "mayor o igual" (outliers altos)
hi <- round(ni / n * 100, 2)                           # frecuencia relativa (%)
todos <- c(abajo, ni, arriba)                          # con las clases abiertas, para acumular
cerradas <- 2:(k + 1)                                  # posiciones de los intervalos normales
Ni_asc <- cumsum(todos)[cerradas]                      # acumulada ascendente
Ni_dsc <- rev(cumsum(rev(todos)))[cerradas]            # acumulada descendente
Hi_asc <- round(cumsum(todos) / n * 100, 2)[cerradas]
Hi_dsc <- round(rev(cumsum(rev(todos))) / n * 100, 2)[cerradas]
TDF <- data.frame(Li, Ls, MC, ni, hi, Ni_asc, Ni_dsc, Hi_asc, Hi_dsc)
knitr::kable(TDF, caption = "Tabla N° 3: distribución de frecuencias")
Tabla N° 3: distribución de frecuencias
Li Ls MC ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
1.4 1.6 1.5 548 2.26 2085 22764 8.58 93.68
1.6 1.8 1.7 940 3.87 3025 22216 12.45 91.42
1.8 2.0 1.9 1608 6.62 4633 21276 19.07 87.55
2.0 2.2 2.1 2772 11.41 7405 19668 30.47 80.93
2.2 2.4 2.3 5957 24.51 13362 16896 54.99 69.53
2.4 2.6 2.5 6825 28.09 20187 10939 83.07 45.01
2.6 2.8 2.7 2197 9.04 22384 4114 92.11 16.93
2.8 3.0 2.9 1279 5.26 23663 1917 97.37 7.89
3.0 3.2 3.1 390 1.60 24053 638 98.98 2.63 
c(abajo = abajo, arriba = arriba)                      # datos en las clases abiertas (no se borra ninguno)
##  abajo arriba 
##   1537    248

10.5. Gráficos: histograma + polígono, ojiva y caja

Histograma con polígono de frecuencias, ojivas ascendente y descendente (se cruzan en la mediana) y diagrama de caja.

etiq <- paste0("[", Li, "-", Ls, "[")                  # texto de cada intervalo
nota <- c(if (abajo > 0) paste0("< ", Li[1], ": ", abajo, " (", round(abajo / n * 100, 2), " %)"),
          if (arriba > 0) paste0("≥ ", Ls[k], ": ", arriba, " (", round(arriba / n * 100, 2), " %)"))
layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE), heights = c(1.15, 1))   # arriba 1 gráfico, abajo 2
par(mar = c(8, 5, 4, 1))
h <- barplot(hi, names.arg = etiq, space = 0, las = 2, col = "#3987e5", border = "white",   # histograma
             ylab = "Porcentaje (%)", cex.names = .75, ylim = c(0, max(hi) * 1.15),
             main = "Gráfica No. 16: Histograma y polígono de frecuencias de log10 Cr en rocas volcánicas máficas (ppm)")
lines(c(h[1] - 1, h, h[length(h)] + 1), c(0, hi, 0), col = "#d95926", lwd = 2, type = "b", pch = 16, cex = .6)  # polígono
mtext(if (length(nota) > 0) paste0("Clases abiertas (fuera del gráfico): ", paste(nota, collapse = " · ")) else "",
      side = 1, line = 6.5, cex = .7, col = "gray30")  # nota bajo el histograma
par(mar = c(5, 5, 4, 1))
plot(Ls, Hi_asc, type = "b", pch = 16, col = "#3987e5", ylim = c(0, 100),            # ojiva ascendente
     xlab = "log10 Cr en rocas volcánicas máficas (ppm)", ylab = "Porcentaje acumulado (%)",
     main = "Gráfica No. 17: Ojiva ascendente y descendente")
lines(Li, Hi_dsc, type = "b", pch = 17, col = "#d95926")                              # ojiva descendente
abline(h = 50, lty = 3, col = "gray50")                # línea del 50 %: donde se cruzan está la mediana
legend("right", c("ascendente (≤)", "descendente (≥)"), col = c("#3987e5", "#d95926"), pch = c(16, 17), bty = "n", cex = .85)
boxplot(x, horizontal = TRUE, col = "#86b6ef", xlab = "log10 Cr en rocas volcánicas máficas (ppm)",       # diagrama de caja
        main = "Gráfica No. 18: Diagrama de caja (todos los datos)")

10.6. Error de agrupación: media y desviación con la TDF vs con los datos

Se comparan la media y la desviación calculadas con la TDF (cada dato = la marca de clase de su intervalo) con las calculadas con los datos; menos de 5 % de error indica que la agrupación representa bien los datos.

xc <- x[x >= Li[1] & x <= Ls[k]]                       # los datos de los intervalos normales
media_agrupada <- sum(MC * ni) / sum(ni)               # cada dato = la marca de clase de su intervalo
desv_agrupada <- sqrt(sum((MC - media_agrupada)^2 * ni) / sum(ni))
c(error_media = abs(media_agrupada - mean(xc)) / mean(xc) * 100,   # en %; menos de 5 % = bien
  error_desviacion = abs(desv_agrupada - sd(xc)) / sd(xc) * 100)
##      error_media error_desviacion 
##       0.04005985       1.23913046

10.7. Indicadores (con los datos sin agrupar, en log10)

Medidas de tendencia central, dispersión y forma con los datos sin agrupar (asimetría y curtosis con el paquete e1071).

ind <- c(n = n, media = mean(x), mediana = median(x), desviacion = sd(x), varianza = var(x),
  CV = sd(x) / mean(x) * 100, asimetria = skewness(x), curtosis = kurtosis(x))

knitr::kable(data.frame(indicador = names(ind), valor = round(ind, 4)), row.names = FALSE,
             caption = "Tabla N° 4: indicadores estadísticos")
Tabla N° 4: indicadores estadísticos
indicador valor
n 24301.0000
media 2.2672
mediana 2.3655
desviacion 0.5030
varianza 0.2531
CV 22.1877
asimetria -1.8952
curtosis 7.4237

11. Conclusión

Se analizaron 37 911 casos con dato. La mediana es 235 ppm; en log10, la escala del análisis, la asimetría es -0,84 (cola hacia los valores bajos) y el coeficiente de variación es 33,52 %. 2138 datos quedaron en las clases abiertas: no se borraron, porque son otras rocas y no errores. El error de agrupación de la media es 0,02 %. El cromo varía en varios órdenes de magnitud (de menos de 1 ppm a cromititas con más de 10 000 ppm), por eso se analiza en log10. En esa escala la distribución es mucho más simétrica.