Esta solución implementa una sugerencia hecha por @Innuo en un comentario a la pregunta:
Puede mantener un subconjunto aleatorio muestreado uniformemente de tamaño 100 o 1000 de todos los datos vistos hasta ahora. Este conjunto y las "vallas" asociadas pueden actualizarse en $O(1)$ tiempo.
Una vez que sabemos cómo mantener este subconjunto, podemos seleccionar cualquier nos gusta estimar la media de una población a partir de dicha muestra. Se trata de un método universal, que no hace ninguna suposición, y que funcionará con cualquier de entrada con una precisión que puede predecirse mediante fórmulas de muestreo estadísticas estándar. (La precisión es inversamente proporcional a la raíz cuadrada del tamaño de la muestra).
Este algoritmo acepta como entrada un flujo de datos $x(t),$ $t=1, 2, \ldots,$ un tamaño de muestra $m$ y emite un flujo de muestras $s(t)$ cada uno de los cuales representa a la población $X(t) = (x(1),x(2), \ldots, x(t))$ . En concreto, para $1\le i \le t$ , $s(i)$ es una muestra aleatoria simple de tamaño $m$ de $X(t)$ (sin sustitución).
Para ello, basta con que cada $m$ -subconjunto de elementos de $\{1,2,\ldots,t\}$ tienen las mismas posibilidades de ser los índices de $x$ en $s(t)$ . Esto implica la posibilidad de que $x(i),$ $1\le i\lt t,$ está en $s(t)$ es igual a $m/t$ proporcionó $t \ge m$ .
Al principio sólo recogemos el flujo hasta $m$ elementos han sido almacenados. En ese momento sólo hay una muestra posible, por lo que la condición de probabilidad se cumple trivialmente.
El algoritmo toma el relevo cuando $t=m+1$ . Inductivamente, supongamos que $s(t)$ es una muestra aleatoria simple de $X(t)$ para $t\gt m$ . Fijado provisionalmente $s(t+1) = s(t)$ . Dejemos que $U(t+1)$ sea una variable aleatoria uniforme (independiente de cualquier variable anterior utilizada para construir $s(t)$ ). Si $U(t+1) \le m/(t+1)$ y luego reemplazar un elemento elegido al azar de $s$ por $x(t+1)$ . ¡Ese es todo el procedimiento!
Claramente $x(t+1)$ tiene probabilidad $m/(t+1)$ de estar en $s(t+1)$ . Además, por la hipótesis de inducción, $x(i)$ tenía probabilidad $m/t$ de estar en $s(t)$ cuando $i\le t$ . Con probabilidad $m/(t+1) \times 1/m$ = $1/(t+1)$ se habrá eliminado de $s(t+1)$ por lo que su probabilidad de permanencia es igual a
$$\frac{m}{t}\left(1 - \frac{1}{t+1}\right) = \frac{m}{t+1},$$
exactamente como se necesita. Por inducción, entonces, todas las probabilidades de inclusión del $x(i)$ en el $s(t)$ son correctas y está claro que no hay una correlación especial entre esas inclusiones. Eso demuestra que el algoritmo es correcto.
La eficiencia del algoritmo es $O(1)$ porque en cada etapa se calculan como máximo dos números aleatorios y como máximo un elemento de una matriz de $m$ se sustituyen los valores. La necesidad de almacenamiento es $O(m)$ .
La estructura de datos de este algoritmo consiste en la muestra $s$ junto con el índice $t$ de la población $X(t)$ que muestre. Inicialmente tomamos $s = X(m)$ y proceder con el algoritmo para $t=m+1, m+2, \ldots.$ Aquí hay un R implementación para actualizar $(s,t)$ con un valor $x$ para producir $(s,t+1)$ . (El argumento n desempeña el papel de $t$ y sample.size es $m$ . El índice $t$ será mantenido por la persona que llama).
update <- function(s, x, n, sample.size) {
if (length(s) < sample.size) {
s <- c(s, x)
} else if (runif(1) <= sample.size / n) {
i <- sample.int(length(s), 1)
s[i] <- x
}
return (s)
}
Para ilustrar y comprobar esto, utilizaré el estimador habitual (no robusto) de la media y compararé la media estimada a partir de $s(t)$ a la media real de $X(t)$ (el conjunto acumulado de datos vistos en cada paso). He elegido un flujo de entrada algo difícil que cambia con bastante suavidad pero que periódicamente sufre saltos dramáticos. El tamaño de la muestra de $m=50$ es bastante pequeño, lo que nos permite ver las fluctuaciones del muestreo en estos gráficos.
n <- 10^3
x <- sapply(1:(7*n), function(t) cos(pi*t/n) + 2*floor((1+t)/n))
n.sample <- 50
s <- x[1:(n.sample-1)]
online <- sapply(n.sample:length(x), function(i) {
s <<- update(s, x[i], i, n.sample)
summary(s)})
actual <- sapply(n.sample:length(x), function(i) summary(x[1:i]))
En este punto online es la secuencia de estimaciones medias producidas por el mantenimiento de esta muestra continua de $50$ valores mientras actual es la secuencia de estimaciones medias producidas a partir de todo los datos disponibles en cada momento. El gráfico muestra los datos (en gris), actual (en negro), y dos aplicaciones independientes de este procedimiento de muestreo (en colores). La concordancia está dentro del error de muestreo esperado:
plot(x, pch=".", col="Gray")
lines(1:dim(actual)[2], actual["Mean", ])
lines(1:dim(online)[2], online["Mean", ], col="Red")
![Figure]()
Para los estimadores robustos de la media, busque en nuestro sitio El atípico y términos relacionados. Entre las posibilidades que vale la pena considerar están las medias winsorizadas y los estimadores M.
