library(ggplot2)
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(scales)
library(patchwork)Análisis de la Teoría del Muestreo (TOS) — Error Fundamental de Muestreo (FSE)
Caso: Mineral de Niobio — Mina Boa Vista, CMOC Brasil
Parámetros del Caso Boa Vista
# ── Constantes calibradas del test de heterogeneidad (Chieregati et al., 2025)
# Ecuación de regresión log-log: IHL = 11.493 * d^3.0976 (R² = 0.9952)
K_cal <- 11.493 # constante de muestreo calibrada [g/cm^alpha]
alpha <- 3.0976 # exponente empírico de dN
# ── Parámetros del protocolo actual (Tabla 1 del artículo)
protocolo_actual <- data.frame(
etapa = c("1. Muestreo primario y splitting",
"2. Trituración",
"3. División",
"4. Pulverización",
"5. Selección muestra analítica"),
ML_g = c(100e6, 6000, 6000, 500, 500),
MS_g = c(6000, 6000, 500, 500, 60),
dN_cm = c(1.70, 0.20, 0.20, 0.0104, 0.0104)
)
# ── Protocolo optimizado (Tabla 2 del artículo)
protocolo_opt <- protocolo_actual
protocolo_opt$MS_g[1] <- 90000 # 90 kg en etapa primaria
protocolo_opt$ML_g[c(2,3)] <- 90000Estimación del FSE
Test de Heterogeneidad — Calibración de K y α
# Datos del test de heterogeneidad (extraídos de Fig. 2 del artículo)
# Fracciones: -17+12.7 mm, -12.7+6.35 mm, -6.35+4.75 mm, -4.75+3.35 mm
d1 <- c(12.7, 6.35, 4.75, 3.35)
d2 <- c(17.0, 12.7, 6.35, 4.75)
# Diámetro nominal según Ec. 5: dN = ((d1³ + d2³)/2)^(1/3) [cm]
dN_cm <- ((( d1^3 + d2^3) / 2)^(1/3)) / 10
# IHL estimados aproximados a partir de la gráfica del artículo (Fig. 2)
IHL_est <- K_cal * dN_cm^alpha
ht_data <- data.frame(dN_cm = dN_cm, IHL = IHL_est,
fraccion = paste0("-", d2, "+", d1, " mm"))
# Curva de regresión suave
d_curve <- seq(0.03, 2.0, length.out = 200)
curva <- data.frame(dN_cm = d_curve,
IHL = K_cal * d_curve^alpha)
ggplot() +
geom_line(data = curva, aes(x = dN_cm, y = IHL),
color = "#c62828", linewidth = 1.0) +
geom_point(data = ht_data, aes(x = dN_cm, y = IHL),
shape = 21, fill = "#1565c0", color = "white",
size = 4, stroke = 1.5) +
geom_text(data = ht_data,
aes(x = dN_cm, y = IHL, label = fraccion),
hjust = -0.1, vjust = 0.5, size = 3, color = "gray30") +
annotate("text", x = 0.1, y = 40,
label = expression(IH[L] == "11.493" ~ d[N]^{3.0976} ~~ (R^2 == 0.9952)),
size = 3.8, color = "#c62828", hjust = 0) +
scale_x_log10(breaks = c(0.05, 0.1, 0.5, 1.0, 2.0),
labels = c("0.05","0.1","0.5","1.0","2.0")) +
scale_y_log10(breaks = c(0.1, 1, 10, 100),
labels = c("0.1","1","10","100")) +
annotation_logticks(sides = "bl", color = "gray60", linewidth = 0.3) +
labs(
title = expression(IH[L] ~ "vs. Tamaño nominal de fragmento — Test de Heterogeneidad"),
subtitle = "Mineral de niobio (pirocloro), Mina Boa Vista — CMOC Brasil",
x = expression(d[N] ~ "(cm)"),
y = expression(IH[L] ~ "(g)"),
caption = "Fuente: reproducido de Chieregati et al. (2025, Fig. 2)"
)
Cálculo del FSE por Etapa del Protocolo
# Función: varianza FSE usando constantes calibradas K y alpha
calc_FSE <- function(MS_g, ML_g, dN_cm, K, a) {
K * dN_cm^a * (1/MS_g - 1/ML_g)
}
# Calcular FSE para protocolo actual y optimizado
calc_protocolo <- function(df, K, a) {
df %>%
mutate(
IHL_g = K * dN_cm^a,
var_FSE = calc_FSE(MS_g, ML_g, dN_cm, K, a),
desv_FSE = sqrt(abs(var_FSE)) * 100
)
}
res_actual <- calc_protocolo(protocolo_actual, K_cal, alpha)
res_opt <- calc_protocolo(protocolo_opt, K_cal, alpha)
total_actual <- sqrt(sum(res_actual$var_FSE, na.rm = TRUE)) * 100
total_opt <- sqrt(sum(res_opt$var_FSE, na.rm = TRUE)) * 100res_actual %>%
select(etapa, ML_g, MS_g, dN_cm, IHL_g, var_FSE, desv_FSE) %>%
mutate(across(where(is.numeric), ~ round(.x, 6))) %>%
rename(
"Etapa" = etapa,
"ML (g)" = ML_g,
"MS (g)" = MS_g,
"dN (cm)" = dN_cm,
"IHL (g)" = IHL_g,
"s²(FSE)" = var_FSE,
"s(FSE) (%)" = desv_FSE
) %>%
kable(align = c("l","r","r","r","r","r","r")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE) %>%
row_spec(0, bold = TRUE, background = "#1F3864", color = "white") %>%
row_spec(nrow(res_actual), bold = TRUE) %>%
footnote(general = paste0("TOTAL s(FSE) = ", round(total_actual, 2), "%"),
general_title = "")| Etapa | ML (g) | MS (g) | dN (cm) | IHL (g) | s²(FSE) | s(FSE) (%) |
|---|---|---|---|---|---|---|
| 1. Muestreo primario y splitting | 1e+08 | 6000 | 1.7000 | 59.466446 | 0.009910 | 9.955139 |
| 2. Trituración | 6e+03 | 6000 | 0.2000 | 0.078579 | 0.000000 | 0.000000 |
| 3. División | 6e+03 | 500 | 0.2000 | 0.078579 | 0.000144 | 1.200253 |
| 4. Pulverización | 5e+02 | 500 | 0.0104 | 0.000008 | 0.000000 | 0.000000 |
| 5. Selección muestra analítica | 5e+02 | 60 | 0.0104 | 0.000008 | 0.000000 | 0.034847 |
| TOTAL s(FSE) = 10.03% |
res_opt %>%
select(etapa, ML_g, MS_g, dN_cm, IHL_g, var_FSE, desv_FSE) %>%
mutate(across(where(is.numeric), ~ round(.x, 6))) %>%
rename(
"Etapa" = etapa,
"ML (g)" = ML_g,
"MS (g)" = MS_g,
"dN (cm)" = dN_cm,
"IHL (g)" = IHL_g,
"s²(FSE)" = var_FSE,
"s(FSE) (%)" = desv_FSE
) %>%
kable(align = c("l","r","r","r","r","r","r")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE) %>%
row_spec(0, bold = TRUE, background = "#2e7d32", color = "white") %>%
row_spec(nrow(res_opt), bold = TRUE) %>%
footnote(general = paste0("TOTAL s(FSE) = ", round(total_opt, 2), "%"),
general_title = "")| Etapa | ML (g) | MS (g) | dN (cm) | IHL (g) | s²(FSE) | s(FSE) (%) |
|---|---|---|---|---|---|---|
| 1. Muestreo primario y splitting | 1e+08 | 90000 | 1.7000 | 59.466446 | 0.000660 | 2.569326 |
| 2. Trituración | 9e+04 | 6000 | 0.2000 | 0.078579 | 0.000012 | 0.349619 |
| 3. División | 9e+04 | 500 | 0.2000 | 0.078579 | 0.000156 | 1.250136 |
| 4. Pulverización | 5e+02 | 500 | 0.0104 | 0.000008 | 0.000000 | 0.000000 |
| 5. Selección muestra analítica | 5e+02 | 60 | 0.0104 | 0.000008 | 0.000000 | 0.034847 |
| TOTAL s(FSE) = 2.88% |
══════════════════════════════════════════════════ IMPACTO DE LA OPTIMIZACIÓN DEL PROTOCOLO══════════════════════════════════════════════════ s(FSE) protocolo actual = 10.03% s(FSE) protocolo optimizado = 2.88% Reducción = 7.15 puntos porcentuales Recomendación Pitard (2013): ≤ 2.5% para control metalúrgico Cumplimiento optimizado: NO CUMPLE══════════════════════════════════════════════════Variograma
# ── Datos del variograma (Tabla 3 del artículo) ───────────────────────────────
variog <- data.frame(
lag_j = 0:24,
tiempo = seq(0, 120, by = 5),
v_j = c(0.0047, 0.0045, 0.0039, 0.0038, 0.0041,
0.0040, 0.0053, 0.0053, 0.0059, 0.0041,
0.0061, 0.0066, 0.0060, 0.0067, 0.0062,
0.0079, 0.0069, 0.0067, 0.0075, 0.0066,
0.0076, 0.0070, 0.0077, 0.0060, 0.0085)
)
# ── 1. Nugget v(0): regresión lineal de los 3 primeros puntos a j = 0 ─────────
# Método de Gy (1998) citado en Chieregati et al. (2025):
# el nugget no se observa directamente; se estima extrapolando a j = 0
fit_nugget <- lm(v_j ~ lag_j, data = variog[1:3, ])
v0_nugget <- as.numeric(predict(fit_nugget,
newdata = data.frame(lag_j = 0)))
# ── 2. Sill: varianza estadística global para j > 0 ───────────────────────────
# Cuando el variograma alcanza el sill, los incrementos ya no tienen
# correlación temporal. Es la referencia para la Bias line del nomograma.
sill <- mean(variog$v_j[variog$lag_j > 0])
# ── 3. RSD de las líneas de referencia para el nomograma ─────────────────────
nugget_rsd <- sqrt(v0_nugget) * 100 # Nugget line [%]
bias_rsd <- sqrt(sill) * 100 # Bias line [%]
cat("══════════════════════════════════════════════════\n")══════════════════════════════════════════════════cat(" PARÁMETROS DEL VARIOGRAMA\n") PARÁMETROS DEL VARIOGRAMAcat("══════════════════════════════════════════════════\n")══════════════════════════════════════════════════cat(sprintf(" Nugget v(0) = %.4f → Nugget line = %.2f%%\n",
v0_nugget, nugget_rsd)) Nugget v(0) = 0.0048 → Nugget line = 6.90%cat(sprintf(" Sill = %.4f → Bias line = %.2f%%\n",
sill, bias_rsd)) Sill = 0.0060 → Bias line = 7.77%cat("══════════════════════════════════════════════════\n")══════════════════════════════════════════════════# ── 4. Funciones auxiliares punto a punto (Gy, 1998; Ec. 13–14)
# S(j) = integral simple del variograma
# S'(j) = integral doble del variograma
n <- nrow(variog)
variog$S_j <- NA_real_
variog$Sp_j <- NA_real_
variog$S_j[1] <- 0
variog$Sp_j[1] <- 0
for (i in 2:n) {
variog$S_j[i] <- variog$S_j[i-1] +
0.5 * (variog$v_j[i-1] + variog$v_j[i])
variog$Sp_j[i] <- variog$Sp_j[i-1] +
0.5 * (variog$S_j[i-1] + variog$S_j[i])
}
# w(j) = S(j)/j y w'(j) = 2·S'(j)/j² (Ec. 9 y 11)
variog <- variog |>
mutate(
w_j = ifelse(lag_j > 0, S_j / lag_j, NA_real_),
wp_j = ifelse(lag_j > 0, 2 * Sp_j / lag_j^2, NA_real_)
)
# ── 5. Función generadora de error W(j) = 2·w(j/2) − w'(j) (Ec. 12)
calc_w_half <- function(j, S_vec) {
j0 <- j / 2
if (j %% 2 == 0) {
# j par: w(j0) = S(j0) / j0 (Ec. 15)
S_vec[j0 + 1] / j0
} else {
# j impar: interpolación lineal en j0 + 0.5 (Ec. 16–17)
j0_int <- floor(j0)
S_half <- S_vec[j0_int + 1] +
0.25 * (S_vec[j0_int + 1] + S_vec[j0_int + 2])
2 * S_half / j
}
}
variog$W_j <- NA_real_
for (i in 2:n) {
w_half <- calc_w_half(variog$lag_j[i], variog$S_j)
variog$W_j[i] <- 2 * w_half - variog$wp_j[i]
}
# ── 6. Varianza y desviación estándar del HFE (Ec. 18) ──────────────────────
# s²(HFE) = W(j) / Q
variog <- variog |>
mutate(
Q = c(NA, 36:13),
var_HFE_rel = ifelse(!is.na(Q) & !is.na(W_j), W_j / Q, NA_real_),
sd_HFE_rel = sqrt(abs(var_HFE_rel)) * 100
)
# ── 7. Tabla resumen HFE (comparar con Tabla 3 del artículo) ─────────────────
kable(
variog |>
filter(lag_j >= 1, lag_j <= 12) |>
select(lag_j, tiempo, Q, v_j, var_HFE_rel, sd_HFE_rel),
digits = 6,
col.names = c("Lag j", "Tiempo (min)", "Q",
"v(j)", "s²(HFE) rel.", "s(HFE) rel. (%)"),
caption = "Tabla 3. Varianza y desviación estándar del HFE por intervalo
de muestreo (primeros 12 lags)"
) |>
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) |>
row_spec(1, background = "#e8f5e9", bold = TRUE) # j = 5 min (óptimo)| Lag j | Tiempo (min) | Q | v(j) | s²(HFE) rel. | s(HFE) rel. (%) |
|---|---|---|---|---|---|
| 1 | 5 | 36 | 0.0045 | 0.000000 | 0.000000 |
| 2 | 10 | 35 | 0.0039 | 0.000134 | 1.158817 |
| 3 | 15 | 34 | 0.0038 | 0.000183 | 1.352195 |
| 4 | 20 | 33 | 0.0041 | 0.000137 | 1.168559 |
| 5 | 25 | 32 | 0.0040 | 0.000222 | 1.488707 |
| 6 | 30 | 31 | 0.0053 | 0.000136 | 1.165898 |
| 7 | 35 | 30 | 0.0053 | 0.000239 | 1.545566 |
| 8 | 40 | 29 | 0.0059 | 0.000138 | 1.175816 |
| 9 | 45 | 28 | 0.0041 | 0.000256 | 1.600333 |
| 10 | 50 | 27 | 0.0061 | 0.000143 | 1.196445 |
| 11 | 55 | 26 | 0.0066 | 0.000278 | 1.668298 |
| 12 | 60 | 25 | 0.0060 | 0.000157 | 1.252054 |
Figura 2. Variograma experimental del flujo de alimentación de la planta Boa Vista (reproducción de Chieregati et al., 2025, Fig. 3). El sill determina la Bias line del nomograma.
# ── 8. Gráfico del variograma ─────────────────────────────────────────────────
# Línea de regresión del nugget (extrapolación a j = 0)
x_reg <- seq(0, 2, length.out = 50)
y_reg <- predict(fit_nugget, newdata = data.frame(lag_j = x_reg))
reg_df <- data.frame(tiempo = x_reg * 5, v_reg = y_reg)
ggplot(variog, aes(x = tiempo, y = v_j)) +
# Sill
geom_hline(yintercept = sill,
color = "firebrick", linewidth = 0.8, linetype = "solid") +
annotate("text", x = 122, y = sill + 0.0035,
label = sprintf("Sill = %.4f\n→ Bias line = %.1f%%",
sill, bias_rsd),
color = "firebrick", size = 2.9, hjust = 1) +
# Regresión para estimar nugget
geom_line(data = reg_df,
aes(x = tiempo, y = v_reg),
color = "#1565c0", linewidth = 0.7,
linetype = "dotted", inherit.aes = FALSE) +
# Nugget estimado en j = 0
geom_point(data = data.frame(tiempo = 0, v_j = v0_nugget),
aes(x = tiempo, y = v_j),
shape = 21, fill = "white", color = "#1565c0",
size = 3.5, stroke = 1.4, inherit.aes = FALSE) +
annotate("text", x = 3, y = v0_nugget - 0.0022,
label = sprintf("Nugget v(0) = %.4f\n→ Nugget line = %.1f%%",
v0_nugget, nugget_rsd),
color = "#1565c0", size = 2.9, hjust = 0) +
# Variograma
geom_line(color = "#1565c0", linewidth = 0.9) +
geom_point(color = "#1565c0", size = 2.2, shape = 16) +
# Intervalo óptimo j = 1 (5 min)
geom_vline(xintercept = 5,
linetype = "dotted", color = "#2e7d32", linewidth = 0.8) +
annotate("text", x = 7, y = 0.0081,
label = sprintf("j = 5 min (óptimo)\ns(HFE) = %.2f%%",
variog$sd_HFE_rel[variog$tiempo == 5]),
color = "#2e7d32", size = 2.9, hjust = 0) +
scale_x_continuous(breaks = seq(0, 120, 15), limits = c(0, 128)) +
scale_y_continuous(limits = c(0, 0.0105),
breaks = seq(0, 0.010, 0.001),
labels = number_format(accuracy = 0.001)) +
labs(
title = "Variograma experimental — Alimentación planta Boa Vista",
subtitle = "Niobio | 36 incrementos × 5 min | CMOC Brasil",
x = "Intervalo de tiempo j (min)",
y = "v(j)",
caption = "Fuente: Chieregati et al. (2025), Fig. 3 y Tabla 3"
) +
theme_bw(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
panel.grid.major = element_line(color = "gray88", linewidth = 0.4),
panel.grid.minor = element_line(color = "gray93", linewidth = 0.2)
)
Nomograma de Muestreo
El nomograma traduce la ecuación de Gy a un gráfico log-log que permite leer directamente la masa mínima de muestra para un error objetivo, o el error esperado para una masa dada:
\[\log(M_S) = \log(K) + \alpha \cdot \log(d_N) - \log(s^2_{\text{obj}})\]
# ── NOMOGRAMA FORMATO PITARD (MS en X, RSD en Y) ─────────────────────────────
# Rango de masa de muestra en el eje X
MS_seq_g <- 10^seq(1, 6, length.out = 400)
MS_seq_kg <- MS_seq_g / 1000
# Tamaños de partícula únicos
tamanos <- data.frame(
dN_mm = c(0.104, 2.0, 17.0),
label = c("0.104 mm", "2 mm", "17 mm")
)
# Curvas del nomograma
nomo_data <- do.call(rbind, lapply(seq_len(nrow(tamanos)), function(i) {
dN_cm <- tamanos$dN_mm[i] / 10
IHL <- K_cal * dN_cm^alpha
RSD <- sqrt(IHL / MS_seq_g) * 100
data.frame(
MS_kg = MS_seq_kg,
RSD = RSD,
tamano = tamanos$label[i]
)
}))
nomo_data$tamano <- factor(
nomo_data$tamano,
levels = c("17 mm", "2 mm", "0.104 mm")
)
# Puntos reales del protocolo actual
protocolo <- data.frame(
etapa = c("P", "T", "D", "Pv", "A"),
MS_kg = c(6.000, 6.000, 0.500, 0.500, 0.060),
dN_mm = c(17.0, 2.00, 2.00, 0.104, 0.104)
) |>
mutate(
dN_cm = dN_mm / 10,
IHL = K_cal * dN_cm^alpha,
RSD = sqrt(IHL / (MS_kg * 1000)) * 100
)
# Líneas de referencia
design_rsd <- 2.5
bias_rsd <- sqrt(sill) * 100 # derivado del variograma
nugget_rsd <- sqrt(v0_nugget) * 100 # derivado del variograma
# Colores
colores_lineas <- c(
"0.104 mm" = "#2e7d32",
"2 mm" = "#6a1b9a",
"17 mm" = "#c62828"
)
ggplot(
nomo_data,
aes(x = MS_kg, y = RSD, color = tamano)
) +
geom_line(linewidth = 1) +
# Design line
geom_hline(
yintercept = design_rsd,
linetype = "dashed",
color = "gray30"
) +
annotate(
"text",
x = 800, y = design_rsd * 1.15,
label = "Design line (2.5%)",
hjust = 1, size = 3.2
) +
# Nugget line — sqrt(v(0)) × 100
geom_hline(
yintercept = nugget_rsd,
linetype = "dotted",
color = "#1565c0"
) +
annotate(
"text",
x = 800, y = nugget_rsd * 4,
label = sprintf("Nugget line (%.1f%%) [v(0) = %.4f]",
nugget_rsd, v0_nugget),
hjust = 1, size = 3.2, color = "#1565c0"
) +
# Bias line — sqrt(sill) × 100
geom_hline(
yintercept = bias_rsd,
linetype = "dashed",
color = "gray30"
) +
annotate(
"text",
x = 800, y = bias_rsd * 2,
label = sprintf("Bias line (%.1f%%) [sill = %.4f]",
bias_rsd, sill),
hjust = 1, size = 3.2
) +
geom_point(
data = protocolo,
aes(x = MS_kg, y = RSD),
inherit.aes = FALSE,
shape = 21,
fill = "black",
color = "black",
size = 2
) +
geom_text(
data = protocolo,
aes(x = MS_kg, y = RSD, label = etapa),
inherit.aes = FALSE,
size = 2.7,
vjust = -0.5,
lineheight = 0.85
) +
geom_text(
data = nomo_data |>
group_by(tamano) |>
filter(MS_kg == max(MS_kg)),
aes(label = tamano),
hjust = -0.1,
size = 3,
show.legend = FALSE
) +
scale_x_log10(
breaks = c(0.01, 0.1, 1, 10, 100, 1000),
labels = c("10 g", "100 g", "1 kg", "10 kg", "100 kg", "1,000 kg"),
expand = expansion(mult = c(0.02, 0.18))
) +
scale_y_log10(
breaks = c(1, 2.5, 10, 20, 100, 1000),
labels = c("1%", "2.5%", "10%", "20%", "100%", "1000%"),
sec.axis = sec_axis(
transform = ~.,
breaks = c(1, 2.5, 10, 20, 100, 1000),
labels = c("1%", "2.5%", "10%", "20%", "100%", "1000%"),
name = "(RSD)"
)
) +
scale_color_manual(values = colores_lineas) +
annotation_logticks(
sides = "bl",
color = "gray60",
linewidth = 0.3
) +
labs(
title = "Nomograma FSE — Mina Boa Vista (Niobio)",
subtitle = bquote(
K == .(round(K_cal, 3)) ~
"|" ~ alpha == .(round(alpha, 4))
),
x = "Sample Mass",
y = "Sampling Variance Increments — s(FSE) (%)",
color = "Particle size (dN)",
caption = paste0(
"Puntos negros: protocolo actual | ",
"Design line: Pitard (2013) | ",
"Nugget line: sqrt(v(0))×100 | ",
"Bias line: sqrt(sill)×100 — Chieregati et al. (2025)"
)
) +
theme_bw(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "bottom",
panel.grid.major = element_line(color = "gray80", linewidth = 0.4),
panel.grid.minor = element_line(color = "gray90", linewidth = 0.2)
)