2015-07-30 1 views
4

У меня есть целое число, которое нужно разделить на ячейки в соответствии с распределением вероятности. Например, если у меня было N=100 объектов, идущих в [0.02, 0.08, 0.16, 0.29, 0.45], тогда вы можете получить [1, 10, 20, 25, 44].Как я могу сгенерировать случайную выборку счетчиков bin, учитывая последовательность вероятностей bin?

import numpy as np 
# sample distribution 
d = np.array([x ** 2 for x in range(1,6)], dtype=float) 
d = d/d.sum() 
dcs = d.cumsum() 
bins = np.zeros(d.shape) 
N = 100 
for roll in np.random.rand(N): 
    # grab the first index that the roll satisfies 
    i = np.where(roll < dcs)[0][0] 
    bins[i] += 1 

В действительности, N и мое число бункеров очень велико, поэтому зацикливания на самом деле не является жизнеспособным вариантом. Есть ли способ, которым я могу векторизовать эту операцию, чтобы ускорить ее?

ответ

3

Вы можете конвертировать PDF в CDF, беря cumsum, использовать это, чтобы определить набор бункеров между 0 и 1, а затем использовать эти контейнеры для вычисления гистограммы в N -длинных случайный однородный вектор:

cdf = np.cumsum([0, 0.02, 0.08, 0.16, 0.29, 0.45])  # leftmost bin edge = 0 
counts, edges = np.histogram(np.random.rand(100), bins=cdf) 

print(counts) 
# [ 4, 8, 16, 30, 42] 
+0

Конечно, я использую термин биннинг и безнадзорность смотреть на гистограмму. Благодаря тонну! – Matt

+0

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

2

Вы можете использовать np.bincount для операции биннинга Alongwith np.searchsorted выполнить эквивалент roll < dcs операции. Вот реализация выполнить эти обещания -

bins = np.bincount(np.searchsorted(dcs,np.random.rand(N),'right')) 

тест выполнения, с использованием заданных параметров -

In [72]: %%timeit 
    ...: for roll in np.random.rand(N): 
    ...:  # grab the first index that the roll satisfies 
    ...:  i = np.where(roll < dcs)[0][0] 
    ...:  bins[i] += 1 
    ...: 
1000 loops, best of 3: 721 µs per loop 

In [73]: %%timeit 
    ...: np.bincount(np.searchsorted(dcs,np.random.rand(N),'right')) 
    ...: 
100000 loops, best of 3: 13.5 µs per loop 
+0

Ваше решение довольно быстро, чем другое, для меньших размеров системы в 10 раз, но медленнее, чем 10 раз медленнее для больших размеров системы. Неплохо, спасибо за альтернативное решение! – Matt