¿Algoritmo eficiente para el muestreo aleatorio de una distribución al tiempo que permite actualizaciones?

Question

Esta es la pregunta que me hicieron hace un tiempo en una entrevista, no pude encontrar la respuesta.

Dadas algunas muestras S1, S2, ... Sn y sus distribuciones de probabilidad (o ponderaciones, como se llame) P1, P2, .. Pn, algoritmo de diseño que elige aleatoriamente la muestra teniendo en cuenta su probabilidad. la solución Vine con es como sigue:

1. Build formación acumulativa de pesos Ci, tales

   C0 = 0; Ci = C [i-1] + Pi.

al mismo tiempo calcular T = P1 + P2 + ... Pn. Se necesita tiempo O (n)

2. Generar número uniformemente aleatorio R = T * aleatorio [0..1]
3. utilizando el algoritmo de búsqueda binaria, vuelva menos i tales Ci> = R. resultado es Si. Toma el tiempo O (logN).

Ahora la pregunta real es: Supongamos que deseo cambiar uno de los pesos Pj iniciales. cómo hacer esto en un tiempo mejor que O (n)? otras estructuras de datos son aceptables, pero el algoritmo de muestreo aleatorio no debería ser peor que O (logN).

Answer 1

Una forma de resolver esto es repensar cómo se construye el árbol de búsqueda binaria que contiene los totales acumulados. En lugar de construir un árbol binario de búsqueda, pensar en tener cada nodo interpretado de la siguiente manera:

• cada nodo almacena una gama de valores que se dedica al propio nodo.
• Los nodos en el subárbol izquierdo representan el muestreo de la distribución de probabilidad justo a la izquierda de ese rango.
• Los nodos en el subárbol derecho representan el muestreo de la distribución de probabilidad justo a la derecha de ese rango.

Por ejemplo, supongamos que nuestros pesos son 3, 2, 2, 2, 1 y 1 para los eventos A, B, C, D, E, F y G. Construimos este árbol binario A, B, C, D, e, F y G:

   D 
      / \ 
      B  F 
     /\ /\ 
     A C E G 

Ahora, anotar el árbol con probabilidades.Dado que A, C, E y G son todas las hojas, damos a cada uno de ellos una masa de probabilidad:

   D 
      / \ 
      B  F 
     /\ /\ 
     A C E G 
     1 1 1 1 

Ahora, mira el árbol de B. B tiene un peso de ser elegido 2, A tiene un peso 3 de ser elegido, y C tiene la probabilidad 2 de ser elegido. Si los normalizamos en el rango [0, 1), A representa 3/7 de la probabilidad y B y C representan cada uno 2/7s. Así tenemos el nodo para B decir que cualquier cosa en el rango [0, 3/7) va al subárbol izquierdo, cualquier cosa en el rango [3/7, 5/7) se asigna a B, y cualquier cosa en el rango [5/7, 1) se asigna al subárbol derecho:

    D 
       / \ 
      B    F 
[0, 3/7)/\ [5/7, 1)/\ 
     A C   E G 
     1 1   1 1 

mismo modo, vamos a procesar F. e tiene un peso 2 de ser elegido, mientras que F y G cada uno tienen un peso probabilidad 1 de ser elegido. Por lo tanto, el subárbol para E representa 1/2 de la masa de probabilidad aquí, el nodo F da cuenta de 1/4 y el subárbol de cuentas G para 1/4. Esto significa que podemos asignar probabilidades como

     D 
        / \ 
      B      F 
[0, 3/7)/\ [5/7, 1) [0, 1/2)/\ [3/4, 1) 
     A C     E G 
     1 1     1 1 

Finalmente, echemos un vistazo a la raíz. El peso combinado del subárbol izquierdo es 3 + 2 + 2 = 7. El peso combinado del subárbol derecho es 2 + 1 + 1 = 4. El peso de D en sí es 2. Por lo tanto, el subárbol izquierdo tiene probabilidad 7/13 de siendo escogido, D tiene una probabilidad 2/13 de ser elegido, y el subárbol correcto tiene una probabilidad 4/13 de ser elegido. De este modo podemos finalizado las probabilidades como

     D 
      [0, 7/13)/ \ [9/13, 1) 
      B      F 
[0, 3/7)/\ [5/7, 1) [0, 1/2)/\ [3/4, 1) 
     A C     E G 
     1 1     1 1 

Para generar un valor aleatorio, sería repetir lo siguiente:

• A partir de la raíz:
  • Elija un valor uniforme-aleatorio en el rango [0, 1).
  • Si está en el rango para el subárbol izquierdo, bajar en ella.
  • Si está en el rango para el subárbol derecho, bajar en ella.
  • De lo contrario, devuelva el valor correspondiente al nodo actual.

Las probabilidades sí se puede determinar de forma recursiva cuando el árbol se construye:

• El probabilidades izquierdo y derecho son 0 para cualquier nodo hoja.
• Si un mismo nodo interior tiene un peso W, su árbol de la izquierda tiene peso total W _L, y su árbol de la derecha tiene peso total W _R, entonces la probabilidad izquierda es (W _L)/(W + W _L + W _R) y la probabilidad derecha es (W _R)/(W + W _L + W _R).

La razón de que esta reformulación sea útil es que nos da una forma de actualizar las probabilidades en el tiempo O (log n) por probabilidad actualizada. En particular, pensemos en qué invariantes van a cambiar si actualizamos el peso de algún nodo en particular. Para simplificar, supongamos que el nodo es una hoja por ahora.Cuando actualizamos el peso del nodo hoja, las probabilidades siguen siendo correctas para el nodo hoja, pero son incorrectas para el nodo que está justo encima, porque el peso de uno de los subárboles de ese nodo ha cambiado. Por lo tanto, podemos (en tiempo O (1)) volver a calcular las probabilidades para el nodo padre simplemente usando la misma fórmula que la anterior. Pero entonces el padre de ese nodo ya no tiene los valores correctos porque uno de sus pesos de subárbol ha cambiado, por lo que también podemos recalcular la probabilidad allí. Este proceso se repite hasta la raíz del árbol, con nosotros haciendo O (1) cálculo por nivel para rectificar los pesos asignados a cada borde. Suponiendo que el árbol está equilibrado, por lo tanto, tenemos que hacer el trabajo total O (log n) para actualizar una probabilidad. La lógica es idéntica si el nodo no es un nodo hoja; simplemente comenzamos en algún lugar del árbol.

En resumen, esto da

tiempo

• O (n) para construir el árbol (utilizando un enfoque de abajo arriba),
• O (log n) tiempo para generar un valor aleatorio, y
• O (log n) tiempo para actualizar cualquier valor.

Espero que ayude!

Answer 2

En lugar de una matriz, almacene la búsqueda estructurada como un árbol binario equilibrado. Cada nodo del árbol debe almacenar el peso total de los elementos que contiene. Dependiendo del valor de R, el procedimiento de búsqueda devuelve el nodo actual o busca a través del subárbol izquierdo o derecho.

Cuando se cambia el peso de un elemento, la actualización de la estructura de búsqueda es una cuestión de ajustar los pesos en la ruta del elemento a la raíz del árbol.

Dado que el árbol está equilibrado, las operaciones de búsqueda y actualización de peso son ambas O (log N).

Answer 3

Para aquellos de ustedes que quieran algo de código, he aquí una implementación de Python:

import numpy 


class DynamicProbDistribution(object): 
    """ Given a set of weighted items, randomly samples an item with probability 
    proportional to its weight. This class also supports fast modification of the 
    distribution, so that changing an item's weight requires O(log N) time. 
    Sampling requires O(log N) time. """ 

    def __init__(self, weights): 
    self.num_weights = len(weights) 
    self.weights = numpy.empty((1+len(weights),), 'float32') 
    self.weights[0] = 0 # Not necessary but easier to read after printing 
    self.weights[1:] = weights 
    self.weight_tree = numpy.zeros((1+len(weights),), 'float32') 
    self.populate_weight_tree() 

    def populate_weight_tree(self): 
    """ The value of every node in the weight tree is equal to the sum of all 
    weights in the subtree rooted at that node. """ 
    i = self.num_weights 
    while i > 0: 
     weight_sum = self.weights[i] 
     twoi = 2*i 
     if twoi < self.num_weights: 
     weight_sum += self.weight_tree[twoi] + self.weight_tree[twoi+1] 
     elif twoi == self.num_weights: 
     weight_sum += self.weights[twoi] 
     self.weight_tree[i] = weight_sum 
     i -= 1 

    def set_weight(self, item_idx, weight): 
    """ Changes the weight of the given item. """ 
    i = item_idx + 1 
    self.weights[i] = weight 
    while i > 0: 
     weight_sum = self.weights[i] 
     twoi = 2*i 
     if twoi < self.num_weights: 
     weight_sum += self.weight_tree[twoi] + self.weight_tree[twoi+1] 
     elif twoi == self.num_weights: 
     weight_sum += self.weights[twoi] 
     self.weight_tree[i] = weight_sum 
     i /= 2 # Only need to modify the parents of this node 

    def sample(self): 
    """ Returns an item index sampled from the distribution. """ 
    i = 1 
    while True: 
     twoi = 2*i 

     if twoi < self.num_weights: 
     # Two children 
     val = numpy.random.random() * self.weight_tree[i] 
     if val < self.weights[i]: 
      # all indices are offset by 1 for fast traversal of the 
      # internal binary tree 
      return i-1 
     elif val < self.weights[i] + self.weight_tree[twoi]: 
      i = twoi # descend into the subtree 
     else: 
      i = twoi + 1 

     elif twoi == self.num_weights: 
     # One child 
     val = numpy.random.random() * self.weight_tree[i] 
     if val < self.weights[i]: 
      return i-1 
     else: 
      i = twoi 

     else: 
     # No children 
     return i-1 


def validate_distribution_results(dpd, weights, samples_per_item=1000): 
    import time 

    bins = numpy.zeros((len(weights),), 'float32') 
    num_samples = samples_per_item * numpy.sum(weights) 

    start = time.time() 
    for i in xrange(num_samples): 
    bins[dpd.sample()] += 1 
    duration = time.time() - start 

    bins *= numpy.sum(weights) 
    bins /= num_samples 

    print "Time to make %s samples: %s" % (num_samples, duration) 

    # These should be very close to each other 
    print "\nWeights:\n", weights 
    print "\nBins:\n", bins 

    sdev_tolerance = 10 # very unlikely to be exceeded 
    tolerance = float(sdev_tolerance)/numpy.sqrt(samples_per_item) 
    print "\nTolerance:\n", tolerance 

    error = numpy.abs(weights - bins) 
    print "\nError:\n", error 

    assert (error < tolerance).all() 


#@test 
def test_DynamicProbDistribution(): 
    # First test that the initial distribution generates valid samples. 

    weights = [2,5,4, 8,3,6, 6,1,3, 4,7,9] 
    dpd = DynamicProbDistribution(weights) 

    validate_distribution_results(dpd, weights) 

    # Now test that we can change the weights and still sample from the 
    # distribution. 

    print "\nChanging weights..." 
    dpd.set_weight(4, 10) 
    weights[4] = 10 
    dpd.set_weight(9, 2) 
    weights[9] = 2 
    dpd.set_weight(5, 4) 
    weights[5] = 4 
    dpd.set_weight(11, 3) 
    weights[11] = 3 

    validate_distribution_results(dpd, weights) 

    print "\nTest passed" 


if __name__ == '__main__': 
    test_DynamicProbDistribution()