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 

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

Вопрос: Может ли кто-то предложить другой способ привлечь в питона Дирихле или указать мне решение использовать новый пробник без перекомпилировать 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() 

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

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

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

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

UPDATE

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

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

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