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

Variable ordinal: clase según SiO2 (ultrabásica < básica < intermedia < ácida). Se recorren los análisis en el orden del sorteo y se cuenta cuántos no son ultrabásicos antes de encontrar uno ultrabásico (el éxito). Ese número de fracasos es la variable aleatoria X.

x <- as.character(muestra[["SiO2 Class"]])
x <- x[!is.na(x) & x != ""]                            # quita vacíos (en el orden del sorteo)
exito <- x == "ultrabasica"                            # TRUE = rocas ultrabásicas
c(casos = length(x), exitos = sum(exito), proporcion = mean(exito))
##       casos      exitos  proporcion 
## 373.0000000  64.0000000   0.1715818
posicion <- which(exito)                               # dónde aparece cada ultrabásica
y <- diff(c(0, posicion)) - 1                          # fracasos antes de cada éxito
Fo <- as.vector(table(factor(y, levels = 0:max(y))))   # cuántas veces salió 0, 1, 2, ...
hi <- round(Fo / sum(Fo) * 100, 2)
TDF <- data.frame(valor = 0:max(y), ni = Fo, hi)
knitr::kable(TDF, caption = "Tabla N° 1: distribución de X = análisis antes de la primera roca ultrabásica")
Tabla N° 1: distribución de X = análisis antes de la primera roca ultrabásica
valor ni hi
0 14 21.88
1 8 12.50
2 9 14.06
3 5 7.81
4 5 7.81
5 4 6.25
6 2 3.12
7 3 4.69
8 2 3.12
9 2 3.12
10 4 6.25
11 0 0.00
12 1 1.56
13 0 0.00
14 1 1.56
15 2 3.12
16 0 0.00
17 0 0.00
18 0 0.00
19 0 0.00
20 0 0.00
21 0 0.00
22 0 0.00
23 1 1.56
24 0 0.00
25 0 0.00
26 0 0.00
27 0 0.00
28 0 0.00
29 0 0.00
30 0 0.00
31 0 0.00
32 1 1.56

3. Gráfica de distribución de cantidad

barplot(hi, names.arg = 0:max(y), col = "#86b6ef", border = "white", las = 1,
        main = "Gráfica N° 1: Análisis revisados antes de una roca ultrabásica", xlab = "Análisis antes de la primera ultrabásica", ylab = "Porcentaje (%)")

4. Conjetura del modelo

X cuenta los fracasos antes del primer éxito, en ensayos independientes con la misma probabilidad p: es el modelo geométrico. p se estima como 1 / (media + 1).

p <- 1 / (mean(y) + 1)                                 # probabilidad de éxito estimada
p
## [1] 0.1763085
valores <- 0:max(y)
P <- dgeom(valores, p)                                        # probabilidad de cada valor según el modelo
P[length(P)] <- 1 - pgeom(max(y) - 1, p)                                # el último valor acumula la cola ("o más")
Fe <- P * length(y)                                    # frecuencia esperada
# Para el chi-cuadrado cada clase debe esperar al menos 5 casos: los valores altos se juntan en una clase "o más"
corte <- max(which(Fe >= 5))
Fo2 <- c(Fo[1:(corte - 1)], sum(Fo[corte:length(Fo)]))
Fe2 <- c(Fe[1:(corte - 1)], sum(Fe[corte:length(Fe)]))
clase <- c(0:(corte - 2), paste0(corte - 1, " o más"))
knitr::kable(data.frame(clase, Fo = Fo2, Fe = round(Fe2, 2)), caption = "Tabla N° 2: frecuencias observadas y esperadas")
Tabla N° 2: frecuencias observadas y esperadas
clase Fo Fe
0 14 11.28
1 8 9.29
2 9 7.66
3 5 6.31
4 o más 28 29.46 
barplot(rbind(Fo2, Fe2), beside = TRUE, names.arg = clase, col = c("#86b6ef", "#d95926"), border = "white", las = 1,
        ylim = c(0, max(Fo2, Fe2) * 1.2),
        main = "Gráfica N° 2: Datos vs modelo geométrico", xlab = "Análisis antes de la primera ultrabásica", ylab = "Cantidad")
legend("topright", c("observado (Fo)", "modelo (Fe)"), fill = c("#86b6ef", "#d95926"), bty = "n")

5. Test de aprobación

lim <- range(0, Fo2 / sum(Fo2) * 100, Fe2 / sum(Fo2) * 100)   # mismos ejes en x e y
plot(Fo2 / sum(Fo2) * 100, Fe2 / sum(Fo2) * 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)

pearson <- cor(Fo2, Fe2) * 100
gl <- length(Fo2) - 1
x2 <- sum((Fo2 - Fe2)^2 / Fe2)
umbral <- qchisq(0.95, gl)
resultado <- ifelse(x2 < umbral, "se acepta el modelo", "se rechaza el modelo")
knitr::kable(data.frame(Variable = "X = análisis antes de la primera roca ultrabásica", 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
X = análisis antes de la primera roca ultrabásica 97.97 1.41 9.49 se acepta el modelo

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

6. Estimación de probabilidades

¿Cuál es la probabilidad de encontrar la primera roca ultrabásica antes de revisar 5 análisis (X ≤ 4)?

p1 <- pgeom(4, p)
p1 * 100
## [1] 62.08397
vp <- 0:15                                       # valores posibles que se dibujan
Pp <- dgeom(vp, p)
barplot(Pp * 100, names.arg = vp, las = 1, border = "white", col = ifelse(vp <= 4, "#d95926", "gray75"),
        main = "Gráfica N° 4: P(X ≤ 4)", xlab = "Análisis antes de la primera ultrabásica", ylab = "Probabilidad (%)")
legend("topright", paste0("P = ", round(p1 * 100, 2), " %"), bty = "n", cex = 1.1)

En promedio, ¿cuántos análisis hay que revisar para encontrar una roca ultrabásica?

c(fracasos_esperados = (1 - p) / p, analisis_a_revisar = 1 / p)
## fracasos_esperados analisis_a_revisar 
##           4.671875           5.671875

7. Intervalos de confianza

Intervalo de confianza del 95 % para la proporción poblacional de rocas ultrabásicas (la p del modelo), con la aproximación normal: p ± z·√(p(1 − p)/n).

p_m <- mean(exito)                                     # proporción en la muestra
z <- qnorm(0.975)
E <- z * sqrt(p_m * (1 - p_m) / length(exito))         # error máximo de estimación
IC <- c(inferior = p_m - E, proporcion = p_m, superior = p_m + E)
round(IC * 100, 2)                                     # en %
##   inferior proporcion   superior 
##      13.33      17.16      20.98

Comprobación con la base completa. La proporción real se conoce: ¿el intervalo la contiene?

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)
todos <- as.character(base[["SiO2 Class"]])
todos <- todos[!is.na(todos) & todos != ""]
p_real <- mean(todos == "ultrabasica")
c(proporcion_real = p_real, dentro_del_intervalo = p_real >= IC[["inferior"]] & p_real <= IC[["superior"]])
##      proporcion_real dentro_del_intervalo 
##            0.1646298            1.0000000

8. Conclusión

El número de análisis antes de encontrar una roca ultrabásica sigue un modelo geométrico con p = 0.176: en promedio hay que revisar 5.7 análisis. El modelo pasa el test (chi-cuadrado 1.41 < umbral 9.49; Pearson 98 %). Con el 95 % de confianza, la proporción poblacional de rocas ultrabásicas está entre 13.33 % y 20.98 %, y la proporción real de la base (16.46 %) cae dentro del intervalo.