Merge pull request #75 from bashtage/update-to-latest-numpy

SYNC: Sync with upstream changes in NumPy

Author: bashtage

randomstate/distributions.c

@@ -501,7 +501,7 @@ double random_beta(aug_state *state, double a, double b)
     if ((a <= 1.0) && (b <= 1.0))
     {
         double U, V, X, Y;
-        /* Use Jonk's algorithm */
+        /* Use Johnk's algorithm */
 
         while (1)
         {

randomstate/randomstate.pyx

@@ -3934,9 +3934,10 @@ cdef class RandomState:
                  0.0, '', CONS_NONE)
 
     # Multivariate distributions:
-    def multivariate_normal(self, mean, cov, size=None, method=__normal_method):
+    def multivariate_normal(self, mean, cov, size=None, check_valid='warn',
+                            tol=1e-8, method=__normal_method):
         """
-        multivariate_normal(mean, cov, size=None, method='bm')
+        multivariate_normal(mean, cov[, size, check_valid, tol])
 
         Draw random samples from a multivariate normal distribution.
 
@@ -3959,6 +3960,10 @@ cdef class RandomState:
             generated, and packed in an `m`-by-`n`-by-`k` arrangement.  Because
             each sample is `N`-dimensional, the output shape is ``(m,n,k,N)``.
             If no shape is specified, a single (`N`-D) sample is returned.
+        check_valid : { 'warn', 'raise', 'ignore' }, optional
+            Behavior when the covariance matrix is not positive semidefinite.
+        tol : float, optional
+            Tolerance when checking the singular values in covariance matrix.
         method : str, optional
             Either 'bm' or 'zig'. 'bm' uses the default Box-Muller transformations
             method.  'zig' uses the much faster Ziggurat method of Marsaglia and Tsang.
@@ -4073,12 +4078,20 @@ cdef class RandomState:
         # not zero. We continue to use the SVD rather than Cholesky in
         # order to preserve current outputs. Note that symmetry has not
         # been checked.
+
         (u, s, v) = svd(cov)
-        neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)
-        if np.any(neg):
-            s[neg] = 0.
-            warnings.warn("covariance is not positive-semidefinite.",
-                          RuntimeWarning)
+
+        if check_valid != 'ignore':
+            if check_valid != 'warn' and check_valid != 'raise':
+                raise ValueError("check_valid must equal 'warn', 'raise', or 'ignore'")
+
+            psd = np.allclose(np.dot(v.T * s, v), cov, rtol=tol, atol=tol)
+            if not psd:
+                if check_valid == 'warn':
+                    warnings.warn("covariance is not positive-semidefinite.",
+                                  RuntimeWarning)
+                else:
+                    raise ValueError("covariance is not positive-semidefinite.")
 
         x = np.dot(x, np.sqrt(s)[:, None] * v)
         x += mean

randomstate/tests/tests_numpy_mt19937.py

@@ -3,7 +3,8 @@
 import numpy as np
 from numpy.testing import (
         TestCase, run_module_suite, assert_, assert_raises, assert_equal,
-        assert_warns, assert_array_equal, assert_array_almost_equal)
+        assert_warns, assert_no_warnings, assert_array_equal,
+        assert_array_almost_equal)
 from numpy.compat import asbytes
 import sys
 import warnings
@@ -614,28 +615,37 @@ def test_multinomial(self):
     def test_multivariate_normal(self):
         mt19937.seed(self.seed)
         mean = (.123456789, 10)
-        # Hmm... not even symmetric.
-        cov = [[1, 0], [1, 0]]
+        cov = [[1, 0], [0, 1]]
         size = (3, 2)
         actual = mt19937.multivariate_normal(mean, cov, size)
-        desired = np.array([[[-1.47027513018564449, 10.],
-                             [-1.65915081534845532, 10.]],
-                            [[-2.29186329304599745, 10.],
-                             [-1.77505606019580053, 10.]],
-                            [[-0.54970369430044119, 10.],
-                             [0.29768848031692957, 10.]]])
+        desired = np.array([[[1.463620246718631, 11.73759122771936 ],
+                             [1.622445133300628, 9.771356667546383]],
+                            [[2.154490787682787, 12.170324946056553],
+                             [1.719909438201865, 9.230548443648306]],
+                            [[0.689515026297799, 9.880729819607714],
+                             [-0.023054015651998, 9.201096623542879]]])
+
         assert_array_almost_equal(actual, desired, decimal=15)
 
         # Check for default size, was raising deprecation warning
         actual = mt19937.multivariate_normal(mean, cov)
-        desired = np.array([-0.79441224511977482, 10.])
+        desired = np.array([0.895289569463708, 9.17180864067987])
         assert_array_almost_equal(actual, desired, decimal=15)
 
-        # Check that non positive-semidefinite covariance raises warning
+        # Check that non positive-semidefinite covariance warns with
+        # RuntimeWarning
         mean = [0, 0]
-        cov = [[1, 1 + 1e-10], [1 + 1e-10, 1]]
+        cov = [[1, 2], [2, 1]]
         assert_warns(RuntimeWarning, mt19937.multivariate_normal, mean, cov)
 
+        # and that it doesn't warn with RuntimeWarning check_valid='ignore'
+        assert_no_warnings(mt19937.multivariate_normal, mean, cov,
+                           check_valid='ignore')
+
+        # and that it raises with RuntimeWarning check_valid='raises'
+        assert_raises(ValueError, mt19937.multivariate_normal, mean, cov,
+                      check_valid='raise')
+
     def test_negative_binomial(self):
         mt19937.seed(self.seed)
         actual = mt19937.negative_binomial(n=100, p=.12345, size=(3, 2))
@@ -810,7 +820,9 @@ def test_uniform_range_bounds(self):
         assert_raises(OverflowError, func, [0], [np.inf])
 
         # (fmax / 1e17) - fmin is within range, so this should not throw
-        mt19937.uniform(low=fmin, high=fmax / 1e17)
+        # account for i386 extended precision DBL_MAX / 1e17 + DBL_MAX >
+        # DBL_MAX by increasing fmin a bit
+        mt19937.uniform(low=np.nextafter(fmin, 1), high=fmax / 1e17)
 
     def test_vonmises(self):
         mt19937.seed(self.seed)