np.random.dirichlet with small parameter: embed future solution in current numpy
Question: 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](https://github.com/numpy/numpy/pull/5872)
and [here](https://github.com/numpy/numpy/issues/5851) 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?


Answer: 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

