np.random.dirichlet с небольшим параметром: встроить будущее решение в текущий numpy

Question

Происходит постоянное обсуждение np.random.dirichlet функция, так как она не работает для небольших параметров:

In [1]: import numpy as np

In [2]: np.random.dirichlet(np.ones(3)*.00001)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-2-464b0fe9c6c4> in <module>()
----> 1 np.random.dirichlet(np.ones(3)*.00001)

mtrand.pyx in mtrand.RandomState.dirichlet (numpy/random/mtrand/mtrand.c:25213)()

mtrand.pyx in mtrand.RandomState.dirichlet (numpy/random/mtrand/mtrand.c:25123)()

ZeroDivisionError: float division

Обсуждение можно прочитать здесь и здесь, и указывает, что это ошибка нормализации. В настоящее время предлагаемое усовершенствование для переключения сэмплеров для небольших параметров не может быть объединено с master of numpy по нескольким причинам.

Вопрос: Может ли кто-нибудь предложить другой способ рисования дирихлетов в python или указать мне решение использовать новый сэмплер без перекомпиляции моего numpy и / или работы над невыпущенной веткой?

Answer 1

Хорошо, давайте попробуем следующее. Вот бета (альфа, бета) вариабельная выборка, которая должна работать для любых небольших чисел.

import math
import random

def sample_beta(alpha, beta):
    x = math.log( random.random() )
    y = math.log( random.random() )

    return x / (x + y*alpha/beta)

# some testing
import matplotlib.pyplot as plt

bins = [0.01 * i for i in range(102)]
plt.hist([sample_beta(0.00001, 0.1) for k in range(10000000)], bins)
plt.show()

Используя его, вы можете попробовать попробовать Dirichlet с помощью бета-вариации, как описано в википедии.

https://en.wikipedia.org/wiki/Dirichlet_distribution

params = [a1, a2, ..., ak]
xs = [sample_beta(params[0], sum(params[1:]))]
for j in range(1,len(params)-1):
    phi = sample_beta(params[j], sum(params[j+1:]))
    xs.append((1-sum(xs)) * phi)
xs.append(1-sum(xs))

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

ОБНОВИТЬ

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

beta(a, b) = gamma(1, a) / (gamma(1, a) + gamma(1, b))

Таким образом, небольшие параметры перемещаются из первого в гамме (если вы производите выборку Дирихле непосредственно через гамма-вариации) во второе. И 1 (один), являющийся первым в гамма-вариациях, означает, что они являются просто экспоненциальным распределением, выбранным как -log(U(0,1)). Пожалуйста, проверьте, если моя математика в порядке, но таким образом выборка может работать