Cómo hacer una muestra aleatoria ponderada de categorías en python

statistics numpy

Question

Dada una lista de tuplas donde cada tupla consiste de una probabilidad y un ítem, me gustaría probar un ítem de acuerdo a su probabilidad. Por ejemplo, proporcione la lista [(.3, ''a''), (.4, ''b''), (.3, ''c'')] Me gustaría probar ''b'' el 40% del tiempo.

¿Cuál es la forma canónica de hacer esto en python?

He observado el módulo aleatorio que no parece tener una función apropiada y at numpy.random que, aunque tiene una función multinomial, no parece devolver los resultados en una buena forma para este problema. Básicamente estoy buscando algo así como mnrnd en matlab.

Muchas gracias.

Gracias por todas las respuestas tan rápido. Para aclarar, no busco explicaciones sobre cómo escribir un esquema de muestreo, sino más bien apuntarme a una forma fácil de tomar muestras de una distribución multinomial dado un conjunto de objetos y pesos, o que me digan que no existe tal función. en una biblioteca estándar, por lo que uno debe escribir uno propio.

Answer 1

¿Cómo crear 3 "a", 4 "b" y 3 "c" en una lista y luego simplemente seleccionar uno al azar. Con suficientes iteraciones obtendrás la probabilidad deseada.

Answer 2

Como nadie usó la función numpy.random.choice , aquí hay una que generará lo que necesita en una única línea compacta:

numpy.random.choice([''a'',''b'',''c''], size = 20, p = [0.3,0.4,0.3])

Answer 3

Creo que la función multinomial es una forma bastante fácil de obtener muestras de una distribución en orden aleatorio. Esto es solo una forma

import numpy
from itertools import izip

def getSamples(input, size):
    probabilities, items = zip(*input)
    sampleCounts = numpy.random.multinomial(size, probabilities)
    samples = numpy.array(tuple(countsToSamples(sampleCounts, items)))
    numpy.random.shuffle(samples)
    return samples

def countsToSamples(counts, items):
    for value, repeats in izip(items, counts):
        for _i in xrange(repeats):
            yield value

Donde las entradas son las especificadas [(.2, ''a''), (.4, ''b''), (.3, ''c'')] y el tamaño es el número de muestras que necesita.

Answer 4

Esto podría hacer lo que quieras:

numpy.array([.3,.4,.3]).cumsum().searchsorted(numpy.random.sample(5))

Answer 5

Esto puede ser de beneficio marginal, pero lo hice de esta manera:

import scipy.stats as sps
N=1000
M3 = sps.multinomial.rvs(1, p = [0.3,0.4,0.3], size=N, random_state=None)
M3a = [ np.where(r==1)[0][0] for r in M3 ] # convert 1-hot encoding to integers

Esto es similar a la respuesta de @ eat.

Answer 6

Existen hacks que puede hacer si, por ejemplo, sus probabilidades se ajustan bien en porcentajes, etc.

Por ejemplo, si está satisfecho con los porcentajes, lo siguiente funcionará (a costa de una gran sobrecarga de memoria):

Pero la forma "real" de hacerlo con probabilidades flotantes arbitrarias es tomar muestras de la distribución acumulativa, después de construirla. Esto es equivalente a subdividir el intervalo de unidad [0,1] en 3 segmentos de línea etiquetados ''a'', ''b'' y ''c''; luego escogiendo un punto al azar en el intervalo de la unidad y viendo qué línea lo segmenta.

#!/usr/bin/python3
def randomCategory(probDict):
    """
        >>> dist = {''a'':.1, ''b'':.2, ''c'':.3, ''d'':.4}

        >>> [randomCategory(dist) for _ in range(5)]
        [''c'', ''c'', ''a'', ''d'', ''c'']

        >>> Counter(randomCategory(dist) for _ in range(10**5))
        Counter({''d'': 40127, ''c'': 29975, ''b'': 19873, ''a'': 10025})
    """
    r = random.random() # range: [0,1)
    total = 0           # range: [0,1]
    for value,prob in probDict.items():
        total += prob
        if total>r:
            return value
    raise Exception(''distribution not normalized: {probs}''.format(probs=probDict))

Uno tiene que tener cuidado con los métodos que devuelven valores incluso si su probabilidad es 0. Afortunadamente este método no lo hace, pero por las dudas, uno podría insertar if prob==0: continue .

Para el registro, aquí está la manera de hackear para hacerlo:

import random

def makeSampler(probDict):
    """
        >>> sampler = makeSampler({''a'':0.3, ''b'':0.4, ''c'':0.3})
        >>> sampler.sample()
        ''a''
        >>> sampler.sample()
        ''c''
    """
    oneHundredElements = sum(([val]*(prob*100) for val,prob in probDict.items()), [])
    def sampler():
        return random.choice(oneHundredElements)
    return sampler

Sin embargo, si no tienes problemas de resolución ... esta es probablemente la manera más rápida posible. =)

Answer 7

No estoy seguro si esta es la manera pitónica de hacer lo que preguntas, pero podrías usar random.sample([''a'',''a'',''a'',''b'',''b'',''b'',''b'',''c'',''c'',''c''],k) donde k es el número de muestras que desea.

Para un método más robusto, biseque el intervalo unitario en secciones en función de la probabilidad acumulada y extraiga de la distribución uniforme (0,1) utilizando random.random (). En este caso, los subintervalos serían (0, .3) (.3, .7) (.7,1). Eliges el elemento en función de en qué subintervalo cae.

Answer 8

Recién inspirado de la sholte muy directa (y correcta) de sholte : voy a demostrar lo fácil que será extenderlo para manejar elementos arbitrarios, como:

In []: s= array([.3, .4, .3]).cumsum().searchsorted(sample(54))
In []: c, _= histogram(s, bins= arange(4))
In []: [item* c[i] for i, item in enumerate(''abc'')]
Out[]: [''aaaaaaaaaaaa'', ''bbbbbbbbbbbbbbbbbbbbbbbbbb'', ''cccccccccccccccc'']

Actualización :
Según los comentarios de phant0m , resulta que se puede implementar una solución aún más sencilla basada en multinomial , como:

In []: s= multinomial(54, [.3, .4, .3])
In []: [item* s[i] for i, item in enumerate(''abc'')]
Out[]: [''aaaaaaaaaaaaaaa'', ''bbbbbbbbbbbbbbbbbbbbbbbbbbb'', ''cccccccccccc'']

En mi humilde opinión, tenemos un buen resumen del empirical cdf y del muestreo multinomial arroja resultados similares. Por lo tanto, en un resumen, recógelo que se adapte mejor a tus propósitos.

Answer 9

import numpy

n = 1000
pairs = [(.3, ''a''), (.3, ''b''), (.4, ''c'')]
probabilities = numpy.random.multinomial(n, zip(*pairs)[0])
result = zip(probabilities, zip(*pairs)[1])
# [(299, ''a''), (299, ''b''), (402, ''c'')]
[x[0] * x[1] for x in result]
# [''aaaaaaaaaa'', ''bbbbbbbbbbbbbbbbbbb'', ''cccccccccccccccccccc'']

¿Cómo exactamente te gustaría recibir los resultados?