Generar números aleatorios con una distribución dada (numérica)

Tengo un archivo con algunas probabilidades para diferentes valores, por ejemplo:

1 0.1
2 0.05
3 0.05
4 0.2
5 0.4
6 0.2

Me gustaría generar números aleatorios usando esta distribución. ¿Existe un módulo existente que maneje esto? Es bastante simple codificar por su cuenta (cree la función de densidad acumulativa, genere un valor aleatorio [0,1] y elija el valor correspondiente), pero parece que esto debería ser un problema común y probablemente alguien ha creado una función / módulo para eso.

Necesito esto porque quiero generar una lista de cumpleaños (que no siguen ninguna distribución en el randommódulo estándar ).

scipy.stats.rv_discretepuede ser lo que quieras Puede proporcionar sus probabilidades a través del valuesparámetro. Luego puede usar el rvs()método del objeto de distribución para generar números aleatorios.

Como señaló Eugene Pakhomov en los comentarios, también puede pasar un pparámetro de palabra clave a numpy.random.choice(), p. Ej.

numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])

Si está utilizando Python 3.6 o superior, puede usarlo random.choices()desde la biblioteca estándar; consulte la respuesta de Mark Dickinson .

Desde Python 3.6, hay una solución para esto en la biblioteca estándar de Python, a saber random.choices.

Ejemplo de uso: configuremos una población y pesos que coincidan con los de la pregunta del OP:

>>> from random import choices
>>> population = [1, 2, 3, 4, 5, 6]
>>> weights = [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]

Ahora choices(population, weights)genera una sola muestra:

>>> choices(population, weights)
4

El argumento opcional de solo palabras clave kpermite solicitar más de una muestra a la vez. Esto es valioso porque hay algo de trabajo preparatorio que random.choicesdebe hacerse cada vez que se llama, antes de generar muestras; Al generar muchas muestras a la vez, solo tenemos que hacer ese trabajo preparatorio una vez. Aquí generamos un millón de muestras y las utilizamos collections.Counterpara verificar que la distribución que obtenemos coincida aproximadamente con los pesos que dimos.

>>> million_samples = choices(population, weights, k=10**6)
>>> from collections import Counter
>>> Counter(million_samples)
Counter({5: 399616, 6: 200387, 4: 200117, 1: 99636, 3: 50219, 2: 50025})

Una ventaja de generar la lista usando CDF es que puede usar la búsqueda binaria. Si bien necesita O (n) tiempo y espacio para el preprocesamiento, puede obtener k números en O (k log n). Como las listas normales de Python son ineficientes, puede usar el arraymódulo.

Si insiste en un espacio constante, puede hacer lo siguiente; O (n) tiempo, O (1) espacio.

def random_distr(l):
    r = random.uniform(0, 1)
    s = 0
    for item, prob in l:
        s += prob
        if s >= r:
            return item
    return item  # Might occur because of floating point inaccuracies

Tal vez sea un poco tarde. Pero puedes usar numpy.random.choice(), pasando el pparámetro:

val = numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])

(OK, sé que estás pidiendo una envoltura retráctil, pero tal vez esas soluciones locales simplemente no fueron lo suficientemente breves para tu gusto. :-)

pdf = [(1, 0.1), (2, 0.05), (3, 0.05), (4, 0.2), (5, 0.4), (6, 0.2)]
cdf = [(i, sum(p for j,p in pdf if j < i)) for i,_ in pdf]
R = max(i for r in [random.random()] for i,c in cdf if c <= r)

Pseudoconfirmé que esto funciona mirando la salida de esta expresión:

sorted(max(i for r in [random.random()] for i,c in cdf if c <= r)
       for _ in range(1000))

Escribí una solución para extraer muestras aleatorias de una distribución continua personalizada .

Necesitaba esto para un caso de uso similar al suyo (es decir, generar fechas aleatorias con una distribución de probabilidad dada).

Solo necesitas la función random_custDisty la línea samples=random_custDist(x0,x1,custDist=custDist,size=1000). El resto es decoración ^^.

import numpy as np

#funtion
def random_custDist(x0,x1,custDist,size=None, nControl=10**6):
    #genearte a list of size random samples, obeying the distribution custDist
    #suggests random samples between x0 and x1 and accepts the suggestion with probability custDist(x)
    #custDist noes not need to be normalized. Add this condition to increase performance. 
    #Best performance for max_{x in [x0,x1]} custDist(x) = 1
    samples=[]
    nLoop=0
    while len(samples)<size and nLoop<nControl:
        x=np.random.uniform(low=x0,high=x1)
        prop=custDist(x)
        assert prop>=0 and prop<=1
        if np.random.uniform(low=0,high=1) <=prop:
            samples += [x]
        nLoop+=1
    return samples

#call
x0=2007
x1=2019
def custDist(x):
    if x<2010:
        return .3
    else:
        return (np.exp(x-2008)-1)/(np.exp(2019-2007)-1)
samples=random_custDist(x0,x1,custDist=custDist,size=1000)
print(samples)

#plot
import matplotlib.pyplot as plt
#hist
bins=np.linspace(x0,x1,int(x1-x0+1))
hist=np.histogram(samples, bins )[0]
hist=hist/np.sum(hist)
plt.bar( (bins[:-1]+bins[1:])/2, hist, width=.96, label='sample distribution')
#dist
grid=np.linspace(x0,x1,100)
discCustDist=np.array([custDist(x) for x in grid]) #distrete version
discCustDist*=1/(grid[1]-grid[0])/np.sum(discCustDist)
plt.plot(grid,discCustDist,label='custom distribustion (custDist)', color='C1', linewidth=4)
#decoration
plt.legend(loc=3,bbox_to_anchor=(1,0))
plt.show()

El rendimiento de esta solución es mejorable con seguridad, pero prefiero la legibilidad.

Haga una lista de artículos, en función de sus weights:

items = [1, 2, 3, 4, 5, 6]
probabilities= [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]
# if the list of probs is normalized (sum(probs) == 1), omit this part
prob = sum(probabilities) # find sum of probs, to normalize them
c = (1.0)/prob # a multiplier to make a list of normalized probs
probabilities = map(lambda x: c*x, probabilities)
print probabilities

ml = max(probabilities, key=lambda x: len(str(x)) - str(x).find('.'))
ml = len(str(ml)) - str(ml).find('.') -1
amounts = [ int(x*(10**ml)) for x in probabilities]
itemsList = list()
for i in range(0, len(items)): # iterate through original items
  itemsList += items[i:i+1]*amounts[i]

# choose from itemsList randomly
print itemsList

Una optimización puede ser normalizar cantidades por el máximo divisor común, para hacer que la lista de objetivos sea más pequeña.

Además, esto podría ser interesante.

Otra respuesta, probablemente más rápido :)

distribution = [(1, 0.2), (2, 0.3), (3, 0.5)]  
# init distribution  
dlist = []  
sumchance = 0  
for value, chance in distribution:  
    sumchance += chance  
    dlist.append((value, sumchance))  
assert sumchance == 1.0 # not good assert because of float equality  

# get random value  
r = random.random()  
# for small distributions use lineair search  
if len(distribution) < 64: # don't know exact speed limit  
    for value, sumchance in dlist:  
        if r < sumchance:  
            return value  
else:  
    # else (not implemented) binary search algorithm  

from __future__ import division
import random
from collections import Counter


def num_gen(num_probs):
    # calculate minimum probability to normalize
    min_prob = min(prob for num, prob in num_probs)
    lst = []
    for num, prob in num_probs:
        # keep appending num to lst, proportional to its probability in the distribution
        for _ in range(int(prob/min_prob)):
            lst.append(num)
    # all elems in lst occur proportional to their distribution probablities
    while True:
        # pick a random index from lst
        ind = random.randint(0, len(lst)-1)
        yield lst[ind]

Verificación:

gen = num_gen([(1, 0.1),
               (2, 0.05),
               (3, 0.05),
               (4, 0.2),
               (5, 0.4),
               (6, 0.2)])
lst = []
times = 10000
for _ in range(times):
    lst.append(next(gen))
# Verify the created distribution:
for item, count in Counter(lst).iteritems():
    print '%d has %f probability' % (item, count/times)

1 has 0.099737 probability
2 has 0.050022 probability
3 has 0.049996 probability 
4 has 0.200154 probability
5 has 0.399791 probability
6 has 0.200300 probability

basado en otras soluciones, genera una distribución acumulativa (como entero o flotante, lo que quiera), luego puede usar bisect para hacerlo más rápido

este es un ejemplo simple (usé enteros aquí)

l=[(20, 'foo'), (60, 'banana'), (10, 'monkey'), (10, 'monkey2')]
def get_cdf(l):
    ret=[]
    c=0
    for i in l: c+=i[0]; ret.append((c, i[1]))
    return ret

def get_random_item(cdf):
    return cdf[bisect.bisect_left(cdf, (random.randint(0, cdf[-1][0]),))][1]

cdf=get_cdf(l)
for i in range(100): print get_random_item(cdf),

la get_cdffunción lo convertiría de 20, 60, 10, 10 en 20, 20 + 60, 20 + 60 + 10, 20 + 60 + 10 + 10

ahora elegimos un número aleatorio hasta 20 + 60 + 10 + 10 usando random.randintluego usamos bisect para obtener el valor real de una manera rápida

es posible que desee echar un vistazo a las distribuciones de muestreo aleatorio NumPy

Ninguna de estas respuestas es particularmente clara o simple.

Aquí hay un método claro y simple que garantiza que funcione.

acumulate_normalize_probabilities toma un diccionario pque asigna símbolos a probabilidades O frecuencias. Produce una lista utilizable de tuplas para hacer la selección.

def accumulate_normalize_values(p):
        pi = p.items() if isinstance(p,dict) else p
        accum_pi = []
        accum = 0
        for i in pi:
                accum_pi.append((i[0],i[1]+accum))
                accum += i[1]
        if accum == 0:
                raise Exception( "You are about to explode the universe. Continue ? Y/N " )
        normed_a = []
        for a in accum_pi:
                normed_a.append((a[0],a[1]*1.0/accum))
        return normed_a

Rendimientos:

>>> accumulate_normalize_values( { 'a': 100, 'b' : 300, 'c' : 400, 'd' : 200  } )
[('a', 0.1), ('c', 0.5), ('b', 0.8), ('d', 1.0)]

Por que funciona

El paso de acumulación convierte cada símbolo en un intervalo entre sí mismo y la probabilidad o frecuencia de los símbolos anteriores (o 0 en el caso del primer símbolo). Estos intervalos se pueden usar para seleccionar (y, por lo tanto, muestrear la distribución proporcionada) simplemente recorriendo la lista hasta que el número aleatorio en el intervalo 0.0 -> 1.0 (preparado anteriormente) sea menor o igual al punto final del intervalo del símbolo actual.

los normalización nos libera de la necesidad de asegurarnos de que todo tenga algún valor. Después de la normalización, el "vector" de probabilidades suma 1.0.

El resto del código para la selección y la generación de una muestra arbitrariamente larga de la distribución se encuentra a continuación:

def select(symbol_intervals,random):
        print symbol_intervals,random
        i = 0
        while random > symbol_intervals[i][1]:
                i += 1
                if i >= len(symbol_intervals):
                        raise Exception( "What did you DO to that poor list?" )
        return symbol_intervals[i][0]


def gen_random(alphabet,length,probabilities=None):
        from random import random
        from itertools import repeat
        if probabilities is None:
                probabilities = dict(zip(alphabet,repeat(1.0)))
        elif len(probabilities) > 0 and isinstance(probabilities[0],(int,long,float)):
                probabilities = dict(zip(alphabet,probabilities)) #ordered
        usable_probabilities = accumulate_normalize_values(probabilities)
        gen = []
        while len(gen) < length:
                gen.append(select(usable_probabilities,random()))
        return gen

Uso:

>>> gen_random (['a','b','c','d'],10,[100,300,400,200])
['d', 'b', 'b', 'a', 'c', 'c', 'b', 'c', 'c', 'c']   #<--- some of the time

Aquí hay una forma más efectiva de hacer esto:

Simplemente llame a la siguiente función con su matriz de 'pesos' (asumiendo los índices como los elementos correspondientes) y el no. de muestras necesarias. Esta función se puede modificar fácilmente para manejar el par ordenado.

Devuelve índices (o elementos) muestreados / seleccionados (con reemplazo) usando sus respectivas probabilidades:

def resample(weights, n):
    beta = 0

    # Caveat: Assign max weight to max*2 for best results
    max_w = max(weights)*2

    # Pick an item uniformly at random, to start with
    current_item = random.randint(0,n-1)
    result = []

    for i in range(n):
        beta += random.uniform(0,max_w)

        while weights[current_item] < beta:
            beta -= weights[current_item]
            current_item = (current_item + 1) % n   # cyclic
        else:
            result.append(current_item)
    return result

Una breve nota sobre el concepto utilizado en el ciclo while. Reducimos el peso del elemento actual de beta acumulativo, que es un valor acumulativo construido de manera uniforme al azar, e incrementamos el índice actual para encontrar el elemento, cuyo peso coincide con el valor de beta.