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

Question

Je dispose d'un fichier contenant des 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

J'aimerais générer des nombres aléatoires à l'aide de cette distribution. Existe-t-il un module qui gère cela ? C'est assez simple à coder soi-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 ce soit un problème courant et que quelqu'un ait probablement créé une fonction/module pour cela.

J'ai besoin de cela parce que je veux générer une liste d'anniversaires (qui ne suivent aucune distribution dans le modèle standard). random ).

Answer 1

J'ai écrit une solution pour le tirage d'é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).

Il suffit d'avoir la fonction random_custDist et la ligne samples=random_custDist(x0,x1,custDist=custDist,size=1000) . Le reste n'est que 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()

La performance de cette solution est certainement améliorable, mais je préfère la lisibilité.

Answer 2

Dressez une liste d'éléments, 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 des cibles.

En outre, cette pourrait être intéressante.

Answer 3

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

Answer 4

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

Answer 5

Sur la base d'autres solutions, vous générez une distribution accumulée (sous forme d'entier ou de flottant, comme vous le souhaitez), puis vous pouvez utiliser bisect pour la rendre rapide.

voici un exemple simple (j'ai utilisé des nombres entiers)

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

les get_cdf convertirait 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.randint Nous utilisons ensuite la bissectrice pour obtenir la valeur réelle de manière rapide