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


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

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

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

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

12 79

12 ответов:

scipy.stats.rv_discrete может быть то, что вы хотите. Вы можете предоставить свои вероятности через

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

пример использования: давайте настроим популяцию и веса, соответствующие тем, которые находятся в вопросе 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]

теперь 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})

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

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

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

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

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

возможно, вы захотите взглянуть на NumPy распределения случайной выборки

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

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

и этой может быть интересно.

другой ответ, вероятно, быстрее :)

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]

проверка:

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

основываясь на других решениях, вы генерируете накопительное распределение (как целое или плавающее, что вам нравится), затем вы можете использовать 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),

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

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

ни один из этих ответов не является особенно ясным или простым.

здесь ясный, простой метод, который гарантированно работает.

accumulate_normalize_probabilities берет словарь p что сопоставляет символы с вероятностями или частот. Он выводит полезный список кортежей, из которых нужно сделать выбор.

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

выходы:

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

почему он работает

в накопление шаг превращает каждый символ в интервал между собой и предыдущими символами вероятности или частоты (или 0 в случае первого символа). Эти интервалы можно использовать для выбора из (и, таким образом, выборки предоставленного распределения), просто шагая по списку, пока случайное число в интервале 0.0 -> 1.0 (подготовленное ранее) не станет меньше или равно конечной точке интервала текущего символа.

The нормализация освобождает нас от необходимость убедиться, что все складывается в некоторую ценность. После нормализации "вектор" вероятностей суммируется до 1,0.

The остальная часть кода для выбора и генерации произвольно длинной выборки из распределения ниже:

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

использование :

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

здесь более эффективным способом для этого:

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

возвращает индексы (или элементы), отобранные/выбранные (с заменой), используя их соответствующие вероятности:

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

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