Генерация случайных чисел с заданным (числовым) распределением

Question

У меня есть файл с некоторыми вероятностями для различных значений, например:

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

Я хотел бы сгенерировать случайные числа, используя это распределение.Существует ли существующий модуль, который обрабатывает это?Довольно просто написать код самостоятельно (создать функцию кумулятивной плотности, сгенерировать случайное значение [0,1] и выбрать соответствующее значение), но, похоже, это должно быть распространенной проблемой, и, возможно, кто-то создал функцию / модуль дляэто.

Мне это нужно, потому что я хочу сгенерировать список дней рождения (которые не следуют ни за каким распространением в стандартном модуле random).

Answer 1

scipy.stats.rv_discrete может быть тем, что вы хотите. Вы можете указать свои вероятности с помощью параметра values. Затем вы можете использовать метод rvs() объекта распределения для генерации случайных чисел.

Как отметил Евгений Пахомов в комментариях, вы также можете передать ключевой параметр p в numpy.random.choice(), например,

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

Если вы используете Python 3.6 или выше, вы можете использовать random.choices() из стандартной библиотеки - см. Ответ от Mark Dickinson .

Answer 2

Начиная с Python 3.6, есть решение для этого в стандартной библиотеке Python, а именно random.choices.

Пример использования: давайте настроим совокупность и веса, соответствующие тем, которые указаны в вопросе ОП:

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

Теперь choices(population, weights) генерирует один образец:

>>> choices(population, weights)
4

Необязательный аргумент только для ключевого слова k позволяет запрашивать более одного образца одновременно. Это ценно, потому что есть некоторая подготовительная работа, которую random.choices должен выполнять каждый раз, когда он вызывается, перед генерацией каких-либо сэмплов; генерируя много образцов одновременно, мы должны выполнить подготовительную работу только один раз. Здесь мы генерируем миллион выборок и используем collections.Counter, чтобы проверить, что полученное нами распределение примерно соответствует весам, которые мы дали.

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

Answer 3

Преимущество создания списка с использованием CDF заключается в том, что вы можете использовать бинарный поиск. В то время как вам нужно O (n) время и пространство для предварительной обработки, вы можете получить k чисел в O (k log n). Так как обычные списки Python неэффективны, вы можете использовать модуль array.

Если вы настаиваете на постоянном месте, вы можете сделать следующее; O (n) время, O (1) пробел.

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

Answer 4

Может быть, уже поздно.Но вы можете использовать numpy.random.choice(), передав параметр p:

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

Answer 5

(Хорошо, я знаю, что вы просите обертку, но, возможно, этих домашних решений не хватило на ваше усмотрение.: -)

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)

Я псевдо-подтвердил, что это работает, глядя на результат выражения:

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

Answer 6

Я написал решение для рисования случайных выборок из пользовательского непрерывного распределения .

Мне это понадобилось для варианта использования, подобного вашему (т. Е. Генерация случайных дат с заданным распределением вероятности).

Вам просто нужна функция random_custDist и строка samples=random_custDist(x0,x1,custDist=custDist,size=1000). Остальное - украшение ^^.

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

Производительность этого решения наверняка невероятна, но я предпочитаю удобочитаемость.

Answer 7

на основе других решений вы генерируете накопительное распределение (в виде целого числа или числа с плавающей запятой, что хотите), затем вы можете использовать bisect, чтобы сделать это быстрым

это простой пример (здесь я использовал целые числа)

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

функция get_cdf теперь преобразует ее из 20, 60, 10, 10 в 20, 20 + 60, 20 + 60 + 10, 20 + 60 + 10 + 10

сейчасмы выбираем случайное число до 20 + 60 + 10 + 10, используя random.randint, затем используем bisect, чтобы быстро получить действительное значение

Answer 8

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]

Проверка:

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 9

Еще один ответ, возможно, быстрее:)

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 10

Составьте список предметов, основываясь на их 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

Оптимизация может заключаться в нормализации сумм по наибольшему общему делителю для уменьшения целевого списка.

Также, это может быть интересно.