Merge pull request #7414 from charris/tweak-corrcoef

Tweak corrcoef

2 files changed, 37 insertions, 8 deletions

diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index a2a24618b..c155babef 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -2461,11 +2461,17 @@ def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
 
     # Handles complex arrays too
     m = np.asarray(m)
+    if m.ndim > 2:
+        raise ValueError("m has more than 2 dimensions")
+
     if y is None:
         dtype = np.result_type(m, np.float64)
     else:
         y = np.asarray(y)
+        if y.ndim > 2:
+            raise ValueError("y has more than 2 dimensions")
         dtype = np.result_type(m, y, np.float64)
+
     X = array(m, ndmin=2, dtype=dtype)
     if rowvar == 0 and X.shape[0] != 1:
         X = X.T
@@ -2589,6 +2595,12 @@ def corrcoef(x, y=None, rowvar=1, bias=np._NoValue, ddof=np._NoValue):
 
     Notes
     -----
+    Due to floating point rounding the resulting array may not be Hermitian,
+    the diagonal elements may not be 1, and the elements may not satisfy the
+    inequality abs(a) <= 1. The real and imaginary parts are clipped to the
+    interval [-1,  1] in an attempt to improve on that situation but is not
+    much help in the complex case.
+
     This function accepts but discards arguments `bias` and `ddof`.  This is
     for backwards compatibility with previous versions of this function.  These
     arguments had no effect on the return values of the function and can be
@@ -2602,13 +2614,21 @@ def corrcoef(x, y=None, rowvar=1, bias=np._NoValue, ddof=np._NoValue):
     c = cov(x, y, rowvar)
     try:
         d = diag(c)
-    except ValueError:  # scalar covariance
+    except ValueError:
+        # scalar covariance
         # nan if incorrect value (nan, inf, 0), 1 otherwise
         return c / c
-    d = sqrt(d)
-    # calculate "c / multiply.outer(d, d)" row-wise ... for memory and speed
-    for i in range(0, d.size):
-        c[i,:] /= (d * d[i])
+    stddev = sqrt(d.real)
+    c /= stddev[:, None]
+    c /= stddev[None, :]
+
+    # Clip real and imaginary parts to [-1, 1].  This does not guarantee
+    # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
+    # excessive work.
+    np.clip(c.real, -1, 1, out=c.real)
+    if np.iscomplexobj(c):
+        np.clip(c.imag, -1, 1, out=c.imag)
+
     return c
 
 
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 91ead0998..e04a497c1 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -1733,8 +1733,13 @@ class TestCorrCoef(TestCase):
                             [[1., -1.], [-1.,  1.]])
 
     def test_simple(self):
-        assert_almost_equal(corrcoef(self.A), self.res1)
-        assert_almost_equal(corrcoef(self.A, self.B), self.res2)
+        tgt1 = corrcoef(self.A)
+        assert_almost_equal(tgt1, self.res1)
+        assert_(np.all(np.abs(tgt1) <= 1.0))
+
+        tgt2 = corrcoef(self.A, self.B)
+        assert_almost_equal(tgt2, self.res2)
+        assert_(np.all(np.abs(tgt2) <= 1.0))
 
     def test_ddof(self):
         # ddof raises DeprecationWarning
@@ -1760,7 +1765,10 @@ class TestCorrCoef(TestCase):
 
     def test_complex(self):
         x = np.array([[1, 2, 3], [1j, 2j, 3j]])
-        assert_allclose(corrcoef(x), np.array([[1., -1.j], [1.j, 1.]]))
+        res = corrcoef(x)
+        tgt = np.array([[1., -1.j], [1.j, 1.]])
+        assert_allclose(res, tgt)
+        assert_(np.all(np.abs(res) <= 1.0))
 
     def test_xy(self):
         x = np.array([[1, 2, 3]])
@@ -1781,6 +1789,7 @@ class TestCorrCoef(TestCase):
         with np.errstate(all='raise'):
             c = corrcoef(x)
         assert_array_almost_equal(c, np.array([[1., -1.], [-1., 1.]]))
+        assert_(np.all(np.abs(c) <= 1.0))
 
 
 class TestCov(TestCase):