Générer des nombres aléatoires avec une distribution (numérique) donnée

J'ai un fichier avec quelques probabilités pour différentes valeurs par exemple:

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

Je voudrais générer des nombres aléatoires en utilisant cette distribution. Existe-t-il un module existant qui gère cela? Il est assez simple de coder vous-même (construire la fonction de densité cumulative, générer une valeur aléatoire [0,1] et choisir la valeur correspondante) mais il semble que cela devrait être un problème courant et probablement quelqu'un a créé une fonction / module pour il.

J'en ai besoin car je veux générer une liste d'anniversaires (qui ne suivent aucune distribution dans le randommodule standard ).

scipy.stats.rv_discretepourrait être ce que vous voulez. Vous pouvez fournir vos probabilités via le valuesparamètre. Vous pouvez ensuite utiliser la rvs()méthode de l'objet de distribution pour générer des nombres aléatoires.

Comme indiqué par Eugene Pakhomov dans les commentaires, vous pouvez également passer un pparamètre de mot - clé à numpy.random.choice(), par exemple

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

Si vous utilisez Python 3.6 ou supérieur, vous pouvez utiliser à random.choices()partir de la bibliothèque standard - voir la réponse de Mark Dickinson .

Depuis Python 3.6, il existe une solution pour cela dans la bibliothèque standard de Python, à savoir random.choices.

Exemple d'utilisation: définissons une population et des poids correspondant à ceux de la question du PO:

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

Génère maintenant choices(population, weights)un seul échantillon:

>>> choices(population, weights)
4

L'argument facultatif de mot-clé uniquement kpermet de demander plus d'un échantillon à la fois. Ceci est précieux car il y a un travail préparatoire random.choicesà faire à chaque fois qu'il est appelé, avant de générer des échantillons; en générant plusieurs échantillons à la fois, nous n'avons à faire ce travail préparatoire qu'une seule fois. Ici, nous générons un million d'échantillons et utilisons collections.Counterpour vérifier que la distribution que nous obtenons correspond approximativement aux poids que nous avons donnés.

>>> 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})

Un avantage de la génération de la liste à l'aide de CDF est que vous pouvez utiliser la recherche binaire. Bien que vous ayez besoin de temps et d'espace O (n) pour le prétraitement, vous pouvez obtenir k nombres dans O (k log n). Les listes Python normales étant inefficaces, vous pouvez utiliser arraymodule.

Si vous insistez sur un espace constant, vous pouvez faire ce qui suit; O (n) temps, O (1) espace.

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

Il est peut-être un peu tard. Mais vous pouvez utiliser numpy.random.choice(), en passant le pparamètre:

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

(OK, je sais que vous demandez du film rétractable, mais peut-être que ces solutions locales n'étaient tout simplement pas assez succinctes à votre goût. :-)

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)

J'ai pseudo-confirmé que cela fonctionne en regardant la sortie de cette expression:

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

J'ai écrit une solution pour tirer des échantillons aléatoires à partir d'une distribution continue personnalisée .

J'en avais besoin pour un cas d'utilisation similaire au vôtre (c'est-à-dire générer des dates aléatoires avec une distribution de probabilité donnée).

Vous avez juste besoin de la fonction random_custDistet de la ligne samples=random_custDist(x0,x1,custDist=custDist,size=1000). Le reste est de la décoration ^^.

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()

Les performances de cette solution sont certes améliorables, mais je préfère la lisibilité.

Faites une liste d'articles, en fonction de leur 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

Une optimisation peut consister à normaliser les montants par le plus grand diviseur commun, afin de réduire la liste cible.

En outre, cela pourrait être intéressant.

Une autre réponse, probablement plus rapide :)

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]

Vérification:

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

basé sur d'autres solutions, vous générez une distribution cumulative (sous forme d'entier ou de flottant comme vous le souhaitez), puis vous pouvez utiliser la bissectrice pour la rendre rapide

ceci est un exemple simple (j'ai utilisé des entiers ici)

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_cdffonction le convertirait de 20, 60, 10, 10 en 20, 20 + 60, 20 + 60 + 10, 20 + 60 + 10 + 10

maintenant nous choisissons un nombre aléatoire jusqu'à 20 + 60 + 10 + 10 en utilisant random.randintpuis nous utilisons la bissectrice pour obtenir la valeur réelle de manière rapide

vous voudrez peut-être jeter un œil aux distributions d'échantillonnage aléatoire NumPy

Aucune de ces réponses n'est particulièrement claire ou simple.

Voici une méthode claire et simple qui garantit son efficacité.

accumulate_normalize_probabilities prend un dictionnaire pqui mappe les symboles aux probabilités OU aux fréquences. Il génère une liste utilisable de tuples à partir de laquelle effectuer la sélection.

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

Rendements:

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

Pourquoi ça marche

L' étape d' accumulation transforme chaque symbole en un intervalle entre lui-même et la probabilité ou la fréquence des symboles précédents (ou 0 dans le cas du premier symbole). Ces intervalles peuvent être utilisés pour sélectionner (et donc échantillonner la distribution fournie) en parcourant simplement la liste jusqu'à ce que le nombre aléatoire dans l'intervalle 0,0 -> 1,0 (préparé plus tôt) soit inférieur ou égal au point final de l'intervalle du symbole actuel.

La normalisation nous libère du besoin de nous assurer que tout a une certaine valeur. Après normalisation, le «vecteur» des probabilités est égal à 1,0.

Le reste du code pour la sélection et la génération d'un échantillon arbitrairement long à partir de la distribution est ci-dessous:

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

Utilisation:

>>> 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

Voici un moyen plus efficace de procéder:

Appelez simplement la fonction suivante avec votre tableau «poids» (en supposant que les indices sont les éléments correspondants) et le no. d'échantillons nécessaires. Cette fonction peut être facilement modifiée pour gérer une paire ordonnée.

Renvoie les index (ou éléments) échantillonnés / sélectionnés (avec remplacement) en utilisant leurs probabilités respectives:

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

Une brève note sur le concept utilisé dans la boucle while. Nous réduisons le poids de l'élément actuel à partir du bêta cumulatif, qui est une valeur cumulée construite uniformément au hasard, et incrémentons l'indice actuel afin de trouver l'élément dont le poids correspond à la valeur de bêta.