Muestreo de depósito

El muestreo de depósito es una familia de algoritmos aleatorios para elegir una muestra aleatoria simple, sin reemplazo, de $k$ elementos de una población de tamaño no conocido $n$ en un solo paso por los elementos. El algoritmo no conoce el tamaño de la población $n$ y, por lo general, es demasiado grande para muestrearla en su totalidad y que todos los $n$ elementos quepan en la memoria principal. La población se revela al algoritmo a lo largo del tiempo y el algoritmo no puede mirar atrás a elementos anteriores. En cualquier momento, el estado actual del algoritmo debe permitir la extracción de una muestra aleatoria simple sin reemplazo de tamaño k sobre la parte de la población vista hasta el momento.

Motivación

Supongamos que vemos una secuencia de elementos, uno a la vez. Queremos mantener diez de estos en la memoria y que se seleccionen al azar en la secuencia. Si conocemos el número total de elementos $n$ y podemos acceder a los elementos arbitrariamente; por lo tanto, la solución es fácil: se seleccionan 10 índices distintos $i$ entre 1 y $n$ con la misma probabilidad y se conservan los $i$ -ésimos elementos. El problema es que no siempre conocemos el $n$ exacto de antemano.

Simple: algoritmo R

Jeffrey Vitter creó un algoritmo simple y popular pero lento, el algoritmo R.^[1]

Inicialice una matriz $R$ indexada de $1$ a $k$ , que contenga los primeros $k$ elementos de la entrada $x_{1},...,x_{k}$ . Este es el depósito .
Por cada nueva entrada $x_{i}$ , genere un número aleatorio $j$ uniformemente en $\{1,...,i\}$ . Si $j\in \{1,...,k\}$ , haga que $R[j]:=x_{i}$ , de lo contrario, descarte $x_{i}$ .
Devuelva $R$ una vez procesadas todas las entradas.

Este algoritmo funciona por inducción en $i\geq k$ .

PruebaCuando i = k el algoritmo R devuelve todas las entradas, proporcionando la base para una demostración por inducción matemática.

Aquí, la hipótesis de inducción es que la probabilidad de que una determinada entrada esté incluida en el depósito justo antes de que se procese la (i + 1)-ésima entrada es k / i y debemos demostrar que la probabilidad de que una entrada determinada esté incluida en el depósito es k / (i + 1) justo después de que se procese la (i + 1)-ésima entrada. Aplique el algoritmo R a la (i + 1)-ésima entrada. La entrada x_i+1 se incluye con probabilidad k / (i + 1) por definición del algoritmo. Para cualquier otra j ∈ {k + 1,…, i + 1}, por la hipótesis de inducción, la probabilidad de que la entrada x_r ∈ {x₁,…, x_i} se incluya en el depósito justo antes de que se procese la (i + 1)-ésima entrada, es k / i.

La probabilidad de que siga incluida en el depósito después de procesar x_i+1 (es decir, de que x_r no sea sustituida por x_i+1) es (k / i)(1 – 1 / (i + 1)). Este último resultado se desprende de la suposición de que el número entero j j se genera uniformemente al azar; una vez que queda claro que de hecho se producirá una sustitución, la probabilidad de que x_r en particular sea sustituido por x_i+1 es (k / (i + 1))(1 / k) = 1 / (i + 1).

Hemos demostrado que la probabilidad de que una entrada nueva entre en el depósito es igual a la probabilidad de que una entrada existente se mantenga en el depósito. Por lo tanto, concluimos por el principio de inducción matemática que el algoritmo R produce efectivamente una muestra aleatoria de las entradas.

Aunque es conceptualmente simple y sencillo de entender, este algoritmo necesita generar un número aleatorio para cada elemento de la entrada, incluidos los elementos que se descartan. Por lo tanto, el tiempo de ejecución asintótico del algoritmo es $O(n)$ . Generar esta cantidad de aleatoriedad y el tiempo de ejecución lineal hace que el algoritmo sea innecesariamente lento cuando la población de entrada es grande.

Este es el Algoritmo R, implementado del siguiente modo:

(* S tiene elementos para muestrear, R contendrá el resultado *)
ReservoirSample(S[1..n], R[1..k])
  // llenar el conjunto de depósitos
  for i := 1 to k
      R[i] := S[i]

  // sustituir elementos con probabilidad gradualmente decreciente
  for i := k+1 to n
    (* randomInteger(a, b) genera un entero uniforme del intervalo inclusivo {a, ..., b} *)
    j := randomInteger(1, i)
    if j <= k
        R[j] := S[i]

(* S tiene elementos para muestrear, R contendrá el resultado *)

MuestraDeDeposito(S[1..n], R[1..k])

// llenar el arreglo del depósito

for i := 1 to k

R[i] := S[i]

// sustituir elementos con probabilidad gradualmente decreciente

for i := k+1 to n

(* randomInteger(a, b) genera un entero uniforme del intervalo inclusivo {a, ..., b} *)

j := randomInteger(1, i)

if j <= k

R[j] := S[i]

Óptimo: algoritmo L

Si generamos $n$ números al azar $u_{1},...,u_{n}\sim U[0,1]$ independientemente, los índices de los $k$ más pequeños es una muestra uniforme de los subconjuntos de k de $\{1,...,n\}$ .

El proceso se puede hacer sin conocer $n$ :

Conserve los $k$ más pequeños de $u_{1},...,u_{i}$ vistos hasta ahora, así como $w_{i}$ , el índice de los más grandes entre ellos. Para cada $u_{i+1}$ nuevo compare con $u_{w_{i}}$ . Si $u_{i+1}<u_{w_{i}}$ , descarte $u_{w_{i}}$ , almacene $u_{i+1}$ y establezca $w_{i+1}$ como el índice del mayor entre ellos. De lo contrario, descarte $u_{i+1}$ y haga $w_{i+1}=w_{i}$ .

Ahora combine esto con el flujo de entradas $x_{1},...,x_{n}$ . Cada vez que alguna $u_{i}$ se acepta, guarde el $x_{i}$ correspondiente. Cada vez que alguna $u_{i}$ se descarta, descarte el $x_{i}$ correspondiente.

Este algoritmo todavía necesita $O(n)$ números aleatorios y toma un tiempo de $O(n)$ . Pero se puede simplificar.

Primera simplificación: no es necesario probar $u_{i+1},u_{i+2},...$ nuevos uno por uno, ya que la probabilidad de que la siguiente aceptación suceda en $u_{i+l}$ es $(1-u_{w_{i}})^{l-1}-(1-u_{w_{i}})^{l}$ , es decir, el intervalo $l$ de aceptación sigue una distribución geométrica.

Segunda simplificación: no es necesario recordar todo el arreglo de los $k$ más pequeños de $u_{1},...,u_{i}$ que se ha visto hasta ahora, sino simplemente $w$ , el mayor de ellos. Esto se basa en tres observaciones:

Cada vez que un $x_{i+1}$ nuevo se selecciona para ir al depósito, se descarta una entrada uniformemente aleatoria en el almacenamiento.
$u_{i+1}\sim U[0,w]$
$w_{i+1}$ tiene la misma distribución que $\max\{u_{1},...,u_{k}\}$ , donde todos $u_{1},...,u_{k}\sim U[0,w]$ de modo independiente. Esto puede ser muestreado con la primera toma de $u\sim U[0,1]$ muestras y tomando $w\cdot u^{\frac {1}{k}}$ .

Este es el algoritmo L,^[2] que se implementa así:

(* S tiene elementos para muestrear, R contendrá el resultado *)
MuestraDeDeposito(S[1..n], R[1..k])
 // llenar el arreglo del depósito
 for i = 1 to k
   R[i] := S[i]

 (* random() genera un número (0,1) aleatorio uniforme *)
 W := exp(log(random())/k)

 while i <= n
   i := i + floor(log(random())/log(1-W)) + 1
   if i <= n
     (* sustituir un elemento aleatorio del depósito por el elemento i *)
     R[randomInteger(1,k)] := S[i] // índice aleatorio comprendido entre 1 y k, ambos inclusive
     W := W * exp(log(random())/k)

Este algoritmo calcula tres números aleatorios para cada elemento que entra a formar parte del depósito y no dedica tiempo a los elementos que no lo hacen. Así que el tiempo de ejecución esperado es $O(k(1+\log(n/k)))$ ,^[2] que es óptimo.^[1] Al mismo tiempo, es fácil de implementar de manera eficiente y no depende de desviaciones aleatorias de distribuciones exóticas o difíciles de calcular.

Con clasificación aleatoria

Si asociamos a cada elemento de la entrada un número aleatorio generado uniformemente, los $k$ elementos con los valores asociados más grandes (o, equivalentemente, más pequeños) forman una muestra aleatoria simple.^[3] Por lo tanto, un muestreo de depósito simple mantiene los $k$ elementos con los valores asociados más grandes en el momento en una cola de prioridad .

(*
S es un flujo de elementos para muestrear
S.Current devuelve el elemento actual del flujo
S.Next avanza el flujo a la siguiente posición
 min-priority-queue permite:
  Count -> número de elementos en la cola de prioridad
  Mínimo -> devuelve el valor de clave mínimo de todos los elementos
  Extract-Min() -> Elimina el elemento con clave mínima
  Insert(key, Item) -> Añade el elemento con la clave especificada
 *)
ReservoirSample(S[1..?])
 H := new min-priority-queue
 while S has data
  r := random()  // aleatorio uniforme entre 0 y 1, exclusivo
  if H.Count < k
   H.Insert(r, S.Current)
  else
   // conservar k elementos con las mayores claves asociadas
   if r > H.Minimum
    H.Extract-Min()
    H.Insert(r, S.Current)
  S.Next
 return items in H

El tiempo de ejecución esperado de este algoritmo es $O(n+k\log k\log(n/k))$ y es relevante en particular porque puede extenderse fácilmente a elementos ponderados.

Muestreo aleatorio ponderado

Este método, también conocido como muestreo secuencial, es incorrecto en el sentido de que no permite obtener probabilidades de inclusión fijadas con anterioridad. Algunas aplicaciones requieren que las probabilidades de muestreo de los elementos estén de acuerdo con los pesos asociados con cada elemento. Por ejemplo, es posible que sea necesario muestrear consultas en un motor de búsqueda, con peso, según la cantidad de veces que se realizaron para que la muestra pueda analizarse con respecto al impacto general en la experiencia del usuario. Sean el peso del elemento $i$ $w_{i}$ y la suma de todos los pesos sea $W$ , hay dos maneras de interpretar los pesos asignados a cada elemento del conjunto:^[4]

En cada ronda, la probabilidad de que cada elemento no seleccionado sea seleccionado en esa ronda sea proporcional a su peso en relación con los pesos de todos los elementos no seleccionados. Si $X$ es la muestra actual, la probabilidad de un elemento $i\notin X$ de ser seleccionado en esa ronda es $\textstyle w_{i}/(W-\sum _{j\in X}{w_{j}})$ .
La probabilidad de que cada elemento sea incluido en la muestra aleatoria es proporcional a su peso relativo, es decir, $w_{i}/W$ . Tenga en cuenta que esta interpretación podría no ser posible en algunos casos, por ejemplo cuando $k=n$ .

Algoritmo A-Res

El siguiente algoritmo fue desarrollado por Efraimidis y Spirakis que usan la interpretación 1:^[5]

(*
 S es un flujo de elementos a muestrear
 S.Actual devuelve el elemento actual del flujo
 S.Peso devuelve el peso del elemento actual del flujo
 S.Next avanza el flujo a la siguiente posición
 El operador de potencia se representa mediante ^
 min-priority-queue permite:
  Count -> número de elementos en la cola de prioridad
  Minimum() -> devuelve el valor de clave mínimo de todos los elementos
  Extract-Min() -> Elimina el elemento con clave mínima
  Insert(key, Item) -> Añade el elemento con la clave especificada
 *)
ReservoirSample(S[1..?])
 H := new min-priority-queue
 while S has data
  r := random() ^ (1/S.Weight)  // random() produce un número aleatorio uniforme en (0,1)
  if H.Count < k
   H.Insert(r, S.Current)
  else
   // conservar los k elementos con las claves asociadas más grandes
   if r > H.Minimum
    H.Extract-Min()
    H.Insert(r, S.Current)
  S.Next
 return items in H

Este algoritmo es idéntico al algoritmo dado en <i>Reservoir Sampling with Random Sort</i> salvo en la generación de las claves de los elementos. El algoritmo equivale a asignar a cada elemento una clave $r^{1/w_{i}}$ donde $r$ es el número aleatorio y se seleccionan los $k$ elementos con las claves más grandes. De manera equivalente, una formulación más estable numéricamente de este algoritmo calcula las claves como $-\ln(r)/w_{i}$ y selecciona los elementos $k$ con las claves más pequeñas.^[6]

Algoritmo A-ExpJ

El siguiente algoritmo es una versión más eficiente de A-Res, también proporcionado por Efraimidis y Spirakis:^[5]

(*
 S es un flujo de elementos a muestrear
 S.Actual devuelve el elemento actual del flujo
 S.Peso devuelve el peso del elemento actual del flujo
 S.Next avanza el flujo a la siguiente posición
 El operador de potencia se representa mediante ^
 min-priority-queue permite:
  Count -> número de elementos en la cola de prioridad
  Mínimo -> clave mínima de cualquier elemento de la cola de prioridad
  Extract-Min() -> Elimina el elemento con clave mínima
  Insert(Key, Item) -> Añade el elemento con la clave especificada
 *)
ReservoirSampleWithJumps(S[1..?])
 H := new min-priority-queue
 while S has data and H.Count < k
  r := random() ^ (1/S.Weight)  // random() produce un número aleatorio uniforme en (0,1)
  H.Insert(r, S.Current)
  S.Next
 X := log(random()) / log(H.Minimum) // esta es la cantidad de peso que hay que saltar
 while S has data
  X := X - S.Weight
  if X <= 0
   t := H.Minimum ^ S.Weight
   r := random(t, 1) ^ (1/S.Weight) // random(x, y) produce un número aleatorio uniforme en (x, y)
  
   H.Extract-Min()
   H.Insert(r, S.Current)

   X := log(random()) / log(H.Minimum)
  S.Next
 return items in H

Este algoritmo sigue las mismas propiedades matemáticas que se utilizan en A-Res, pero en vez de calcular la clave para cada elemento y verificar si ese elemento debe insertarse, calcula un salto exponencial al siguiente elemento que se insertará. Esto evita tener que crear variantes aleatorias para cada elemento, lo que puede resultar costoso. El número de variables aleatorias requeridas se reduce de $O(n)$ a $O(k\log(n/k))$ , donde $k$ es el tamaño del depósito y $n$ es el número de elementos en la secuencia.^[5]

Algoritmo A-Chao

El siguiente algoritmo fue desarrollado por MT Chao quien utiliza la interpretación 2:^[7]

(*
 S tiene elementos para muestrear, R contendrá el resultado
 S[i].Weight contiene el peso de cada elemento
 *)
WeightedReservoir-Chao(S[1..n], R[1..k])
 WSum := 0
 // llenar el arreglo del depósito
 for i := 1 to k
   R[i] := S[i]
   WSum := WSum + S[i].Weight
 for i := k+1 to n
  WSum := WSum + S[i].Weight
  p := S[i].Weight / WSum // probabilidad para este elemento
  j := random();     // aleatorio uniforme entre 0 y 1
  if j <= p        // seleccionar elemento según probabilidad
    R[randomInteger(1,k)] := S[i] //selección uniforme en depósito para sustitución

Para cada elemento, se calcula el peso relativo y se usa para decidir aleatoriamente si el elemento se agregará al depósito. Si se selecciona el elemento, uno de los elementos existentes en el depósito se selecciona uniformemente y se reemplaza con el elemento nuevo. La clave aquí es que, si las probabilidades de todos los elementos en el depósito ya son proporcionales a sus pesos, al seleccionar uniformemente qué elemento reemplazar, las probabilidades de todos los elementos siguen siendo proporcionales a su peso después del reemplazo.

Tenga en cuenta que Chao no especifica cómo muestrear los primeros k elementos. Simplemente asume que tenemos alguna otra forma de recogerlos en proporción a su peso. Chao: «Supongamos que tenemos un plan de muestreo de tamaño fijo con respecto a S_k en el momento A ; tal que su probabilidad de inclusión de primer orden de X_t es π(k; i) ».

Algoritmo A-Chao con saltos

De modo similar a los otros algoritmos, se puede calcular un peso aleatorio j y restar los valores de masa de probabilidad de los elementos, haciéndolos a un lado; mientras que j > 0, reduciendo la cantidad de números aleatorios que deben generarse.^[4]

(*
 S tiene elementos para muestrear, R contendrá el resultado
 S[i].Weight contiene el peso de cada elemento
 *)
WeightedReservoir-Chao(S[1..n], R[1..k])
 WSum := 0
 // llenar el arreglo del depósito
 for i := 1 to k
   R[i] := S[i]
   WSum := WSum + S[i].Weight
 j := random() // aleatorio uniforme entre 0 y 1
 pNone := 1 // probabilidad de que no se haya seleccionado ningún elemento hasta el momento (en este salto)
 for i := k+1 to n
  WSum := WSum + S[i].Weight
  p := S[i].Weight / WSum // probabilidad para este elemento
  j -= p * pNone
  pNone := pNone * (1 - p)
  if j <= 0
    R[randomInteger(1,k)] := S[i] //selección uniforme en el depósito para la sustitución
    j = random()
    pNone := 1

Relación con el algoritmo de Fisher-Yates

Supongamos que se quiere tomar $k$ cartas al azar de una baraja. Un método natural sería barajar las cartas y luego tomar las $k$ cartas superiores. En el caso general, el barajado también debe funcionar incluso si el número de cartas no se conoce de antemano, una condición que se cumple con la versión de adentro hacia afuera del algoritmo Fisher-Yates :^[8]

(* S tiene la entrada, R contendrá la permutación de salida *)
Shuffle(S[1..n], R[1..n])
 R[1] := S[1]
 for i from 2 to n do
   j := randomInteger(1, i) // intervalo inclusivo
   R[i] := R[j]
   R[j] := S[i]

Tenga en cuenta que aunque el resto de las cartas se barajan, solo las primeras $k$ importan en el contexto presente. Por lo tanto, el arreglo $R$ solo necesita rastrear las cartas en las primeras $k$ posiciones mientras realiza el barajado, lo que reduce la cantidad de memoria necesaria. Truncando $R$ a la longitud $k$ , el algoritmo se modifica en consecuencia:

(* S tiene elementos para muestrear, R contendrá el resultado *)
ReservoirSample(S[1..n], R[1..k])
 R[1] := S[1]
 for i from 2 to k do
   j := randomInteger(1, i) // intervalo inclusivo
   R[i] := R[j]
   R[j] := S[i]
 for i from k + 1 to n do
   j := randomInteger(1, i) // intervalo inclusivo
   if (j <= k)
     R[j] := S[i]

Dado que el orden de las primeras $k$ cartas no importa, el primer bucle se puede eliminar y $R$ se puede inicializar para ser los primeros $k$ elementos de la entrada. Esto produce el algoritmo R.

Propiedades estadísticas

La probabilidad de selección de los métodos de depósito se analiza en Chao (1982)^[7] y Tillé (2006).^[9] Mientras que la probabilidad de selección de primer orden es igual a $k/n$ (o, en el caso del procedimiento de Chao, a un conjunto arbitrario de probabilidades desiguales), las probabilidades de selección de segundo orden dependen del orden en que se clasifican los registros en el depósito original. El problema se resuelve con el método de muestreo de cubos de Deville y Tillé (2004).^[10]

Limitaciones

El muestreo de depósito supone que la muestra deseada cabe en la memoria principal, lo que a menudo implica que $k$ es una constante independiente de $n$ . En las aplicaciones en las que nos gustaría seleccionar un gran subconjunto de la lista de entrada (digamos un tercero, es decir, $k=n/3$ ), es necesario adoptar otros métodos. Se han propuesto algunas implementaciones distribuidas para resolver este problema.^[11]

Referencias

↑ ^a ^b Vitter, Jeffrey S. (1 de marzo de 1985). «Random sampling with a reservoir». ACM Transactions on Mathematical Software 11 (1): 37-57. doi:10.1145/3147.3165.
↑ ^a ^b Li, Kim-Hung (4 de diciembre de 1994). «Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n)))». ACM Transactions on Mathematical Software 20 (4): 481-493. doi:10.1145/198429.198435.
↑ Fan, C.; Muller, M.E.; Rezucha, I. (1962). «Development of sampling plans by using sequential (item by item) selection techniques and digital computers». Journal of the American Statistical Association 57 (298): 387-402. doi:10.1080/01621459.1962.10480667.
↑ ^a ^b Efraimidis, Pavlos S. (2015). «Weighted Random Sampling over Data Streams». Algorithms, Probability, Networks, and Games. Lecture Notes in Computer Science 9295: 183-195. ISBN 978-3-319-24023-7. arXiv:1012.0256. doi:10.1007/978-3-319-24024-4_12.
↑ ^a ^b ^c Efraimidis, Pavlos S.; Spirakis, Paul G. (16 de marzo de 2006). «Weighted random sampling with a reservoir». Information Processing Letters 97 (5): 181-185. doi:10.1016/j.ipl.2005.11.003.
↑ Arratia, Richard (2002). «On the amount of dependence in the prime factorization of a uniform random integer». En Bela Bollobas, ed. Contemporary Combinatorics 10: 29-91. ISBN 978-3-642-07660-2.
↑ ^a ^b Chao, M. T. (1982). «A general purpose unequal probability sampling plan». Biometrika 69 (3): 653-656. doi:10.1093/biomet/69.3.653.
↑ National Research Council (2013). Frontiers in Massive Data Analysis (en inglés). The National Academies Press. p. 121. ISBN 978-0-309-28781-4.
↑ Tillé, Yves (2006). Sampling Algorithms. Springer. ISBN 978-0-387-30814-2.
↑ Deville, Jean-Claude; Tillé, Yves (2004). «Efficient balanced sampling: The cube method». Biometrika 91 (4): 893-912. doi:10.1093/biomet/91.4.893.
↑ Reservoir Sampling in MapReduce

[vitter-1] Vitter, Jeffrey S. (1 de marzo de 1985). «Random sampling with a reservoir». ACM Transactions on Mathematical Software 11 (1): 37-57. doi:10.1145/3147.3165.

[li-2] Li, Kim-Hung (4 de diciembre de 1994). «Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n)))». ACM Transactions on Mathematical Software 20 (4): 481-493. doi:10.1145/198429.198435.

[fan-3] Fan, C.; Muller, M.E.; Rezucha, I. (1962). «Development of sampling plans by using sequential (item by item) selection techniques and digital computers». Journal of the American Statistical Association 57 (298): 387-402. doi:10.1080/01621459.1962.10480667.

[efraimidis-4] Efraimidis, Pavlos S. (2015). «Weighted Random Sampling over Data Streams». Algorithms, Probability, Networks, and Games. Lecture Notes in Computer Science 9295: 183-195. ISBN 978-3-319-24023-7. arXiv:1012.0256. doi:10.1007/978-3-319-24024-4_12.

[EfrSpi-5] Efraimidis, Pavlos S.; Spirakis, Paul G. (16 de marzo de 2006). «Weighted random sampling with a reservoir». Information Processing Letters 97 (5): 181-185. doi:10.1016/j.ipl.2005.11.003.

[6] Arratia, Richard (2002). «On the amount of dependence in the prime factorization of a uniform random integer». En Bela Bollobas, ed. Contemporary Combinatorics 10: 29-91. ISBN 978-3-642-07660-2.

[:1-7] Chao, M. T. (1982). «A general purpose unequal probability sampling plan». Biometrika 69 (3): 653-656. doi:10.1093/biomet/69.3.653.

[8] National Research Council (2013). Frontiers in Massive Data Analysis (en inglés). The National Academies Press. p. 121. ISBN 978-0-309-28781-4.

[9] Tillé, Yves (2006). Sampling Algorithms. Springer. ISBN 978-0-387-30814-2.

[10] Deville, Jean-Claude; Tillé, Yves (2004). «Efficient balanced sampling: The cube method». Biometrika 91 (4): 893-912. doi:10.1093/biomet/91.4.893.

[11] Reservoir Sampling in MapReduce

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]