Ответ 1
Прежде всего, я считаю, что проблема, с которой вы сталкиваетесь, связана с неправильной нормализацией ваших вероятностей. Эта строка неверна:
a = np.exp(l) / scipy.misc.logsumexp(l)
Вы делите вероятность на логарифмическую вероятность, которая не имеет смысла. Вместо этого вы, вероятно, захотите
a = np.exp(l - scipy.misc.logsumexp(l))
Если вы это сделаете, вы найдете a = [1, 0], и ваш многокомпонентный сэмплер работает как ожидалось до точности с плавающей запятой во второй вероятности.
Решение для малых N: гистограммы
Тем не менее, если вам все еще нужна более высокая точность и производительность, это не так важно, так как вы могли бы добиться прогресса, это реализовать мультиномиальный сэмплер с нуля, а затем модифицировать его для более высокой точности.
Многочленная функция NumPy реализована в Cython и по существу выполняет цикл над несколькими биномиальными образцами и объединяет их в многочленную выборку. и вы можете назвать это следующим образом:
np.random.multinomial(10, [0.1, 0.2, 0.7])
# [0, 1, 9]
(Обратите внимание, что точные выходные значения здесь и ниже являются случайными и будут меняться от вызова к вызову).
Другой способ, которым вы могли бы реализовать мультиномиальный сэмплер, - генерировать N равномерных случайных значений, а затем вычислить гистограмму с ячейками, определенными кумулятивными вероятностями:
def multinomial(N, p):
rand = np.random.uniform(size=N)
p_cuml = np.cumsum(np.hstack([[0], p]))
p_cuml /= p_cuml[-1]
return np.histogram(rand, bins=p_cuml)[0]
multinomial(10, [0.1, 0.2, 0.7])
# [1, 1, 8]
С учетом этого метода мы можем думать о том, чтобы делать вещи с большей точностью, сохраняя все в лог-пространстве. Основной трюк заключается в том, чтобы понять, что журнал равномерных случайных отклонений эквивалентен отрицательному экспоненциальному случайному отклонению, и поэтому вы можете сделать все выше, не покидая пространства журнала:
def multinomial_log(N, logp):
log_rand = -np.random.exponential(size=N)
logp_cuml = np.logaddexp.accumulate(np.hstack([[-np.inf], logp]))
logp_cuml -= logp_cuml[-1]
return np.histogram(log_rand, bins=logp_cuml)[0]
multinomial_log(10, np.log([0.1, 0.2, 0.7]))
# [1, 2, 7]
Полученные многолинейные ничьи будут поддерживать точность даже при очень малых значениях в массиве p.
К сожалению, эти решения на основе гистограммы будут намного медленнее, чем функция native numpy.multinomial, поэтому, если производительность является проблемой, вам может понадобиться другой подход. Один из вариантов - адаптировать приведенный выше Cython-код для работы в лог-пространстве, используя аналогичные математические трюки, которые я использовал здесь.
Решение для большого N: приближение Пуассона
Проблема с приведенным выше решением состоит в том, что при увеличении N он становится очень медленным.
Я подумал об этом и понял там более эффективный путь вперед, несмотря на то, что np.random.multinomial не удалось получить вероятности меньше 1E-16 или так.
Вот пример этого сбоя: на 64-битной машине это всегда будет давать ноль для первой записи из-за способа реализации кода, когда на самом деле он должен дать что-то около 10:
np.random.multinomial(1E18, [1E-17, 1])
# array([ 0, 1000000000000000000])
Если вы вникнете в источник, вы можете проследить эту проблему с помощью биномиальной функции, на которой построена многочлена. Внутренний код cython делает следующее:
def multinomial_basic(N, p, size=None):
results = np.array([np.random.binomial(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_basic(1E18, [1E-17, 1])
# array([ 0, 1000000000000000000])
Проблема заключается в том, что функция binomial дросселируется при очень малых значениях p - это потому, что алгоритм вычисляет значение (1 - p), поэтому значение p ограничено точностью с плавающей точкой.
Итак, что мы можем сделать? Ну, оказывается, что при малых значениях р распределение Пуассона является чрезвычайно хорошим приближением биномиального распределения, и реализация не имеет этих проблем. Таким образом, мы можем построить надежную многочленную функцию, основанную на надежном биномиальном сэмплере, который переключается на пробоотборник Пуассона при малых значениях p:
def binomial_robust(N, p, size=None):
if p < 1E-7:
return np.random.poisson(N * p, size)
else:
return np.random.binomial(N, p, size)
def multinomial_robust(N, p, size=None):
results = np.array([binomial_robust(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_robust(1E18, [1E-17, 1])
array([ 12, 999999999999999988])
Первая запись отлична от нуля и около 10, как и ожидалось! Обратите внимание: мы не можем использовать N больше, чем 1E18, потому что он переполнит длинное целое число.
Но мы можем подтвердить, что наш подход работает для меньших вероятностей с использованием параметра size и усреднения по результатам:
p = [1E-23, 1E-22, 1E-21, 1E-20, 1]
size = int(1E6)
multinomial_robust(1E18, p, size).mean(0)
# array([ 1.70000000e-05, 9.00000000e-05, 9.76000000e-04,
# 1.00620000e-02, 1.00000000e+18])
Мы видим, что даже при этих очень малых вероятностях многочленные значения появляются в правильной пропорции. В результате получается очень надежное и очень быстрое приближение к многочленному распределению при малых p.