Monitoreo de E. coli en rios del Ecuador: simulacion de una campaña de monitoreo
Luis Guevara, Rodolfo Tixi, Juan Paredes, Jonathan Almeida
2026-10-09
El caso
Parto del articulo de Vinueza y colaboradores (2021), Determining the microbial and chemical contamination in Ecuador’s main rivers, Scientific Reports. Revisaron 12 rios y en todos E. coli y los coliformes totales pasaron el limite. Los mas cargados fueron Zamora (2.50 x 10^4 UFC/100 mL), Machangara (2.25 x 10^4) y Esmeraldas (2.00 x 10^4). El mas bajo que citan es el Coca, con 5.00 x 10^3.
Enlace: https://www.nature.com/articles/s41598-021-96926-z
Con eso se armo una campaña de monitoreo de un rio del mismo tipo. El articulo no ajusta una distribución, así que los parametros los propusimos a partir de lo que si reporta.
La variable discreta es el numero de puntos que superan el limite de E. coli en una campana de 20 puntos. La tome binomial con n = 20 y p = 0.80, porque los 12 rios quedaron por encima del limite y, dentro de un rio ya sucio, lo normal es que la mayoria de puntos tampoco cumpla. No pusimos 1 porque en el mismo cauce puede haber un punto mas limpio. La esperanza es 20 * 0.80 = 16 puntos.
La continua es el logaritmo en base 10 de la concentracion de E. coli, en UFC/100 mL. La tomamos normal con media 4.10 y desviacion 0.25. El 4.10 equivale a unos 1.26 x 10^4 UFC/100 mL, entre el promedio amazonico del articulo (7.24 x 10^3) y el andino (1.88 x 10^4). La desviacion no esta en el paper: la propuse para cubrir el rango citado, de 3.70 a 4.40 en logaritmo.
Como se ven las dos leyes
n_campana <- 20
p <- 0.80
esperanza_x <- n_campana * p
x <- 0:n_campana
par(mfrow = c(1, 2))
plot(x, dbinom(x, size = n_campana, prob = p), type = "h", col = "blue", lwd = 2,
main = "Masa de probabilidad", xlab = "Puntos que no cumplen", ylab = "P(X = x)")
points(x, dbinom(x, size = n_campana, prob = p), pch = 16, col = "blue")
plot(x, pbinom(x, size = n_campana, prob = p), type = "s", col = "red", lwd = 2,
main = "Acumulada", xlab = "Puntos que no cumplen", ylab = "F(x)")
Masa y acumulada de los puntos que no cumplen
par(mfrow = c(1, 1))La masa se junta alrededor de 16. En una campana de este rio lo esperable es que fallen unos 16 de 20 puntos.
media <- 4.10
desv <- 0.25
par(mfrow = c(1, 2))
curve(dnorm(x, mean = media, sd = desv), from = 3.2, to = 5.0, col = "blue", lwd = 2,
main = "Densidad", xlab = "log10 de E. coli", ylab = "f(y)")
abline(v = media, lty = 2, col = "gray40")
curve(pnorm(x, mean = media, sd = desv), from = 3.2, to = 5.0, col = "red", lwd = 2,
main = "Acumulada", xlab = "log10 de E. coli", ylab = "F(y)")
abline(v = media, lty = 2, col = "gray40")
Densidad y acumulada del log10 de E. coli
par(mfrow = c(1, 1))Por que simular
El articulo trae 12 rios, no una campaña de miles de puntos. Simular sirve para ver, antes de ir al campo, como se mueve el promedio si se tomaran 10 puntos o 10000, y en que momento el azar deja de cambiar la lectura.
Hay restricciones. La binomial supone que cada punto cumple o no por su cuenta y con la misma probabilidad. En un rio dos puntos cercanos se parecen. La normal supone que el logaritmo no tiene una cola mas larga que la campana. Y el costo manda: 10 puntos caben en una salida, 10000 no.
Abajo se uso el mismo esquema de los ejercicios de clase.
rbinom y rnorm generan la muestra,
replicate repite el escenario, cumsum arma el
promedio acumulado y la linea roja es la esperanza. Los cortes son 10,
50, 100, 1000 y 10000.
Puntos que no cumplen
tamanos <- c(10, 50, 100, 1000, 10000)
par(mfrow = c(5, 2))
for (B in tamanos) {
muestra <- rbinom(B, size = n_campana, prob = p)
prom_acum <- cumsum(muestra) / (1:B)
hist(muestra, breaks = seq(-0.5, 20.5, by = 1), col = "gray", border = "white",
main = paste("Muestra, B =", B),
xlab = "Puntos que no cumplen", ylab = "Frecuencia")
abline(v = esperanza_x, col = "red", lty = 2, lwd = 2)
plot(prom_acum, type = "l", col = "darkgreen", lwd = 2, log = "x",
main = paste("Convergencia, B =", B),
xlab = "Numero de campanas", ylab = "Promedio acumulado")
abline(h = esperanza_x, col = "red", lty = 2)
}
Discreta: histograma de la muestra y convergencia del promedio. Linea roja en 16
par(mfrow = c(1, 1))set.seed(4321)
for (B in tamanos) {
muestra <- rbinom(B, size = n_campana, prob = p)
print(paste("B =", B, "| promedio =", round(mean(muestra), 4),
"| esperanza =", esperanza_x))
}## [1] "B = 10 | promedio = 15.7 | esperanza = 16"
## [1] "B = 50 | promedio = 16.2 | esperanza = 16"
## [1] "B = 100 | promedio = 15.73 | esperanza = 16"
## [1] "B = 1000 | promedio = 15.9 | esperanza = 16"
## [1] "B = 10000 | promedio = 15.9818 | esperanza = 16"Con B = 10 el histograma es un puñado de barras y la linea verde sube y baja. Una salida de 10 campanas puede dejar el promedio en 14 o en 18 sin que el rio haya cambiado. No tomaria una decision de alerta con eso.
Con B = 50 ya se ve el monton cerca de 16, pero la convergencia todavia zigzaguea. Sirve como primera visita.
Con B = 100 el histograma se parece a la masa y el promedio se queda cerca de 16, con idas y vueltas chicas. Es un tamano trabajable. El azar no se ha apagado.
Con B = 1000 la verde se pega a la roja y ya no se despega. Aqui veo la estabilizacion: el azar deja de producir saltos que cambien la lectura de los 16 puntos.
Con B = 10000 el promedio y la esperanza practicamente coinciden. Confirma lo de 1000. En un rio real ese tamaño casi no se justifica.
Para ver el intervalo, repetimos 5000 veces una campaña de muestreo de 100 puntos y marcamos el 95%.
n <- 100
B <- 5000
promedios_x <- replicate(B, mean(rbinom(n, size = n_campana, prob = p)))
media_sim <- mean(promedios_x)
se_sim <- sd(promedios_x)
z_critico <- qnorm(1 - 0.05 / 2)
ic_inf <- media_sim - z_critico * se_sim
ic_sup <- media_sim + z_critico * se_sim
print(paste("Promedio simulado =", round(media_sim, 4)))## [1] "Promedio simulado = 15.9985"print(paste("Error estandar =", round(se_sim, 4)))## [1] "Error estandar = 0.179"print(paste("IC 95% =", round(ic_inf, 4), "a", round(ic_sup, 4)))## [1] "IC 95% = 15.6477 a 16.3493"hist(promedios_x, breaks = 30, prob = TRUE, col = "gray", border = "white",
main = "Promedio de puntos que no cumplen",
xlab = "Promedio en 100 campanas", ylab = "Densidad")
abline(v = c(ic_inf, ic_sup), col = "blue", lty = 2, lwd = 2)
curve(dnorm(x, mean = esperanza_x, sd = sd(promedios_x)), add = TRUE,
col = "darkred", lwd = 2)
5000 promedios de campanas de 100. Lineas azules: IC 95%
El intervalo no dice que una campaña suelta caiga entre esos dos numeros. Dice que el promedio de 100 campanas, en el 95% de los escenarios, cae ahi, pegado a 16. La curva roja es la normal del teorema del limite central: aunque cada campaña sea binomial, el promedio ya se ve como campaña.
prom_acum <- cumsum(promedios_x) / (1:B)
plot(prom_acum, type = "l", col = "darkgreen", lwd = 2, log = "x",
main = "Convergencia de la perdida de cumplimiento",
xlab = "Numero de escenarios", ylab = "Promedio acumulado")
abline(h = esperanza_x, col = "red", lty = 2)
Convergencia de los 5000 escenarios. Linea roja en 16
log10 de E. coli
par(mfrow = c(5, 2))
for (B in tamanos) {
muestra <- rnorm(B, mean = media, sd = desv)
prom_acum <- cumsum(muestra) / (1:B)
hist(muestra, col = "gray", border = "white",
main = paste("Muestra, B =", B),
xlab = "log10 de E. coli", ylab = "Frecuencia")
abline(v = media, col = "red", lty = 2, lwd = 2)
plot(prom_acum, type = "l", col = "darkgreen", lwd = 2, log = "x",
main = paste("Convergencia, B =", B),
xlab = "Numero de muestras", ylab = "Promedio acumulado")
abline(h = media, col = "red", lty = 2)
}
Continua: histograma de la muestra y convergencia del promedio. Linea roja en 4.10
par(mfrow = c(1, 1))set.seed(4321)
for (B in tamanos) {
muestra <- rnorm(B, mean = media, sd = desv)
print(paste("B =", B, "| promedio =", round(mean(muestra), 4),
"| esperanza =", media))
}## [1] "B = 10 | promedio = 4.1703 | esperanza = 4.1"
## [1] "B = 50 | promedio = 4.1493 | esperanza = 4.1"
## [1] "B = 100 | promedio = 4.0997 | esperanza = 4.1"
## [1] "B = 1000 | promedio = 4.0973 | esperanza = 4.1"
## [1] "B = 10000 | promedio = 4.0997 | esperanza = 4.1"Con B = 10 el histograma no tiene forma de campaña y el promedio puede alejarse una o dos decimas de 4.10. En concentracion eso es bastante: de 4.00 a 4.20 se pasa de 10 000 a unos 16 000 UFC/100 mL. Una campaña corta puede pintar el rio mas limpio o mas sucio de lo que es.
Con B = 50 la forma se insinua y el promedio se acerca, pero una decima todavia se mueve. No alcanza para comparar dos fechas.
Con B = 100 el histograma ya parece la densidad y la verde anda junto a 4.10, con un temblor fino. Es el primer tamano con el que yo armarian un informe de rutina.
Con B = 1000 la linea queda quieta sobre 4.10. Igual que en la discreta, aqui se estabiliza: el azar ya no cambia la concentracion tipica que uno reportaria.
Con B = 10000 el promedio y el 4.10 se confunden. Sirve para ver el limite de la simulacion, no como meta de muestreo.
El intervalo lo armamos, con 150 kits alla y aqui con 150 muestras, repetido 5000 veces.
n <- 150
B <- 5000
promedios_y <- replicate(B, mean(rnorm(n, mean = media, sd = desv)))
p_simulado <- mean(promedios_y)
z_critico2 <- qnorm(1 - 0.05 / 2)
se_y <- desv / sqrt(n)
ic_inf_y <- p_simulado - z_critico2 * se_y
ic_sup_y <- p_simulado + z_critico2 * se_y
print(paste("Promedio simulado =", round(p_simulado, 4)))## [1] "Promedio simulado = 4.1001"print(paste("IC 95% =", round(ic_inf_y, 4), "a", round(ic_sup_y, 4)))## [1] "IC 95% = 4.0601 a 4.1401"hist(promedios_y, breaks = 30, prob = TRUE, col = "gray", border = "white",
main = "Promedio del log10 de E. coli",
xlab = "Promedio en 150 muestras", ylab = "Densidad")
abline(v = c(ic_inf_y, ic_sup_y), lty = 2, col = "red", lwd = 2)
curve(dnorm(x, mean = media, sd = se_y), add = TRUE, col = "blue", lwd = 2)
5000 promedios de 150 muestras. Lineas rojas: IC 95%
El intervalo queda apretado alrededor de 4.10. Una muestra de 150 ya no se va lejos de la concentracion tipica. La curva azul es la normal que predice el teorema del limite central. Aqui el punto de partida ya era normal, asi que el promedio tambien lo es y solo se encoge.
prom_acum_y <- cumsum(promedios_y) / (1:B)
plot(prom_acum_y, type = "l", col = "darkgreen", lwd = 2, log = "x",
main = "Convergencia de la concentracion tipica",
xlab = "Numero de escenarios", ylab = "Promedio acumulado")
abline(h = media, col = "red", lty = 2)
Convergencia de los 5000 escenarios. Linea roja en 4.10
Que se lleva uno de los graficos
En las dos variables el corte en el que el grafico se aquieta es B = 1000. En 10, 50 y 100 todavia se ve la fluctuacion. En 1000 el promedio deja de producir saltos que importen frente a 16 puntos o frente a 4.10. El de 10000 solo confirma, y en campo casi no se justifica.
La trayectoria verde es la ley de los grandes numeros: el promedio de lo observado se acerca a la esperanza cuando entran mas datos, no cuando la muestra es de una sola salida. En este rio, el numero de puntos que no cumplen y la concentracion tipica de E. coli solo se leen con calma si la campaña no es corta.
El teorema del limite central esta en los histogramas de los promedios repetidos. Aunque cada campana discreta sea binomial, el promedio de muchas campanas se acomoda en una campana, mas angosta cuando el tamano crece. En la continua el punto de partida ya es normal y el promedio solo se encoge. Por eso las lineas del intervalo y la curva teorica casi se pisan cuando la muestra ya no es de 10.
La lectura practica: para un rio como los del articulo, una campana de 10 no alcanza para decir que bajo la contaminacion. Alrededor de 1000 el promedio ya no se mueve por casualidad.