There is an ongoing discussion about the current np.random.dirichlet function, as it does not function for small parameters:
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
The discussion can be read here and here and points out, that it is a normalization error. Currently the proposed enhancement to switch samplers for small parameters cannot be merged into master of numpy for several reasons.
Question: Can someone suggest a different way to draw dirichlets in python or point me to a solution to use the new sampler without recompile my numpy and/or working on an unreleased branch?
Ok, lets try the following. Here is Beta(alpha,beta) variate sampling which should work for any small numbers.
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()
Using it, you might try to sample Dirichlet via Beta variate as described in the wikipedia
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))
If it works, it could be optimized to have all partial sums precomputed.
UPDATE
Sampling above relies on the fact, that Dirichlet could be sampled via beta variate, and that is better (but slower) choice if case of small parameters. In turn, beta variate could be sampled as pair of gamma variates:
beta(a, b) = gamma(1, a) / (gamma(1, a) + gamma(1, b))
So small parameters are moved from being first in gamma (if you sample Dirichlet directly via gamma variates) to being second. And 1 (one) being first in gamma-variates means they are just exponential distribution, sampled as -log(U(0,1)). Please check if my math is ok, but that way sampling might work