Блог

Допустим, у меня есть собственный пользовательский вектор распределения numpy p .

Тогда p удовлетворяет следующему:

 np.ndim(p) == 1 amp; np.sum(p) == 1 amp; np.all(p >= 0)
 

С помощью этого вектора я могу легко выбрать число [0, p.shape) с помощью np.random.choice(np.arange(len(p)), p=p)

В случае, когда у меня много таких p s, у меня есть матрица (с dim 2) P , которая удовлетворяет:

 np.sum(P[:,i]) == 1   # for all i in P.shape[1]
np.all(P >= 0)
 

Затем я хочу выбрать P.shape[1] числа в диапазоне от 0 до P.shape[0] с вероятностью P.

Например, следующий код:

 P = np.array([[0.2, 0.3],
              [0.5, 0.7],
              [0.3, 0]])
x = np.random.choice(np.arange(P.shape[0], P[:,0]))
y = np.random.choice(np.arange(P.shape[0], P[:,1]))
 

произведет мою волю ( x=0 in 0.2 , x=1 in 0.5 и x=2 in 0.3 и y=0 in 0.3 , y=1 in 0.7 ).

В моем случае P много столбцов, и я хочу выбрать все за один раз.

Конечно, я могу сделать это в цикле for, например:

 random_values = np.empty(P.shape[1])
arange_arr = np.arange(P.shape[0])
for i in range(P.shape[1]):
    random_values[i] = np.random.choice(arange_arr, p=P[:,i])
 

Пытаюсь найти какой-нибудь изящный элегантный способ сделать это.

Вы могли бы сделать что-то вроде этого:

 P = np.array([[0.2, 0.3],
              [0.5, 0.7],
              [0.3, 0]])
P_upper = np.cumsum(P, axis=0)
P_lower = np.concatenate((np.zeros((1, P.shape[1])), P_upper[:-1, :]), axis=0)
 

Это создает набор ячеек, которые вы можете оцифровать. Теперь генерируйте случайные числа от 0 до 1:

 r = np.random.rand(10, P.shape[1])
 

Есть несколько способов назначить данные в правильные ячейки. Быстрый и относительно неэффективный способ — использовать логическую маску:

 mask = (r[None, ...] >= P_lower[:, None, :]) amp; (r[None, ...] < P_upper[:, None, :])
result = np.argmax(mask, axis=0)
 

Более эффективный, но более сложный способ — добавить смещение к каждому столбцу и применить np.digitize или np.searchsorted к результату:

 offset = np.arange(P.shape[1])
ind = np.searchsorted((P_upper   offset).ravel('F'), (r   offset).ravel('F')).reshape(r.shape, order='F')
result = ind - offset * P.shape[0]
 

TL; DR

 def multi_sample(p, n):
    ps = np.cumsum(p, axis=0)
    r = np.random.rand(n, ps.shape[1])
    offset = np.arange(P.shape[1])
    ind = np.searchsorted((P_upper   offset).ravel('F'), (r   offset).ravel('F')).reshape(r.shape, order='F')
    return ind - offset * P.shape[0]