Multivariate floating point precision fix (#2470)

* fixed multivariate

* multivariate fix for pvals > 1

* updated multivariate with hack

* bugfix

* float64 patch to multivariate's multinomial

* reinstated p.squeeze(), accidentally took it out

* updated based on failed test

* further minor fixes to multinomial

* Force push to restart Travis

Author: ericmjl

pymc3/distributions/multivariate.py

@@ -503,18 +503,22 @@ def __init__(self, n, p, *args, **kwargs):
         self.mode = tt.cast(tround(self.mean), 'int32')
 
     def _random(self, n, p, size=None):
+        original_dtype = p.dtype
+        # Set float type to float64 for numpy. This change is related to numpy issue #8317 (https://github.com/numpy/numpy/issues/8317)
+        p = p.astype('float64')
+        # Now, re-normalize all of the values in float64 precision. This is done inside the conditionals
         if size == p.shape:
             size = None
-
         if p.ndim == 1:
+            p = p / p.sum()
             randnum = np.random.multinomial(n, p.squeeze(), size=size)
         elif p.ndim == 2:
+            p = p / p.sum(axis=1, keepdims=True)
             randnum = np.asarray([np.random.multinomial(n, pp, size=size) for pp in p])
         else:
             raise ValueError('Outcome probabilities must be 1- or 2-dimensional '
                              '(supplied `p` has {} dimensions)'.format(p.ndim))
-
-        return randnum
+        return randnum.astype(original_dtype)
 
     def random(self, point=None, size=None):
         n, p = draw_values([self.n, self.p], point=point)