Muestreo de depósito

El muestreo de depósito es una familia de algoritmos aleatorios para elegir una muestra aleatoria simple, sin reemplazo, de  elementos de una población de tamaño no conocido  en un solo paso por los elementos. El algoritmo no conoce el tamaño de la población  y, por lo general, es demasiado grande para muestrearla en su totalidad y que todos los  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  y podemos acceder a los elementos arbitrariamente; por lo tanto, la solución es fácil: se seleccionan 10 índices distintos  entre 1 y  con la misma probabilidad y se conservan los  -ésimos elementos. El problema es que no siempre conocemos el  exacto de antemano.

Simple: algoritmo R 
Jeffrey Vitter creó un algoritmo simple y popular pero lento, el algoritmo R.Inicialice una matriz  indexada de  a , que contenga los primeros  elementos de la entrada  . Este es el depósito .

 Por cada nueva entrada , genere un número aleatorio  uniformemente en  . Si , haga que , de lo contrario, descarte  .

 Devuelva  una vez procesadas todas las entradas.Este algoritmo funciona por inducción en  .
 

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 xi+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 xr ∈ {x1,…, xi} 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 xi+1 (es decir, de que xr no sea sustituida por xi+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 xr en particular sea sustituido por xi+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 . 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úmeros al azar  independientemente, los índices de los  más pequeños es una muestra uniforme de los subconjuntos de k de  .

El proceso se puede hacer sin conocer :Conserve los   más pequeños de  vistos hasta ahora, así como , el índice de los más grandes entre ellos. Para cada  nuevo compare con . Si , descarte , almacene  y establezca  como el índice del mayor entre ellos. De lo contrario, descarte  y haga  .Ahora combine esto con el flujo de entradas  . Cada vez que alguna  se acepta, guarde el  correspondiente. Cada vez que alguna  se descarta, descarte el  correspondiente.

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

Primera simplificación: no es necesario probar  nuevos uno por uno, ya que la probabilidad de que la siguiente aceptación suceda en  es , es decir, el intervalo  de aceptación sigue una distribución geométrica.

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

 Cada vez que un  nuevo se selecciona para ir al depósito, se descarta una entrada uniformemente aleatoria en el almacenamiento.
 
  tiene la misma distribución que , donde todos  de modo independiente. Esto puede ser muestreado con la primera toma de muestras y tomando  .

Este es el algoritmo L, 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 , que es óptimo. 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  elementos con los valores asociados más grandes (o, equivalentemente, más pequeños) forman una muestra aleatoria simple. Por lo tanto, un muestreo de depósito simple mantiene los  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 HEl tiempo de ejecución esperado de este algoritmo es  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   y la suma de todos los pesos sea , hay dos maneras de interpretar los pesos asignados a cada elemento del conjunto:

 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  es la muestra actual, la probabilidad de un elemento  de ser seleccionado en esa ronda es  .
 La probabilidad de que cada elemento sea incluido en la muestra aleatoria es proporcional a su peso relativo, es decir, . Tenga en cuenta que esta interpretación podría no ser posible en algunos casos, por ejemplo cuando.

Algoritmo A-Res 
El siguiente algoritmo fue desarrollado por Efraimidis y Spirakis que usan la interpretación 1:(*
 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 HEste 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  donde  es el número aleatorio y se seleccionan los  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  y selecciona los elementos  con las claves más pequeñas.

Algoritmo A-ExpJ 
El siguiente algoritmo es una versión más eficiente de A-Res, también proporcionado por Efraimidis y Spirakis:(*
 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 HEste 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  a , donde  es el tamaño del depósito y  es el número de elementos en la secuencia.

Algoritmo A-Chao 
El siguiente algoritmo fue desarrollado por MT Chao quien utiliza la interpretación 2:(*
 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ónPara 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.(*
 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  cartas al azar de una baraja. Un método natural sería barajar las cartas y luego tomar las  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 :(* 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  importan en el contexto presente. Por lo tanto, el arreglo  solo necesita rastrear las cartas en las primeras  posiciones mientras realiza el barajado, lo que reduce la cantidad de memoria necesaria. Truncando  a la longitud , 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  cartas no importa, el primer bucle se puede eliminar y  se puede inicializar para ser los primeros  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) y Tillé (2006). Mientras que la probabilidad de selección de primer orden es igual a  (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).

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

Referencias 

Algoritmos probabilistas
Análisis de algoritmos
Algoritmos
Algoritmos de aprendizaje automático