Repository: numpy/numpy
Pull request: BUG: Fix numpy.random.dirichlet returns NaN for small 'alpha' parameters.
<!-- Please be sure you are following the instructions in the dev guidelines
http://www.numpy.org/devdocs/dev/development_workflow.html
-->

<!-- We'd appreciate it if your commit message is properly formatted
http://www.numpy.org/devdocs/dev/development_workflow.html#writing-the-commit-message
-->
BUG: Fix numpy.random.dirichlet returns NaN for small 'alpha' parameters.

The PR introduces a second code branch into the 'dirichlet' method of 'Generator'. This branch will be executed whenever the maximum of all 'alpha' parameters for the dirichlet distribution is smaller than one. In this case the random vector will be generated by the stick breaking approach, cf. https://en.wikipedia.org/wiki/Dirichlet_process. 

When the 'alpha' smallness condition is not triggered, i.e. 'alpha.max() >= 1.' the existing code path is executed. 

Closes #5851.

Problem statement:
Bug in np.random.dirichlet for small alpha parameters
Hi,

I encountered a bug when using np.random.dirichlet with small alpha parameters. Call and traceback are below.

``` python
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-86-73c2067e20c1> in <module>()
----> 1 np.random.dirichlet([0.0001, 0.0, 0.0001])

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

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

ZeroDivisionError: float division
```

I am using numpy-1.9.1. 

I believe this is a floating point issue, the distribution has almost all of its mass very close to either (1, 0, 0) or (0, 0, 1) .

The 'float division' error already occurs for larger values<1, e.g. 0.001. 

It is likely that this occurs because of the Dirichlet distribution is usually sampled via the Gamma distribution followed by normalization. If all values returned from Gamma sampling are zero than a float division error occurs. 

In addition 

``` python
np.random.beta(0.0001, 0.0001)
```

produces 'nan' most of the time, while it should be alternating 'almost always' between (1, 0) and 
(0, 1) 

`python scipy.special.betainc(0.0001, 0.0001, 1e-50) = 0.494..`

It might not be able to fix that in the current algorithmic framework but maybe it is possible to discourage/prevent users from supplying too small parameters.

Wow, this wasn't supposed to become such a long post. Thanks to anyone reading/considering this issue.

Hints from issue discussion:
Verified for current master.

There is a vintage issue (#2089) also related to Dirichlet parameters. It seems that all that code could use a rethink...

This looks fixable with logarithms.

something like this, except in C

``` python
from __future__ import print_function
import numpy as np

def _jhonk_naive(a, b):
    while True:
        u, v = np.random.rand(2)
        x, y = np.power((u, v), (1/a, 1/b))
        if x + y <= 1:
            return x / (x + y)

def _jhonk(a, b):
    while True:
        u, v = np.random.rand(2)
        logx = np.log(u) / a
        logy = np.log(v) / b
        m = np.maximum(logx, logy)
        logxm = logx - m
        logym = logy - m
        logxpym = np.log(np.exp(logxm) + np.exp(logym))
        lse_denom = m + logxpym
        if lse_denom <= 0:
            return np.exp(logxm - logxpym)

def jhonk(a, b, n, naive=False):
    rvs = _jhonk_naive if naive else _jhonk
    return np.array([rvs(a, b) for _ in range(n)])

print(np.random.beta(0.001, 0.001, 10), '\n')
print(jhonk(0.001, 0.001, 10, naive=True), '\n')
print(jhonk(0.001, 0.001, 10, naive=False), '\n')
print(np.histogram(jhonk(0.001, 0.001, 10000, naive=False))[0], '\n')
print(jhonk(1, 1, 10, naive=False), '\n')
print(np.histogram(jhonk(1, 1, 10000, naive=False))[0], '\n')
```

```
$ python check-beta-sample.py 
[  0.  nan  nan   1.   1.  nan   0.  nan   1.  nan] 

check-beta-sample.py:9: RuntimeWarning: invalid value encountered in double_scalars
  return x / (x + y)
[  1.   1.   0.  nan   0.   1.   0.   0.  nan  nan] 

[  1.00000000e+000   1.00000000e+000   1.00000000e+000   0.00000000e+000
   1.00000000e+000   2.66868323e-085   1.00000000e+000   1.00000000e+000
   2.89780592e-157   1.00000000e+000] 

[5001    3    3    3    2    1    1    0    3 4983] 

[ 0.35706113  0.54406512  0.35814279  0.14493494  0.96739129  0.09188022
  0.13299241  0.39228852  0.4331817   0.31667118] 

[ 973 1069  977 1005  989 1027  939 1018 1010  993] 
```

**Edit**: This would be for the beta distribution.

As far as I can tell this is a very nice solution to the problem. The naive branch of the code will be much faster though. It might be good to define values a_min, b_min so that the log-space code is only executed whenever a<a_min and b<b_min. 

I think the probability for x to be smaller than some tolerance, tol is tol^a, so that the probability for x and y to be jointly smaller than some tolerance is tol^(a+b). 

For tol=1e-323 a= 0.1 and b=1e-15 the probability to find both x, y<tol is still smaller than 10^{-33}. 

A test with a=0.1, b=0.0001 for N=1e9 samples with the naive approach, did not results in any nans,

``` python
In [9]: x = np.random.beta(0.1, 0.0001, size=(1e4,1e5))

In [10]: np.any(np.isnan(x))
Out[10]: False
```

A similar test with a=0.01, b=0.001 results in nan values already for N=1e4. Based on this observation and the above probability argument it might be safe to switch to the log-space algorithm whenever max(a, b)<0.1. 

There is still a tiny chance that the code will produce nan, but I believe it is very small so that it might not be worth to sacrifice the efficiency of the naive approach by using it for all a, b<1.

Does your proposal change what values you get from a RNG using a given seed?
On May 8, 2015 3:14 AM, "Benjamin Trendelkamp-Schroer" <
notifications@github.com> wrote:

> As far as I can tell this is a very nice solution to the problem. The
> naive branch of the code will be much faster though. It might be good to
> define values a_min, b_min so that the log-space code is only executed
> whenever a<a_min and b<b_min.
> 
> I think the probability for x to be smaller than some tolerance, tol is
> tol^a, so that the probability for x and y to be jointly smaller than some
> tolerance is tol^(a+b).
> 
> For tol=1e-323 a= 0.1 and b=1e-15 the probability to find both x, y<tol is
> still smaller than 10^{-33}.
> 
> A test with a=0.1, b=0.0001 for N=1e9 samples with the naive approach, did
> not results in any nans,
> 
> In [9]: x = np.random.beta(0.1, 0.0001, size=(1e4,1e5))
> 
> In [10]: np.any(np.isnan(x))
> Out[10]: False
> 
> A similar test with a=0.01, b=0.001 results in nan values already for
> N=1e4. Based on this observation and the above probability argument it
> might be safe to switch to the log-space algorithm whenever max(a, b)<0.1.
> 
> There is still a chance that the code will produce nan, but I believe it
> is very small so that it might not be worth to sacrifice the efficiency of
> the naive approach by using it for all a, b<1.
> 
> —
> Reply to this email directly or view it on GitHub
> https://github.com/numpy/numpy/issues/5851#issuecomment-100183473.

That's soce nice math indeed!

I think it would be better to have the certainty, not a high probabilistic expectation, that NaNs are not going to be produced. Does the following have any loophole in it?

```
while(1) {
    double u = rk_double(state);
    double v = rk_double(state);
    double x = pow(u, 1.0/a);
    double y = pow(v, 1.0/b);

    if (x + y <= 1) {
        if (x + y > 0) { /* Should we perhaps check against some epsilon instead? */
            return x / (x + y);
        }
        else {
            double logx = log(u) / a;
            double logy = log(v) / b;
            double logm = logx > logy ? logx : logy;
            logx -= logm;
            logy -= logm;

            return exp(logx - log(exp(logx) + exp(logy)));
        }
    }
}
```

The only uncertain cases come if either `u` or `v` happen to be exactly zero. If infinities are properly generated and handled by `log` and `exp`, then I think it would work. The only funky case would be if _both_ `u` and `v` are zero, which probably has no practical importance, but I believe we would end up doing `-inf + inf` and therefore have a genuine `NaN` generated.

Answering Nathaniel's question, I don't think the RNG state would be changed.

@jaimefrio Your implementation looks good to me. It has the advantage of more clearly preserving the number of iterations (important for Nathaniel's concern for RNG state backwards compatibility), and of not introducing additional transcendental function overhead when there are no `nan`s.  Want to make a PR?

Do you want to send a PR and become a numpy contributor, @trendelkampschroer? I'll be happy to steal @argriffing's thunder if you don't...

I think you should go ahead with it, since you proposed an elegant solution. I can not take credit for anything more than pointing out the problem.

I am quite happy to know that there is such a nice fix, I am not sure if the same strategy could be applied to the Dirichlet case, though. The current approach of generating via the Gamma distribution will definitely suffer from the same problem (zero denominator in the normalization).

Commit message: Add tests demonstrating slowness for beta and dirichlet generators for small alpha (and beta) values

Patch to review:
diff --git a/numpy/random/_generator.pyx b/numpy/random/_generator.pyx
index b6c222cc0720..584bb42d8338 100644
--- a/numpy/random/_generator.pyx
+++ b/numpy/random/_generator.pyx
@@ -4124,8 +4124,9 @@ cdef class Generator:
         # return val
 
         cdef np.npy_intp k, totsize, i, j
-        cdef np.ndarray alpha_arr, val_arr
+        cdef np.ndarray alpha_arr, val_arr, alpha_csum_arr
         cdef double *alpha_data
+        cdef double *alpha_csum_data
         cdef double *val_data
         cdef double acc, invacc
 
@@ -4136,6 +4137,11 @@ cdef class Generator:
             raise ValueError('alpha <= 0')
         alpha_data = <double*>np.PyArray_DATA(alpha_arr)
 
+        alpha_csum = np.ascontiguousarray(np.cumsum(alpha[::-1])[::-1])
+        alpha_csum_arr = <np.ndarray>np.PyArray_FROM_OTF(
+            alpha_csum, np.NPY_DOUBLE, np.NPY_ALIGNED | np.NPY_ARRAY_C_CONTIGUOUS)
+        alpha_csum_data = <double*>np.PyArray_DATA(alpha_csum_arr)
+
         if size is None:
             shape = (k,)
         else:
@@ -4150,17 +4156,33 @@ cdef class Generator:
 
         i = 0
         totsize = np.PyArray_SIZE(val_arr)
-        with self.lock, nogil:
-            while i < totsize:
-                acc = 0.0
-                for j in range(k):
-                    val_data[i+j] = random_standard_gamma(&self._bitgen,
-                                                              alpha_data[j])
-                    acc = acc + val_data[i + j]
-                invacc = 1/acc
-                for j in range(k):
-                    val_data[i + j] = val_data[i + j] * invacc
-                i = i + k
+
+        # Small alpha case: Use stick-breaking approach with beta random variates
+        if alpha_arr.max() < 1.0:
+            with self.lock, nogil:
+                while i < totsize:
+                    acc = 1.
+                    for j in range(k-1):
+                        val_data[i+j] = acc * random_beta(&self._bitgen,
+                                                          alpha_data[j],
+                                                          alpha_csum_data[j+1])
+                        acc *= (1. - val_data[i+j])
+                    val_data[i+j+1] = acc
+                    i = i + k
+
+        # Non-small alpha case: Perform unit normalisation of a sum of gamma random variates
+        else:
+            with self.lock, nogil:
+                while i < totsize:
+                    acc = 0.0
+                    for j in range(k):
+                        val_data[i+j] = random_standard_gamma(&self._bitgen,
+                                                                  alpha_data[j])
+                        acc = acc + val_data[i + j]
+                    invacc = 1/acc
+                    for j in range(k):
+                        val_data[i + j] = val_data[i + j] * invacc
+                    i = i + k
 
         return diric
 
diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py
index d835f16bd7b2..3bf96ff89c8f 100644
--- a/numpy/random/tests/test_generator_mt19937.py
+++ b/numpy/random/tests/test_generator_mt19937.py
@@ -987,6 +987,15 @@ def test_beta(self):
              [2.830254190078299e-04, 1.744709918330393e-01]])
         assert_array_almost_equal(actual, desired, decimal=15)
 
+    def test_beta_small_a_and_b(self):
+        eps = 1.0e-12  # 1.0e-11 -> runtime x50; 1.0e-12 -> runtime x200
+        a = eps
+        b = eps * 1.0e-3
+
+        random = Generator(MT19937(self.seed))
+        actual = random.beta(a, b)
+        assert actual == 1.
+
     def test_binomial(self):
         random = Generator(MT19937(self.seed))
         actual = random.binomial(100.123, .456, size=(3, 2))
@@ -1054,6 +1063,29 @@ def test_dirichlet_alpha_non_contiguous(self):
                                   size=(3, 2))
         assert_array_almost_equal(non_contig, contig)
 
+    def test_dirichlet_small_alpha(self):
+        eps = 1.0e-9  # 1.0e-10 -> runtime x 10; 1e-11 -> runtime x 200, etc.
+        alpha = eps * np.array([1., 1.0e-3])
+        random = Generator(MT19937(self.seed))
+        actual = random.dirichlet(alpha, size=(3, 2))
+        expected = np.array([
+            [
+                [1., 0.],
+                [1., 0.]
+            ],
+
+            [
+                [1., 0.],
+                [1., 0.]
+            ],
+
+            [
+                [1., 0.],
+                [1., 0.]
+            ]
+        ])
+        assert_array_almost_equal(actual, expected, decimal=15)
+
     def test_exponential(self):
         random = Generator(MT19937(self.seed))
         actual = random.exponential(1.1234, size=(3, 2))


Code review comments:
diff_hunk: @@ -4150,17 +4156,33 @@ cdef class Generator:
 
         i = 0
         totsize = np.PyArray_SIZE(val_arr)
-        with self.lock, nogil:
-            while i < totsize:
-                acc = 0.0
-                for j in range(k):
-                    val_data[i+j] = random_standard_gamma(&self._bitgen,
-                                                              alpha_data[j])
-                    acc = acc + val_data[i + j]
-                invacc = 1/acc
-                for j in range(k):
-                    val_data[i + j] = val_data[i + j] * invacc
-                i = i + k
+
+        # Small alpha case: Use stick-breaking approach with beta random variates
+        if alpha_arr.max() < 1.0:
+            with self.lock, nogil:
+                while i < totsize:
+                    acc = 1.
+                    for j in range(k-1):
+                        val_data[i+j] = acc * random_beta(&self._bitgen,
+                                                          alpha_data[j],
+                                                          alpha_csum_data[j+1])
+                        acc *= (1. - val_data[i+j])
path: numpy/random/_generator.pyx
text: @user1:
It looks like you want `acc` to "accumulate" the products of 1 minus the beta random variates, but you can't use `val_data[i+j]` for that, because it includes the product of `acc` and the random beta (not just the random beta).  You can make this work if you save the random beta variate in its own variable.  The body of the `while` loop then looks like this:

                    acc = 1.
                    for j in range(k - 1):
                        v = random_beta(&self._bitgen,
                                        alpha_data[j], alpha_csum_data[j + 1])
                        val_data[i + j] = acc * v
                        acc *= (1. - v)
                    val_data[i + j + 1] = acc
                    i = i + k

where `v` is `cdef`'d to be type `double` at the beginning of the function.

@author:
Thanks for the hint. I have pushed the suggested changes. This should restore the correct distribution.

@user1:
Thanks, the update looks good.

I'm currently running a bunch of tests to see if anything breaks. :)