Merge pull request #3534 from charris/nan-stat-functions

Add nanmean, nanvar, and nanstd functions.

11 files changed, 1507 insertions, 524 deletions

diff --git a/numpy/core/_methods.py b/numpy/core/_methods.py
index 66a05e868..27e445a54 100644
--- a/numpy/core/_methods.py
+++ b/numpy/core/_methods.py
@@ -5,9 +5,12 @@ and the Python code for the NumPy-namespace function
 """
 from __future__ import division, absolute_import, print_function
 
+import warnings
+
 from numpy.core import multiarray as mu
 from numpy.core import umath as um
 from numpy.core.numeric import asanyarray
+from numpy.core import numerictypes as nt
 
 def _amax(a, axis=None, out=None, keepdims=False):
     return um.maximum.reduce(a, axis=axis,
@@ -46,58 +49,65 @@ def _count_reduce_items(arr, axis):
 def _mean(a, axis=None, dtype=None, out=None, keepdims=False):
     arr = asanyarray(a)
 
-    # Upgrade bool, unsigned int, and int to float64
-    if dtype is None and arr.dtype.kind in ['b','u','i']:
-        ret = um.add.reduce(arr, axis=axis, dtype='f8',
-                            out=out, keepdims=keepdims)
-    else:
-        ret = um.add.reduce(arr, axis=axis, dtype=dtype,
-                            out=out, keepdims=keepdims)
     rcount = _count_reduce_items(arr, axis)
+    # Make this warning show up first
+    if rcount == 0:
+        warnings.warn("Mean of empty slice.", RuntimeWarning)
+
+
+    # Cast bool, unsigned int, and int to float64 by default
+    if dtype is None and issubclass(arr.dtype.type, (nt.integer, nt.bool_)):
+        dtype = mu.dtype('f8')
+
+    ret = um.add.reduce(arr, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
     if isinstance(ret, mu.ndarray):
-        ret = um.true_divide(ret, rcount,
-                        out=ret, casting='unsafe', subok=False)
+        ret = um.true_divide(
+                ret, rcount, out=ret, casting='unsafe', subok=False)
     else:
-        ret = ret / float(rcount)
+        ret = ret.dtype.type(ret / rcount)
+
     return ret
 
-def _var(a, axis=None, dtype=None, out=None, ddof=0,
-                            keepdims=False):
+def _var(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     arr = asanyarray(a)
 
-    # First compute the mean, saving 'rcount' for reuse later
-    if dtype is None and arr.dtype.kind in ['b','u','i']:
-        arrmean = um.add.reduce(arr, axis=axis, dtype='f8', keepdims=True)
-    else:
-        arrmean = um.add.reduce(arr, axis=axis, dtype=dtype, keepdims=True)
     rcount = _count_reduce_items(arr, axis)
+    # Make this warning show up on top.
+    if ddof >= rcount:
+        warnings.warn("Degrees of freedom <= 0 for slice", RuntimeWarning)
+
+    # Cast bool, unsigned int, and int to float64 by default
+    if dtype is None and issubclass(arr.dtype.type, (nt.integer, nt.bool_)):
+        dtype = mu.dtype('f8')
+
+    # Compute the mean.
+    # Note that if dtype is not of inexact type then arraymean will
+    # not be either.
+    arrmean = um.add.reduce(arr, axis=axis, dtype=dtype, keepdims=True)
     if isinstance(arrmean, mu.ndarray):
-        arrmean = um.true_divide(arrmean, rcount,
-                            out=arrmean, casting='unsafe', subok=False)
+        arrmean = um.true_divide(
+                arrmean, rcount, out=arrmean, casting='unsafe', subok=False)
     else:
-        arrmean = arrmean / float(rcount)
+        arrmean = arrmean.dtype.type(arrmean / rcount)
 
-    # arr - arrmean
+    # Compute sum of squared deviations from mean
+    # Note that x may not be inexact
     x = arr - arrmean
-
-    # (arr - arrmean) ** 2
-    if arr.dtype.kind == 'c':
+    if issubclass(arr.dtype.type, nt.complexfloating):
         x = um.multiply(x, um.conjugate(x), out=x).real
     else:
         x = um.multiply(x, x, out=x)
-
-    # add.reduce((arr - arrmean) ** 2, axis)
     ret = um.add.reduce(x, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
 
-    # add.reduce((arr - arrmean) ** 2, axis) / (n - ddof)
-    if not keepdims and isinstance(rcount, mu.ndarray):
-        rcount = rcount.squeeze(axis=axis)
-    rcount -= ddof
+    # Compute degrees of freedom and make sure it is not negative.
+    rcount = max([rcount - ddof, 0])
+
+    # divide by degrees of freedom
     if isinstance(ret, mu.ndarray):
-        ret = um.true_divide(ret, rcount,
-                        out=ret, casting='unsafe', subok=False)
+        ret = um.true_divide(
+                ret, rcount, out=ret, casting='unsafe', subok=False)
     else:
-        ret = ret / float(rcount)
+        ret = ret.dtype.type(ret / rcount)
 
     return ret
 
@@ -108,6 +118,7 @@ def _std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     if isinstance(ret, mu.ndarray):
         ret = um.sqrt(ret, out=ret)
     else:
-        ret = um.sqrt(ret)
+        ret = ret.dtype.type(um.sqrt(ret))
 
     return ret
+
diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py
index 9d7c226fc..cb1c4ff05 100644
--- a/numpy/core/fromnumeric.py
+++ b/numpy/core/fromnumeric.py
@@ -14,17 +14,16 @@ from . import _methods
 _dt_ = nt.sctype2char
 
 
-# functions that are now methods
-__all__ = ['take', 'reshape', 'choose', 'repeat', 'put',
-           'swapaxes', 'transpose', 'sort', 'argsort', 'partition', 'argpartition',
-           'argmax', 'argmin',
-           'searchsorted', 'alen',
-           'resize', 'diagonal', 'trace', 'ravel', 'nonzero', 'shape',
-           'compress', 'clip', 'sum', 'product', 'prod', 'sometrue', 'alltrue',
-           'any', 'all', 'cumsum', 'cumproduct', 'cumprod', 'ptp', 'ndim',
-           'rank', 'size', 'around', 'round_', 'mean', 'std', 'var', 'squeeze',
-           'amax', 'amin',
-          ]
+# functions that are methods
+__all__ = [
+        'alen', 'all', 'alltrue', 'amax', 'amin', 'any', 'argmax',
+        'argmin', 'argpartition', 'argsort', 'around', 'choose', 'clip',
+        'compress', 'cumprod', 'cumproduct', 'cumsum', 'diagonal', 'mean',
+        'ndim', 'nonzero', 'partition', 'prod', 'product', 'ptp', 'put',
+        'rank', 'ravel', 'repeat', 'reshape', 'resize', 'round_',
+        'searchsorted', 'shape', 'size', 'sometrue', 'sort', 'squeeze',
+        'std', 'sum', 'swapaxes', 'take', 'trace', 'transpose', 'var',
+        ]
 
 
 try:
@@ -63,9 +62,9 @@ def take(a, indices, axis=None, out=None, mode='raise'):
         The source array.
     indices : array_like
         The indices of the values to extract.
-        
+
         .. versionadded:: 1.8.0
-        
+
         Also allow scalars for indices.
     axis : int, optional
         The axis over which to select values. By default, the flattened
@@ -2665,6 +2664,7 @@ def mean(a, axis=None, dtype=None, out=None, keepdims=False):
     See Also
     --------
     average : Weighted average
+    std, var, nanmean, nanstd, nanvar
 
     Notes
     -----
@@ -2711,7 +2711,6 @@ def mean(a, axis=None, dtype=None, out=None, keepdims=False):
     return _methods._mean(a, axis=axis, dtype=dtype,
                             out=out, keepdims=keepdims)
 
-
 def std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     """
     Compute the standard deviation along the specified axis.
@@ -2752,7 +2751,7 @@ def std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
 
     See Also
     --------
-    var, mean
+    var, mean, nanmean, nanstd, nanvar
     numpy.doc.ufuncs : Section "Output arguments"
 
     Notes
@@ -2855,8 +2854,7 @@ def var(a, axis=None, dtype=None, out=None, ddof=0,
 
     See Also
     --------
-    std : Standard deviation
-    mean : Average
+    std , mean, nanmean, nanstd, nanvar
     numpy.doc.ufuncs : Section "Output arguments"
 
     Notes
@@ -2915,3 +2913,4 @@ def var(a, axis=None, dtype=None, out=None, ddof=0,
 
     return _methods._var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
                                 keepdims=keepdims)
+
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index b7235d05f..d37bcd331 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -2399,7 +2399,178 @@ class TestView(TestCase):
         assert_array_equal(y, z)
         assert_array_equal(y, [67305985, 134678021])
 
+def _mean(a, **args):
+    return a.mean(**args)
+
+def _var(a, **args):
+    return a.var(**args)
+
+def _std(a, **args):
+    return a.std(**args)
+
 class TestStats(TestCase):
+
+    funcs = [_mean, _var, _std]
+
+    def setUp(self):
+        np.random.seed(range(3))
+        self.rmat = np.random.random((4, 5))
+        self.cmat = self.rmat + 1j * self.rmat
+
+    def test_keepdims(self):
+        mat = np.eye(3)
+        for f in self.funcs:
+            for axis in [0, 1]:
+                res = f(mat, axis=axis, keepdims=True)
+                assert_(res.ndim == mat.ndim)
+                assert_(res.shape[axis] == 1)
+            for axis in [None]:
+                res = f(mat, axis=axis, keepdims=True)
+                assert_(res.shape == (1, 1))
+
+    def test_out(self):
+        mat = np.eye(3)
+        for f in self.funcs:
+            out = np.zeros(3)
+            tgt = f(mat, axis=1)
+            res = f(mat, axis=1, out=out)
+            assert_almost_equal(res, out)
+            assert_almost_equal(res, tgt)
+        out = np.empty(2)
+        assert_raises(ValueError, f, mat, axis=1, out=out)
+        out = np.empty((2, 2))
+        assert_raises(ValueError, f, mat, axis=1, out=out)
+
+    def test_dtype_from_input(self):
+        icodes = np.typecodes['AllInteger']
+        fcodes = np.typecodes['AllFloat']
+
+        # integer types
+        for f in self.funcs:
+            for c in icodes:
+                mat = np.eye(3, dtype=c)
+                tgt = np.float64
+                res = f(mat, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                res = f(mat, axis=None).dtype.type
+                assert_(res is tgt)
+        # mean for float types
+        for f in [_mean]:
+            for c in fcodes:
+                mat = np.eye(3, dtype=c)
+                tgt = mat.dtype.type
+                res = f(mat, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                res = f(mat, axis=None).dtype.type
+                assert_(res is tgt)
+        # var, std for float types
+        for f in [_var, _std]:
+            for c in fcodes:
+                mat = np.eye(3, dtype=c)
+                tgt = mat.real.dtype.type
+                res = f(mat, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                res = f(mat, axis=None).dtype.type
+                assert_(res is tgt)
+
+    def test_dtype_from_dtype(self):
+        icodes = np.typecodes['AllInteger']
+        fcodes = np.typecodes['AllFloat']
+        mat = np.eye(3)
+
+        # stats for integer types
+        # fixme:
+        # this needs definition as there are lots places along the line
+        # where type casting may take place.
+        #for f in self.funcs:
+            #for c in icodes:
+                #tgt = np.dtype(c).type
+                #res = f(mat, axis=1, dtype=c).dtype.type
+                #assert_(res is tgt)
+                ## scalar case
+                #res = f(mat, axis=None, dtype=c).dtype.type
+                #assert_(res is tgt)
+
+        # stats for float types
+        for f in self.funcs:
+            for c in fcodes:
+                tgt = np.dtype(c).type
+                res = f(mat, axis=1, dtype=c).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                res = f(mat, axis=None, dtype=c).dtype.type
+                assert_(res is tgt)
+
+    def test_ddof(self):
+        for f in [_var]:
+            for ddof in range(3):
+                dim = self.rmat.shape[1]
+                tgt = f(self.rmat, axis=1) * dim
+                res = f(self.rmat, axis=1, ddof=ddof) * (dim - ddof)
+        for f in [_std]:
+            for ddof in range(3):
+                dim = self.rmat.shape[1]
+                tgt = f(self.rmat, axis=1) * np.sqrt(dim)
+                res = f(self.rmat, axis=1, ddof=ddof) * np.sqrt(dim - ddof)
+                assert_almost_equal(res, tgt)
+                assert_almost_equal(res, tgt)
+
+    def test_ddof_too_big(self):
+        dim = self.rmat.shape[1]
+        for f in [_var, _std]:
+            for ddof in range(dim, dim + 2):
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    res = f(self.rmat, axis=1, ddof=ddof)
+                    assert_(not (res < 0).any())
+                    assert_(len(w) > 0)
+                    assert_(issubclass(w[0].category, RuntimeWarning))
+
+    def test_empty(self):
+        A = np.zeros((0,3))
+        for f in self.funcs:
+            for axis in [0, None]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_(np.isnan(f(A, axis=axis)).all())
+                    assert_(len(w) > 0)
+                    assert_(issubclass(w[0].category, RuntimeWarning))
+            for axis in [1]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_equal(f(A, axis=axis), np.zeros([]))
+
+    def test_mean_values(self):
+        for mat in [self.rmat, self.cmat]:
+            for axis in [0, 1]:
+                tgt = mat.sum(axis=axis)
+                res = _mean(mat, axis=axis) * mat.shape[axis]
+                assert_almost_equal(res, tgt)
+            for axis in [None]:
+                tgt = mat.sum(axis=axis)
+                res = _mean(mat, axis=axis) * np.prod(mat.shape)
+                assert_almost_equal(res, tgt)
+
+    def test_var_values(self):
+        for mat in [self.rmat, self.cmat]:
+            for axis in [0, 1, None]:
+                msqr = _mean(mat * mat.conj(), axis=axis)
+                mean = _mean(mat, axis=axis)
+                tgt = msqr - mean * mean.conj()
+                res = _var(mat, axis=axis)
+                assert_almost_equal(res, tgt)
+
+    def test_std_values(self):
+        for mat in [self.rmat, self.cmat]:
+            for axis in [0, 1, None]:
+                tgt = np.sqrt(_var(mat, axis=axis))
+                res = _std(mat, axis=axis)
+                assert_almost_equal(res, tgt)
+
+
     def test_subclass(self):
         class TestArray(np.ndarray):
             def __new__(cls, data, info):
@@ -2409,6 +2580,7 @@ class TestStats(TestCase):
                 return result
             def __array_finalize__(self, obj):
                 self.info = getattr(obj, "info", '')
+
         dat = TestArray([[1,2,3,4],[5,6,7,8]], 'jubba')
         res = dat.mean(1)
         assert_(res.info == dat.info)
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 782ddd687..341884f04 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -3,6 +3,7 @@ from __future__ import division, absolute_import, print_function
 import sys
 import platform
 from decimal import Decimal
+import warnings
 
 import numpy as np
 from numpy.core import *
@@ -177,18 +178,33 @@ class TestNonarrayArgs(TestCase):
         assert_(all(mean(A,0) == array([2.5,3.5,4.5])))
         assert_(all(mean(A,1) == array([2.,5.])))
 
+        with warnings.catch_warnings(record=True) as w:
+            warnings.filterwarnings('always', '', RuntimeWarning)
+            assert_(isnan(mean([])))
+            assert_(w[0].category is RuntimeWarning)
+
     def test_std(self):
         A = [[1,2,3],[4,5,6]]
         assert_almost_equal(std(A), 1.707825127659933)
         assert_almost_equal(std(A,0), array([1.5, 1.5, 1.5]))
         assert_almost_equal(std(A,1), array([0.81649658, 0.81649658]))
 
+        with warnings.catch_warnings(record=True) as w:
+            warnings.filterwarnings('always', '', RuntimeWarning)
+            assert_(isnan(std([])))
+            assert_(w[0].category is RuntimeWarning)
+
     def test_var(self):
         A = [[1,2,3],[4,5,6]]
         assert_almost_equal(var(A), 2.9166666666666665)
         assert_almost_equal(var(A,0), array([2.25, 2.25, 2.25]))
         assert_almost_equal(var(A,1), array([0.66666667, 0.66666667]))
 
+        with warnings.catch_warnings(record=True) as w:
+            warnings.filterwarnings('always', '', RuntimeWarning)
+            assert_(isnan(var([])))
+            assert_(w[0].category is RuntimeWarning)
+
 
 class TestBoolScalar(TestCase):
     def test_logical(self):
@@ -1447,7 +1463,6 @@ class TestIsclose(object):
         assert_array_equal(x, array([inf, 1]))
         assert_array_equal(y, array([0, inf]))
 
-
 class TestStdVar(TestCase):
     def setUp(self):
         self.A = array([1,-1,1,-1])
@@ -1469,7 +1484,6 @@ class TestStdVar(TestCase):
         assert_almost_equal(std(self.A,ddof=2)**2,
                             self.real_var*len(self.A)/float(len(self.A)-2))
 
-
 class TestStdVarComplex(TestCase):
     def test_basic(self):
         A = array([1,1.j,-1,-1.j])
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index 740957215..86b0e48a2 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -1449,14 +1449,17 @@ class TestRegression(TestCase):
                       and not issubclass(x, np.timedelta64))]
         a = np.array([], dtypes[0])
         failures = []
-        for x in dtypes:
-            b = a.astype(x)
-            for y in dtypes:
-                c = a.astype(y)
-                try:
-                    np.dot(b, c)
-                except TypeError as e:
-                    failures.append((x, y))
+        # ignore complex warnings
+        with warnings.catch_warnings():
+            warnings.simplefilter('ignore', np.ComplexWarning)
+            for x in dtypes:
+                b = a.astype(x)
+                for y in dtypes:
+                    c = a.astype(y)
+                    try:
+                        np.dot(b, c)
+                    except TypeError as e:
+                        failures.append((x, y))
         if failures:
             raise AssertionError("Failures: %r" % failures)
 
diff --git a/numpy/lib/__init__.py b/numpy/lib/__init__.py
index 12acae95b..64a8550c6 100644
--- a/numpy/lib/__init__.py
+++ b/numpy/lib/__init__.py
@@ -1,11 +1,14 @@
 from __future__ import division, absolute_import, print_function
 
+import math
+
 from .info import __doc__
 from numpy.version import version as __version__
 
 from .type_check import *
 from .index_tricks import *
 from .function_base import *
+from .nanfunctions import *
 from .shape_base import *
 from .stride_tricks import *
 from .twodim_base import *
@@ -18,7 +21,6 @@ from .utils import *
 from .arraysetops import *
 from .npyio import *
 from .financial import *
-import math
 from .arrayterator import *
 from .arraypad import *
 
@@ -36,6 +38,7 @@ __all__ += utils.__all__
 __all__ += arraysetops.__all__
 __all__ += npyio.__all__
 __all__ += financial.__all__
+__all__ += nanfunctions.__all__
 
 from numpy.testing import Tester
 test = Tester().test
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index b3b5ef735..9336579af 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -1,15 +1,14 @@
 from __future__ import division, absolute_import, print_function
 
 __docformat__ = "restructuredtext en"
-__all__ = ['select', 'piecewise', 'trim_zeros', 'copy', 'iterable',
-           'percentile', 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex',
-           'disp', 'extract', 'place', 'nansum', 'nanmax', 'nanargmax',
-           'nanargmin', 'nanmin', 'vectorize', 'asarray_chkfinite', 'average',
-           'histogram', 'histogramdd', 'bincount', 'digitize', 'cov',
-           'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning',
-           'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc',
-           'add_docstring', 'meshgrid', 'delete', 'insert', 'append', 'interp',
-           'add_newdoc_ufunc']
+__all__ = [
+    'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
+    'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp',
+    'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
+    'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef',
+    'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
+    'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
+    'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc']
 
 import warnings
 import types
@@ -1361,331 +1360,6 @@ def place(arr, mask, vals):
     """
     return _insert(arr, mask, vals)
 
-def _nanop(op, fill, a, axis=None):
-    """
-    General operation on arrays with not-a-number values.
-
-    Parameters
-    ----------
-    op : callable
-        Operation to perform.
-    fill : float
-        NaN values are set to fill before doing the operation.
-    a : array-like
-        Input array.
-    axis : {int, None}, optional
-        Axis along which the operation is computed.
-        By default the input is flattened.
-
-    Returns
-    -------
-    y : {ndarray, scalar}
-        Processed data.
-
-    """
-    y = array(a, subok=True)
-
-    # We only need to take care of NaN's in floating point arrays
-    dt = y.dtype
-    if np.issubdtype(dt, np.integer) or np.issubdtype(dt, np.bool_):
-        return op(y, axis=axis)
-
-    mask = isnan(a)
-    # y[mask] = fill
-    # We can't use fancy indexing here as it'll mess w/ MaskedArrays
-    # Instead, let's fill the array directly...
-    np.copyto(y, fill, where=mask)
-    res = op(y, axis=axis)
-    mask_all_along_axis = mask.all(axis=axis)
-
-    # Along some axes, only nan's were encountered.  As such, any values
-    # calculated along that axis should be set to nan.
-    if mask_all_along_axis.any():
-        if np.isscalar(res):
-            res = np.nan
-        else:
-            res[mask_all_along_axis] = np.nan
-
-    return res
-
-def nansum(a, axis=None):
-    """
-    Return the sum of array elements over a given axis treating
-    Not a Numbers (NaNs) as zero.
-
-    Parameters
-    ----------
-    a : array_like
-        Array containing numbers whose sum is desired. If `a` is not an
-        array, a conversion is attempted.
-    axis : int, optional
-        Axis along which the sum is computed. The default is to compute
-        the sum of the flattened array.
-
-    Returns
-    -------
-    y : ndarray
-        An array with the same shape as a, with the specified axis removed.
-        If a is a 0-d array, or if axis is None, a scalar is returned with
-        the same dtype as `a`.
-
-    See Also
-    --------
-    numpy.sum : Sum across array including Not a Numbers.
-    isnan : Shows which elements are Not a Number (NaN).
-    isfinite: Shows which elements are not: Not a Number, positive and
-             negative infinity
-
-    Notes
-    -----
-    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
-    (IEEE 754). This means that Not a Number is not equivalent to infinity.
-    If positive or negative infinity are present the result is positive or
-    negative infinity. But if both positive and negative infinity are present,
-    the result is Not A Number (NaN).
-
-    Arithmetic is modular when using integer types (all elements of `a` must
-    be finite i.e. no elements that are NaNs, positive infinity and negative
-    infinity because NaNs are floating point types), and no error is raised
-    on overflow.
-
-
-    Examples
-    --------
-    >>> np.nansum(1)
-    1
-    >>> np.nansum([1])
-    1
-    >>> np.nansum([1, np.nan])
-    1.0
-    >>> a = np.array([[1, 1], [1, np.nan]])
-    >>> np.nansum(a)
-    3.0
-    >>> np.nansum(a, axis=0)
-    array([ 2.,  1.])
-
-    When positive infinity and negative infinity are present
-
-    >>> np.nansum([1, np.nan, np.inf])
-    inf
-    >>> np.nansum([1, np.nan, np.NINF])
-    -inf
-    >>> np.nansum([1, np.nan, np.inf, np.NINF])
-    nan
-
-    """
-    return _nanop(np.sum, 0, a, axis)
-
-def nanmin(a, axis=None):
-    """
-    Return the minimum of an array or minimum along an axis, ignoring any NaNs.
-
-    Parameters
-    ----------
-    a : array_like
-        Array containing numbers whose minimum is desired. If `a` is not
-        an array, a conversion is attempted.
-    axis : int, optional
-        Axis along which the minimum is computed. The default is to compute
-        the minimum of the flattened array.
-
-    Returns
-    -------
-    nanmin : ndarray
-        An array with the same shape as `a`, with the specified axis removed.
-        If `a` is a 0-d array, or if axis is None, an ndarray scalar is
-        returned.  The same dtype as `a` is returned.
-
-    See Also
-    --------
-    nanmax :
-        The maximum value of an array along a given axis, ignoring any NaNs.
-    amin :
-        The minimum value of an array along a given axis, propagating any NaNs.
-    fmin :
-        Element-wise minimum of two arrays, ignoring any NaNs.
-    minimum :
-        Element-wise minimum of two arrays, propagating any NaNs.
-    isnan :
-        Shows which elements are Not a Number (NaN).
-    isfinite:
-        Shows which elements are neither NaN nor infinity.
-        
-    amax, fmax, maximum
-
-    Notes
-    -----
-    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
-    (IEEE 754). This means that Not a Number is not equivalent to infinity.
-    Positive infinity is treated as a very large number and negative infinity
-    is treated as a very small (i.e. negative) number.
-
-    If the input has a integer type the function is equivalent to np.min.
-
-    Examples
-    --------
-    >>> a = np.array([[1, 2], [3, np.nan]])
-    >>> np.nanmin(a)
-    1.0
-    >>> np.nanmin(a, axis=0)
-    array([ 1.,  2.])
-    >>> np.nanmin(a, axis=1)
-    array([ 1.,  3.])
-
-    When positive infinity and negative infinity are present:
-
-    >>> np.nanmin([1, 2, np.nan, np.inf])
-    1.0
-    >>> np.nanmin([1, 2, np.nan, np.NINF])
-    -inf
-
-    """
-    a = np.asanyarray(a)
-    if axis is not None:
-        return np.fmin.reduce(a, axis)
-    else:
-        return np.fmin.reduce(a.flat)
-
-def nanargmin(a, axis=None):
-    """
-    Return indices of the minimum values over an axis, ignoring NaNs.
-
-    Parameters
-    ----------
-    a : array_like
-        Input data.
-    axis : int, optional
-        Axis along which to operate.  By default flattened input is used.
-
-    Returns
-    -------
-    index_array : ndarray
-        An array of indices or a single index value.
-
-    See Also
-    --------
-    argmin, nanargmax
-
-    Examples
-    --------
-    >>> a = np.array([[np.nan, 4], [2, 3]])
-    >>> np.argmin(a)
-    0
-    >>> np.nanargmin(a)
-    2
-    >>> np.nanargmin(a, axis=0)
-    array([1, 1])
-    >>> np.nanargmin(a, axis=1)
-    array([1, 0])
-
-    """
-    return _nanop(np.argmin, np.inf, a, axis)
-
-def nanmax(a, axis=None):
-    """
-    Return the maximum of an array or maximum along an axis, ignoring any NaNs.
-
-    Parameters
-    ----------
-    a : array_like
-        Array containing numbers whose maximum is desired. If `a` is not
-        an array, a conversion is attempted.
-    axis : int, optional
-        Axis along which the maximum is computed. The default is to compute
-        the maximum of the flattened array.
-
-    Returns
-    -------
-    nanmax : ndarray
-        An array with the same shape as `a`, with the specified axis removed.
-        If `a` is a 0-d array, or if axis is None, an ndarray scalar is
-        returned.  The same dtype as `a` is returned.
-
-    See Also
-    --------
-    nanmin :
-        The minimum value of an array along a given axis, ignoring any NaNs.
-    amax :
-        The maximum value of an array along a given axis, propagating any NaNs.
-    fmax :
-        Element-wise maximum of two arrays, ignoring any NaNs.
-    maximum :
-        Element-wise maximum of two arrays, propagating any NaNs.
-    isnan :
-        Shows which elements are Not a Number (NaN).
-    isfinite:
-        Shows which elements are neither NaN nor infinity.
-        
-    amin, fmin, minimum
-
-    Notes
-    -----
-    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
-    (IEEE 754). This means that Not a Number is not equivalent to infinity.
-    Positive infinity is treated as a very large number and negative infinity
-    is treated as a very small (i.e. negative) number.
-
-    If the input has a integer type the function is equivalent to np.max.
-
-    Examples
-    --------
-    >>> a = np.array([[1, 2], [3, np.nan]])
-    >>> np.nanmax(a)
-    3.0
-    >>> np.nanmax(a, axis=0)
-    array([ 3.,  2.])
-    >>> np.nanmax(a, axis=1)
-    array([ 2.,  3.])
-
-    When positive infinity and negative infinity are present:
-
-    >>> np.nanmax([1, 2, np.nan, np.NINF])
-    2.0
-    >>> np.nanmax([1, 2, np.nan, np.inf])
-    inf
-
-    """
-    a = np.asanyarray(a)
-    if axis is not None:
-        return np.fmax.reduce(a, axis)
-    else:
-        return np.fmax.reduce(a.flat)
-
-def nanargmax(a, axis=None):
-    """
-    Return indices of the maximum values over an axis, ignoring NaNs.
-
-    Parameters
-    ----------
-    a : array_like
-        Input data.
-    axis : int, optional
-        Axis along which to operate.  By default flattened input is used.
-
-    Returns
-    -------
-    index_array : ndarray
-        An array of indices or a single index value.
-
-    See Also
-    --------
-    argmax, nanargmin
-
-    Examples
-    --------
-    >>> a = np.array([[np.nan, 4], [2, 3]])
-    >>> np.argmax(a)
-    0
-    >>> np.nanargmax(a)
-    1
-    >>> np.nanargmax(a, axis=0)
-    array([1, 0])
-    >>> np.nanargmax(a, axis=1)
-    array([1, 1])
-
-    """
-    return _nanop(np.argmax, -np.inf, a, axis)
-
 def disp(mesg, device=None, linefeed=True):
     """
     Display a message on a device.
diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py
new file mode 100644
index 000000000..81f5aee2e
--- /dev/null
+++ b/numpy/lib/nanfunctions.py
@@ -0,0 +1,812 @@
+"""
+Functions that ignore NaN.
+
+Functions
+---------
+
+- `nanmin` -- minimum non-NaN value
+- `nanmax` -- maximum non-NaN value
+- `nanargmin` -- index of minimum non-NaN value
+- `nanargmax` -- index of maximum non-NaN value
+- `nansum` -- sum of non-NaN values
+- `nanmean` -- mean of non-NaN values
+- `nanvar` -- variance of non-NaN values
+- `nanstd` -- standard deviation of non-NaN values
+
+Classes
+-------
+- `NanWarning` -- Warning raised by nanfunctions
+
+"""
+from __future__ import division, absolute_import, print_function
+
+import warnings
+import numpy as np
+
+__all__ = [
+    'nansum', 'nanmax', 'nanmin', 'nanargmax', 'nanargmin', 'nanmean',
+    'nanvar', 'nanstd', 'NanWarning'
+    ]
+
+class NanWarning(RuntimeWarning): pass
+
+
+def _replace_nan(a, val):
+    """
+    If `a` is of inexact type, make a copy of `a`, replace NaNs with
+    the `val` value, and return the copy together with a boolean mask
+    marking the locations where NaNs were present. If `a` is not of
+    inexact type, do nothing and return `a` together with a mask of None.
+
+    Parameters
+    ----------
+    a : array-like
+        Input array.
+    val : float
+        NaN values are set to val before doing the operation.
+
+    Returns
+    -------
+    y : ndarray
+        If `a` is of inexact type, return a copy of `a` with the NaNs
+        replaced by the fill value, otherwise return `a`.
+    mask: {bool, None}
+        If `a` is of inexact type, return a boolean mask marking locations of
+        NaNs, otherwise return None.
+
+    """
+    is_new = not isinstance(a, np.ndarray)
+    if is_new:
+        a = np.array(a)
+    if not issubclass(a.dtype.type, np.inexact):
+        return a, None
+    if not is_new:
+        # need copy
+        a = np.array(a, subok=True)
+
+    mask = np.isnan(a)
+    np.copyto(a, val, where=mask)
+    return a, mask
+
+
+def _copyto(a, val, mask):
+    """
+    Replace values in `a` with NaN where `mask` is True.  This differs from
+    copyto in that it will deal with the case where `a` is a numpy scalar.
+
+    Parameters
+    ----------
+    a : ndarray or numpy scalar
+        Array or numpy scalar some of whose values are to be replaced
+        by val.
+    val : numpy scalar
+        Value used a replacement.
+    mask : ndarray, scalar
+        Boolean array. Where True the corresponding element of `a` is
+        replaced by `val`. Broadcasts.
+
+    Returns
+    -------
+    res : ndarray, scalar
+        Array with elements replaced or scalar `val`.
+
+    """
+    if isinstance(a, np.ndarray):
+        np.copyto(a, val, where=mask, casting='unsafe')
+    else:
+        a = a.dtype.type(val)
+    return a
+
+
+def _divide_by_count(a, b, out=None):
+    """
+    Compute a/b ignoring invalid results. If `a` is an array the division
+    is done in place. If `a` is a scalar, then its type is preserved in the
+    output. If out is None, then then a is used instead so that the
+    division is in place.
+
+    Parameters
+    ----------
+    a : {ndarray, numpy scalar}
+        Numerator. Expected to be of inexact type but not checked.
+    b : {ndarray, numpy scalar}
+        Denominator.
+    out : ndarray, optional
+        Alternate output array in which to place the result.  The default
+        is ``None``; if provided, it must have the same shape as the
+        expected output, but the type will be cast if necessary.
+
+    Returns
+    -------
+    ret : {ndarray, numpy scalar}
+        The return value is a/b. If `a` was an ndarray the division is done
+        in place. If `a` is a numpy scalar, the division preserves its type.
+
+    """
+    with np.errstate(invalid='ignore'):
+        if isinstance(a, np.ndarray):
+            if out is None:
+                return np.divide(a, b, out=a, casting='unsafe')
+            else:
+                return np.divide(a, b, out=out, casting='unsafe')
+        else:
+            if out is None:
+                return a.dtype.type(a / b)
+            else:
+                # This is questionable, but currently a numpy scalar can
+                # be output to a zero dimensional array.
+                return np.divide(a, b, out=out, casting='unsafe')
+
+
+def nanmin(a, axis=None, out=None, keepdims=False):
+    """
+    Return the minimum of an array or minimum along an axis, ignoring any
+    NaNs.
+
+    Parameters
+    ----------
+    a : array_like
+        Array containing numbers whose minimum is desired. If `a` is not
+        an array, a conversion is attempted.
+    axis : int, optional
+        Axis along which the minimum is computed. The default is to compute
+        the minimum of the flattened array.
+    out : ndarray, optional
+        Alternate output array in which to place the result.  The default
+        is ``None``; if provided, it must have the same shape as the
+        expected output, but the type will be cast if necessary.
+        See `doc.ufuncs` for details.
+
+        .. versionadded:: 1.8.0
+    keepdims : bool, optional
+        If this is set to True, the axes which are reduced are left
+        in the result as dimensions with size one. With this option,
+        the result will broadcast correctly against the original `a`.
+
+        .. versionadded:: 1.8.0
+
+    Returns
+    -------
+    nanmin : ndarray
+        An array with the same shape as `a`, with the specified axis removed.
+        If `a` is a 0-d array, or if axis is None, an ndarray scalar is
+        returned.  The same dtype as `a` is returned.
+
+    See Also
+    --------
+    nanmax :
+        The maximum value of an array along a given axis, ignoring any NaNs.
+    amin :
+        The minimum value of an array along a given axis, propagating any NaNs.
+    fmin :
+        Element-wise minimum of two arrays, ignoring any NaNs.
+    minimum :
+        Element-wise minimum of two arrays, propagating any NaNs.
+    isnan :
+        Shows which elements are Not a Number (NaN).
+    isfinite:
+        Shows which elements are neither NaN nor infinity.
+
+    amax, fmax, maximum
+
+    Notes
+    -----
+    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
+    (IEEE 754). This means that Not a Number is not equivalent to infinity.
+    Positive infinity is treated as a very large number and negative infinity
+    is treated as a very small (i.e. negative) number.
+
+    If the input has a integer type the function is equivalent to np.min.
+
+    Examples
+    --------
+    >>> a = np.array([[1, 2], [3, np.nan]])
+    >>> np.nanmin(a)
+    1.0
+    >>> np.nanmin(a, axis=0)
+    array([ 1.,  2.])
+    >>> np.nanmin(a, axis=1)
+    array([ 1.,  3.])
+
+    When positive infinity and negative infinity are present:
+
+    >>> np.nanmin([1, 2, np.nan, np.inf])
+    1.0
+    >>> np.nanmin([1, 2, np.nan, np.NINF])
+    -inf
+
+    """
+    return np.fmin.reduce(a, axis=axis, out=out, keepdims=keepdims)
+
+
+def nanmax(a, axis=None, out=None, keepdims=False):
+    """
+    Return the maximum of an array or maximum along an axis, ignoring any NaNs.
+
+    Parameters
+    ----------
+    a : array_like
+        Array containing numbers whose maximum is desired. If `a` is not
+        an array, a conversion is attempted.
+    axis : int, optional
+        Axis along which the maximum is computed. The default is to compute
+        the maximum of the flattened array.
+    out : ndarray, optional
+        Alternate output array in which to place the result.  The default
+        is ``None``; if provided, it must have the same shape as the
+        expected output, but the type will be cast if necessary.
+        See `doc.ufuncs` for details.
+
+        .. versionadded:: 1.8.0
+    keepdims : bool, optional
+        If this is set to True, the axes which are reduced are left
+        in the result as dimensions with size one. With this option,
+        the result will broadcast correctly against the original `a`.
+
+        .. versionadded:: 1.8.0
+
+    Returns
+    -------
+    nanmax : ndarray
+        An array with the same shape as `a`, with the specified axis removed.
+        If `a` is a 0-d array, or if axis is None, an ndarray scalar is
+        returned.  The same dtype as `a` is returned.
+
+    See Also
+    --------
+    nanmin :
+        The minimum value of an array along a given axis, ignoring any NaNs.
+    amax :
+        The maximum value of an array along a given axis, propagating any NaNs.
+    fmax :
+        Element-wise maximum of two arrays, ignoring any NaNs.
+    maximum :
+        Element-wise maximum of two arrays, propagating any NaNs.
+    isnan :
+        Shows which elements are Not a Number (NaN).
+    isfinite:
+        Shows which elements are neither NaN nor infinity.
+
+    amin, fmin, minimum
+
+    Notes
+    -----
+    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
+    (IEEE 754). This means that Not a Number is not equivalent to infinity.
+    Positive infinity is treated as a very large number and negative infinity
+    is treated as a very small (i.e. negative) number.
+
+    If the input has a integer type the function is equivalent to np.max.
+
+    Examples
+    --------
+    >>> a = np.array([[1, 2], [3, np.nan]])
+    >>> np.nanmax(a)
+    3.0
+    >>> np.nanmax(a, axis=0)
+    array([ 3.,  2.])
+    >>> np.nanmax(a, axis=1)
+    array([ 2.,  3.])
+
+    When positive infinity and negative infinity are present:
+
+    >>> np.nanmax([1, 2, np.nan, np.NINF])
+    2.0
+    >>> np.nanmax([1, 2, np.nan, np.inf])
+    inf
+
+    """
+    return np.fmax.reduce(a, axis=axis, out=out, keepdims=keepdims)
+
+
+def nanargmin(a, axis=None):
+    """
+    Return the indices of the minimum values in the specified axis ignoring
+    NaNs. For all-NaN slices, the negative number ``np.iinfo('intp').min``
+    is returned. It is platform dependent. Warning: the results cannot be
+    trusted if a slice contains only NaNs and Infs.
+
+    Parameters
+    ----------
+    a : array_like
+        Input data.
+    axis : int, optional
+        Axis along which to operate.  By default flattened input is used.
+
+    Returns
+    -------
+    index_array : ndarray
+        An array of indices or a single index value.
+
+    See Also
+    --------
+    argmin, nanargmax
+
+    Examples
+    --------
+    >>> a = np.array([[np.nan, 4], [2, 3]])
+    >>> np.argmin(a)
+    0
+    >>> np.nanargmin(a)
+    2
+    >>> np.nanargmin(a, axis=0)
+    array([1, 1])
+    >>> np.nanargmin(a, axis=1)
+    array([1, 0])
+
+    """
+    a, mask = _replace_nan(a, np.inf)
+    if mask is None:
+        return np.argmin(a, axis)
+    # May later want to do something special for all nan slices.
+    mask = mask.all(axis=axis)
+    ind = np.argmin(a, axis)
+    if mask.any():
+        warnings.warn("All NaN axis detected.", NanWarning)
+        ind =_copyto(ind, np.iinfo(np.intp).min, mask)
+    return ind
+
+
+def nanargmax(a, axis=None):
+    """
+    Return the indices of the maximum values in the specified axis ignoring
+    NaNs. For  all-NaN slices, the negative number ``np.iinfo('intp').min``
+    is returned. It is platform dependent. Warning: the results cannot be
+    trusted if a slice contains only NaNs and -Infs.
+
+
+    Parameters
+    ----------
+    a : array_like
+        Input data.
+    axis : int, optional
+        Axis along which to operate.  By default flattened input is used.
+
+    Returns
+    -------
+    index_array : ndarray
+        An array of indices or a single index value.
+
+    See Also
+    --------
+    argmax, nanargmin
+
+    Examples
+    --------
+    >>> a = np.array([[np.nan, 4], [2, 3]])
+    >>> np.argmax(a)
+    0
+    >>> np.nanargmax(a)
+    1
+    >>> np.nanargmax(a, axis=0)
+    array([1, 0])
+    >>> np.nanargmax(a, axis=1)
+    array([1, 1])
+
+    """
+    a, mask = _replace_nan(a, -np.inf)
+    if mask is None:
+        return np.argmax(a, axis)
+    # May later want to do something special for all nan slices.
+    mask = mask.all(axis=axis)
+    ind = np.argmax(a, axis)
+    if mask.any():
+        warnings.warn("All NaN axis detected.", NanWarning)
+        ind = _copyto(ind, np.iinfo(np.intp).min, mask)
+    return ind
+
+
+def nansum(a, axis=None, dtype=None, out=None, keepdims=0):
+    """
+    Return the sum of array elements over a given axis treating
+    Not a Numbers (NaNs) as zero.
+
+    FutureWarning: In Numpy versions <= 1.8 Nan is returned for slices that
+    are all-NaN or empty. In later versions zero will be returned.
+
+    Parameters
+    ----------
+    a : array_like
+        Array containing numbers whose sum is desired. If `a` is not an
+        array, a conversion is attempted.
+    axis : int, optional
+        Axis along which the sum is computed. The default is to compute
+        the sum of the flattened array.
+    dtype : data-type, optional
+        Type to use in computing the sum.  For integer inputs, the default
+        is the same as `int64`. For inexact inputs, it must be inexact.
+
+        .. versionadded:: 1.8.0
+    out : ndarray, optional
+        Alternate output array in which to place the result.  The default
+        is ``None``. If provided, it must have the same shape as the
+        expected output, but the type will be cast if necessary.
+        See `doc.ufuncs` for details. The casting of NaN to integer can
+        yield unexpected results.
+
+        .. versionadded:: 1.8.0
+    keepdims : bool, optional
+        If True, the axes which are reduced are left in the result as
+        dimensions with size one. With this option, the result will
+        broadcast correctly against the original `arr`.
+
+        .. versionadded:: 1.8.0
+
+    Returns
+    -------
+    y : ndarray or numpy scalar
+
+    See Also
+    --------
+    numpy.sum : Sum across array propagating NaNs.
+    isnan : Show which elements are NaN.
+    isfinite: Show which elements are not NaN or +/-inf.
+
+    Notes
+    -----
+    Numpy uses the IEEE Standard for Binary Floating-Point for Arithmetic
+    (IEEE 754). This means that Not a Number is not equivalent to infinity.
+    If positive or negative infinity are present the result is positive or
+    negative infinity. But if both positive and negative infinity are present,
+    the result is Not A Number (NaN).
+
+    Arithmetic is modular when using integer types (all elements of `a` must
+    be finite i.e. no elements that are NaNs, positive infinity and negative
+    infinity because NaNs are floating point types), and no error is raised
+    on overflow.
+
+
+    Examples
+    --------
+    >>> np.nansum(1)
+    1
+    >>> np.nansum([1])
+    1
+    >>> np.nansum([1, np.nan])
+    1.0
+    >>> a = np.array([[1, 1], [1, np.nan]])
+    >>> np.nansum(a)
+    3.0
+    >>> np.nansum(a, axis=0)
+    array([ 2.,  1.])
+    >>> np.nansum([1, np.nan, np.inf])
+    inf
+    >>> np.nansum([1, np.nan, np.NINF])
+    -inf
+    >>> np.nansum([1, np.nan, np.inf, -np.inf]) # both +/- infinity present
+    nan
+
+    """
+    a, mask = _replace_nan(a, 0)
+    # In version 1.9 uncomment the following line and delete the rest.
+    #return a.sum(axis, dtype, out, keepdims)
+    warnings.warn("In Numpy 1.9 the sum along empty slices will be zero.",
+            FutureWarning)
+
+    if mask is None:
+        return a.sum(axis, dtype, out, keepdims)
+    mask = mask.all(axis, keepdims=keepdims)
+    tot = np.add.reduce(a, axis, dtype, out, keepdims)
+    if mask.any():
+        tot = _copyto(tot, np.nan, mask)
+    return tot
+
+
+def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
+    """
+    Compute the arithmetic mean along the specified axis, ignoring NaNs.
+
+    Returns the average of the array elements.  The average is taken over
+    the flattened array by default, otherwise over the specified axis.
+    `float64` intermediate and return values are used for integer inputs.
+
+    For all-NaN slices, NaN is returned and a `NanWarning` is raised.
+
+    .. versionadded:: 1.8.0
+
+    Parameters
+    ----------
+    a : array_like
+        Array containing numbers whose mean is desired. If `a` is not an
+        array, a conversion is attempted.
+    axis : int, optional
+        Axis along which the means are computed. The default is to compute
+        the mean of the flattened array.
+    dtype : data-type, optional
+        Type to use in computing the mean.  For integer inputs, the default
+        is `float64`; for inexact inputs, it is the same as the
+        input dtype.
+    out : ndarray, optional
+        Alternate output array in which to place the result.  The default
+        is ``None``; if provided, it must have the same shape as the
+        expected output, but the type will be cast if necessary.
+        See `doc.ufuncs` for details.
+    keepdims : bool, optional
+        If this is set to True, the axes which are reduced are left
+        in the result as dimensions with size one. With this option,
+        the result will broadcast correctly against the original `arr`.
+
+    Returns
+    -------
+    m : ndarray, see dtype parameter above
+        If `out=None`, returns a new array containing the mean values,
+        otherwise a reference to the output array is returned. Nan is
+        returned for slices that contain only NaNs.
+
+    See Also
+    --------
+    average : Weighted average
+    mean : Arithmetic mean taken while not ignoring NaNs
+    var, nanvar
+
+    Notes
+    -----
+    The arithmetic mean is the sum of the non-NaN elements along the axis
+    divided by the number of non-NaN elements.
+
+    Note that for floating-point input, the mean is computed using the
+    same precision the input has.  Depending on the input data, this can
+    cause the results to be inaccurate, especially for `float32`.
+    Specifying a higher-precision accumulator using the `dtype` keyword
+    can alleviate this issue.
+
+    Examples
+    --------
+    >>> a = np.array([[1, np.nan], [3, 4]])
+    >>> np.nanmean(a)
+    2.6666666666666665
+    >>> np.nanmean(a, axis=0)
+    array([ 2.,  4.])
+    >>> np.nanmean(a, axis=1)
+    array([ 1.,  3.5])
+
+    """
+    arr, mask = _replace_nan(a, 0)
+    if mask is None:
+        return np.mean(arr, axis, dtype=dtype, out=out, keepdims=keepdims)
+
+    if dtype is not None:
+        dtype = np.dtype(dtype)
+    if dtype is not None and not issubclass(dtype.type, np.inexact):
+        raise TypeError("If a is inexact, then dtype must be inexact")
+    if out is not None and not issubclass(out.dtype.type, np.inexact):
+        raise TypeError("If a is inexact, then out must be inexact")
+
+    # The warning context speeds things up.
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        cnt = np.add.reduce(~mask, axis, dtype=np.intp, keepdims=keepdims)
+        tot = np.add.reduce(arr, axis, dtype=dtype, out=out, keepdims=keepdims)
+        avg = _divide_by_count(tot, cnt, out=out)
+
+    isbad = (cnt == 0)
+    if isbad.any():
+        warnings.warn("Mean of empty slice", NanWarning)
+        # NaN is the only possible bad value, so no further
+        # action is needed to handle bad results.
+    return avg
+
+
+def nanvar(a, axis=None, dtype=None, out=None, ddof=0,
+           keepdims=False):
+    """
+    Compute the variance along the specified axis, while ignoring NaNs.
+
+    Returns the variance of the array elements, a measure of the spread of a
+    distribution.  The variance is computed for the flattened array by
+    default, otherwise over the specified axis.
+
+    For all-NaN slices, NaN is returned and a `NanWarning` is raised.
+
+    .. versionadded:: 1.8.0
+
+    Parameters
+    ----------
+    a : array_like
+        Array containing numbers whose variance is desired.  If `a` is not an
+        array, a conversion is attempted.
+    axis : int, optional
+        Axis along which the variance is computed.  The default is to compute
+        the variance of the flattened array.
+    dtype : data-type, optional
+        Type to use in computing the variance.  For arrays of integer type
+        the default is `float32`; for arrays of float types it is the same as
+        the array type.
+    out : ndarray, optional
+        Alternate output array in which to place the result.  It must have
+        the same shape as the expected output, but the type is cast if
+        necessary.
+    ddof : int, optional
+        "Delta Degrees of Freedom": the divisor used in the calculation is
+        ``N - ddof``, where ``N`` represents the number of non-NaN
+        elements. By default `ddof` is zero.
+    keepdims : bool, optional
+        If this is set to True, the axes which are reduced are left
+        in the result as dimensions with size one. With this option,
+        the result will broadcast correctly against the original `arr`.
+
+    Returns
+    -------
+    variance : ndarray, see dtype parameter above
+        If `out` is None, return a new array containing the variance,
+        otherwise return a reference to the output array. If ddof is >= the
+        number of non-NaN elements in a slice or the slice contains only
+        NaNs, then the result for that slice is NaN.
+
+    See Also
+    --------
+    std : Standard deviation
+    mean : Average
+    var : Variance while not ignoring NaNs
+    nanstd, nanmean
+    numpy.doc.ufuncs : Section "Output arguments"
+
+    Notes
+    -----
+    The variance is the average of the squared deviations from the mean,
+    i.e.,  ``var = mean(abs(x - x.mean())**2)``.
+
+    The mean is normally calculated as ``x.sum() / N``, where ``N = len(x)``.
+    If, however, `ddof` is specified, the divisor ``N - ddof`` is used
+    instead.  In standard statistical practice, ``ddof=1`` provides an
+    unbiased estimator of the variance of a hypothetical infinite population.
+    ``ddof=0`` provides a maximum likelihood estimate of the variance for
+    normally distributed variables.
+
+    Note that for complex numbers, the absolute value is taken before
+    squaring, so that the result is always real and nonnegative.
+
+    For floating-point input, the variance is computed using the same
+    precision the input has.  Depending on the input data, this can cause
+    the results to be inaccurate, especially for `float32` (see example
+    below).  Specifying a higher-accuracy accumulator using the ``dtype``
+    keyword can alleviate this issue.
+
+    Examples
+    --------
+    >>> a = np.array([[1, np.nan], [3, 4]])
+    >>> np.var(a)
+    1.5555555555555554
+    >>> np.nanvar(a, axis=0)
+    array([ 1.,  0.])
+    >>> np.nanvar(a, axis=1)
+    array([ 0.,  0.25])
+
+    """
+    arr, mask = _replace_nan(a, 0)
+    if mask is None:
+        return np.var(arr, axis, dtype=dtype, out=out, keepdims=keepdims)
+
+    if dtype is not None:
+        dtype = np.dtype(dtype)
+    if dtype is not None and not issubclass(dtype.type, np.inexact):
+        raise TypeError("If a is inexact, then dtype must be inexact")
+    if out is not None and not issubclass(out.dtype.type, np.inexact):
+        raise TypeError("If a is inexact, then out must be inexact")
+
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+
+        # Compute mean
+        cnt = np.add.reduce(~mask, axis, dtype=np.intp, keepdims=True)
+        tot = np.add.reduce(arr, axis, dtype=dtype, keepdims=True)
+        avg = np.divide(tot, cnt, out=tot)
+
+        # Compute squared deviation from mean.
+        x = arr - avg
+        np.copyto(x, 0, where=mask)
+        if issubclass(arr.dtype.type, np.complexfloating):
+            sqr = np.multiply(x, x.conj(), out=x).real
+        else:
+            sqr = np.multiply(x, x, out=x)
+
+        # adjust cnt.
+        if not keepdims:
+            cnt = cnt.squeeze(axis)
+        cnt -= ddof
+
+        # Compute variance.
+        var = np.add.reduce(sqr, axis, dtype=dtype, out=out, keepdims=keepdims)
+        var = _divide_by_count(var, cnt)
+
+    isbad = (cnt <= 0)
+    if isbad.any():
+        warnings.warn("Degrees of freedom <= 0 for slice.", NanWarning)
+        # NaN, inf, or negative numbers are all possible bad
+        # values, so explicitly replace them with NaN.
+        var = _copyto(var, np.nan, isbad)
+    return var
+
+
+def nanstd(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
+    """
+    Compute the standard deviation along the specified axis, while
+    ignoring NaNs.
+
+    Returns the standard deviation, a measure of the spread of a distribution,
+    of the non-NaN array elements. The standard deviation is computed for the
+    flattened array by default, otherwise over the specified axis.
+
+    For all-NaN slices, NaN is returned and a `NanWarning` is raised.
+
+    .. versionadded:: 1.8.0
+
+    Parameters
+    ----------
+    a : array_like
+        Calculate the standard deviation of the non-NaN values.
+    axis : int, optional
+        Axis along which the standard deviation is computed. The default is
+        to compute the standard deviation of the flattened array.
+    dtype : dtype, optional
+        Type to use in computing the standard deviation. For arrays of
+        integer type the default is float64, for arrays of float types it is
+        the same as the array type.
+    out : ndarray, optional
+        Alternative output array in which to place the result. It must have
+        the same shape as the expected output but the type (of the calculated
+        values) will be cast if necessary.
+    ddof : int, optional
+        Means Delta Degrees of Freedom.  The divisor used in calculations
+        is ``N - ddof``, where ``N`` represents the number of non-NaN
+        elements.  By default `ddof` is zero.
+    keepdims : bool, optional
+        If this is set to True, the axes which are reduced are left
+        in the result as dimensions with size one. With this option,
+        the result will broadcast correctly against the original `arr`.
+
+    Returns
+    -------
+    standard_deviation : ndarray, see dtype parameter above.
+        If `out` is None, return a new array containing the standard
+        deviation, otherwise return a reference to the output array. If
+        ddof is >= the number of non-NaN elements in a slice or the slice
+        contains only NaNs, then the result for that slice is NaN.
+
+    See Also
+    --------
+    var, mean, std
+    nanvar, nanmean
+    numpy.doc.ufuncs : Section "Output arguments"
+
+    Notes
+    -----
+    The standard deviation is the square root of the average of the squared
+    deviations from the mean, i.e., ``std = sqrt(mean(abs(x - x.mean())**2))``.
+
+    The average squared deviation is normally calculated as
+    ``x.sum() / N``, where ``N = len(x)``.  If, however, `ddof` is specified,
+    the divisor ``N - ddof`` is used instead. In standard statistical
+    practice, ``ddof=1`` provides an unbiased estimator of the variance
+    of the infinite population. ``ddof=0`` provides a maximum likelihood
+    estimate of the variance for normally distributed variables. The
+    standard deviation computed in this function is the square root of
+    the estimated variance, so even with ``ddof=1``, it will not be an
+    unbiased estimate of the standard deviation per se.
+
+    Note that, for complex numbers, `std` takes the absolute
+    value before squaring, so that the result is always real and nonnegative.
+
+    For floating-point input, the *std* is computed using the same
+    precision the input has. Depending on the input data, this can cause
+    the results to be inaccurate, especially for float32 (see example below).
+    Specifying a higher-accuracy accumulator using the `dtype` keyword can
+    alleviate this issue.
+
+    Examples
+    --------
+    >>> a = np.array([[1, np.nan], [3, 4]])
+    >>> np.nanstd(a)
+    1.247219128924647
+    >>> np.nanstd(a, axis=0)
+    array([ 1.,  0.])
+    >>> np.nanstd(a, axis=1)
+    array([ 0.,  0.5])
+
+    """
+    var = nanvar(a, axis, dtype, out, ddof, keepdims)
+    if isinstance(var, np.ndarray):
+        std = np.sqrt(var, out=var)
+    else:
+        std = var.dtype.type(np.sqrt(var))
+    return std
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 814743442..cf303993b 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -3,10 +3,10 @@ from __future__ import division, absolute_import, print_function
 import warnings
 import numpy as np
 from numpy.testing import (
-        run_module_suite, TestCase, assert_, assert_equal,
-        assert_array_equal, assert_almost_equal, assert_array_almost_equal,
-        assert_raises, assert_allclose, assert_array_max_ulp, assert_warns
-        )
+    run_module_suite, TestCase, assert_, assert_equal, assert_array_equal,
+    assert_almost_equal, assert_array_almost_equal, assert_raises,
+    assert_allclose, assert_array_max_ulp, assert_warns
+    )
 from numpy.random import rand
 from numpy.lib import *
 from numpy.compat import long
@@ -1111,127 +1111,6 @@ class TestCheckFinite(TestCase):
         assert_(a.dtype == np.float64)
 
 
-class TestNaNFuncts(TestCase):
-    def setUp(self):
-        self.A = np.array([[[ np.nan, 0.01319214, 0.01620964],
-                         [ 0.11704017, np.nan, 0.75157887],
-                         [ 0.28333658, 0.1630199 , np.nan       ]],
-                        [[ 0.59541557, np.nan, 0.37910852],
-                         [ np.nan, 0.87964135, np.nan       ],
-                         [ 0.70543747, np.nan, 0.34306596]],
-                        [[ 0.72687499, 0.91084584, np.nan       ],
-                         [ 0.84386844, 0.38944762, 0.23913896],
-                         [ np.nan, 0.37068164, 0.33850425]]])
-
-    def test_nansum(self):
-        assert_almost_equal(nansum(self.A), 8.0664079100000006)
-        assert_almost_equal(nansum(self.A, 0),
-                            np.array([[ 1.32229056, 0.92403798, 0.39531816],
-                                   [ 0.96090861, 1.26908897, 0.99071783],
-                                   [ 0.98877405, 0.53370154, 0.68157021]]))
-        assert_almost_equal(nansum(self.A, 1),
-                            np.array([[ 0.40037675, 0.17621204, 0.76778851],
-                                   [ 1.30085304, 0.87964135, 0.72217448],
-                                   [ 1.57074343, 1.6709751 , 0.57764321]]))
-        assert_almost_equal(nansum(self.A, 2),
-                            np.array([[ 0.02940178, 0.86861904, 0.44635648],
-                                   [ 0.97452409, 0.87964135, 1.04850343],
-                                   [ 1.63772083, 1.47245502, 0.70918589]]))
-
-    def test_nanmin(self):
-        assert_almost_equal(nanmin(self.A), 0.01319214)
-        assert_almost_equal(nanmin(self.A, 0),
-                            np.array([[ 0.59541557, 0.01319214, 0.01620964],
-                                   [ 0.11704017, 0.38944762, 0.23913896],
-                                   [ 0.28333658, 0.1630199 , 0.33850425]]))
-        assert_almost_equal(nanmin(self.A, 1),
-                            np.array([[ 0.11704017, 0.01319214, 0.01620964],
-                                   [ 0.59541557, 0.87964135, 0.34306596],
-                                   [ 0.72687499, 0.37068164, 0.23913896]]))
-        assert_almost_equal(nanmin(self.A, 2),
-                            np.array([[ 0.01319214, 0.11704017, 0.1630199 ],
-                                   [ 0.37910852, 0.87964135, 0.34306596],
-                                   [ 0.72687499, 0.23913896, 0.33850425]]))
-        assert_(np.isnan(nanmin([np.nan, np.nan])))
-
-    def test_nanargmin(self):
-        assert_almost_equal(nanargmin(self.A), 1)
-        assert_almost_equal(nanargmin(self.A, 0),
-                            np.array([[1, 0, 0],
-                                   [0, 2, 2],
-                                   [0, 0, 2]]))
-        assert_almost_equal(nanargmin(self.A, 1),
-                            np.array([[1, 0, 0],
-                                   [0, 1, 2],
-                                   [0, 2, 1]]))
-        assert_almost_equal(nanargmin(self.A, 2),
-                            np.array([[1, 0, 1],
-                                   [2, 1, 2],
-                                   [0, 2, 2]]))
-
-    def test_nanmax(self):
-        assert_almost_equal(nanmax(self.A), 0.91084584000000002)
-        assert_almost_equal(nanmax(self.A, 0),
-                            np.array([[ 0.72687499, 0.91084584, 0.37910852],
-                                   [ 0.84386844, 0.87964135, 0.75157887],
-                                   [ 0.70543747, 0.37068164, 0.34306596]]))
-        assert_almost_equal(nanmax(self.A, 1),
-                            np.array([[ 0.28333658, 0.1630199 , 0.75157887],
-                                   [ 0.70543747, 0.87964135, 0.37910852],
-                                   [ 0.84386844, 0.91084584, 0.33850425]]))
-        assert_almost_equal(nanmax(self.A, 2),
-                            np.array([[ 0.01620964, 0.75157887, 0.28333658],
-                                   [ 0.59541557, 0.87964135, 0.70543747],
-                                   [ 0.91084584, 0.84386844, 0.37068164]]))
-        assert_(np.isnan(nanmax([np.nan, np.nan])))
-
-    def test_nanmin_allnan_on_axis(self):
-        assert_array_equal(np.isnan(nanmin([[np.nan] * 2] * 3, axis=1)),
-                     [True, True, True])
-
-    def test_nanmin_masked(self):
-        a = np.ma.fix_invalid([[2, 1, 3, np.nan], [5, 2, 3, np.nan]])
-        ctrl_mask = a._mask.copy()
-        test = np.nanmin(a, axis=1)
-        assert_equal(test, [1, 2])
-        assert_equal(a._mask, ctrl_mask)
-        assert_equal(np.isinf(a), np.zeros((2, 4), dtype=bool))
-
-
-class TestNanFunctsIntTypes(TestCase):
-
-    int_types = (
-            np.int8, np.int16, np.int32, np.int64, np.uint8,
-            np.uint16, np.uint32, np.uint64)
-
-    def setUp(self, *args, **kwargs):
-        self.A = np.array([127, 39,  93,  87, 46])
-
-    def integer_arrays(self):
-        for dtype in self.int_types:
-            yield self.A.astype(dtype)
-
-    def test_nanmin(self):
-        min_value = min(self.A)
-        for A in self.integer_arrays():
-            assert_equal(nanmin(A), min_value)
-
-    def test_nanmax(self):
-        max_value = max(self.A)
-        for A in self.integer_arrays():
-            assert_equal(nanmax(A), max_value)
-
-    def test_nanargmin(self):
-        min_arg = np.argmin(self.A)
-        for A in self.integer_arrays():
-            assert_equal(nanargmin(A), min_arg)
-
-    def test_nanargmax(self):
-        max_arg = np.argmax(self.A)
-        for A in self.integer_arrays():
-            assert_equal(nanargmax(A), max_arg)
-
-
 class TestCorrCoef(TestCase):
     A = np.array([[ 0.15391142, 0.18045767, 0.14197213],
                [ 0.70461506, 0.96474128, 0.27906989],
@@ -1278,7 +1157,7 @@ class TestCov(TestCase):
         assert_equal(cov(np.array([]).reshape(0, 2)).shape, (0, 2))
 
 
-class Test_i0(TestCase):
+class Test_I0(TestCase):
     def test_simple(self):
         assert_almost_equal(i0(0.5), np.array(1.0634833707413234))
         A = np.array([ 0.49842636, 0.6969809 , 0.22011976, 0.0155549])
@@ -1596,7 +1475,5 @@ class TestAdd_newdoc_ufunc(TestCase):
         assert_raises(TypeError, add_newdoc_ufunc, np.add, 3)
 
 
-
-
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/tests/test_nanfunctions.py b/numpy/lib/tests/test_nanfunctions.py
new file mode 100644
index 000000000..93a5ef855
--- /dev/null
+++ b/numpy/lib/tests/test_nanfunctions.py
@@ -0,0 +1,417 @@
+from __future__ import division, absolute_import, print_function
+
+import warnings
+
+import numpy as np
+from numpy.testing import (
+    run_module_suite, TestCase, assert_, assert_equal, assert_almost_equal,
+    assert_raises
+    )
+from numpy.lib import (
+    nansum, nanmax, nanargmax, nanargmin, nanmin, nanmean, nanvar, nanstd,
+    NanWarning
+    )
+
+
+_ndat = np.array(
+        [[ 0.6244, np.nan, 0.2692,  0.0116, np.nan, 0.1170],
+         [ 0.5351, 0.9403, np.nan,  0.2100, 0.4759, 0.2833],
+         [ np.nan, np.nan, np.nan,  0.1042, np.nan, 0.5954],
+         [ 0.161 , np.nan, np.nan,  0.1859, 0.3146, np.nan]]
+        )
+
+# rows of _ndat with nans removed
+_rdat = [
+        np.array([ 0.6244, 0.2692, 0.0116, 0.1170]),
+        np.array([ 0.5351, 0.9403, 0.2100, 0.4759, 0.2833]),
+        np.array([ 0.1042, 0.5954]),
+        np.array([ 0.1610, 0.1859, 0.3146])
+       ]
+
+
+class TestNanFunctions_MinMax(TestCase):
+
+    nanfuncs = [nanmin, nanmax]
+    stdfuncs = [np.min, np.max]
+
+    def test_mutation(self):
+        # Check that passed array is not modified.
+        ndat = _ndat.copy()
+        for f in self.nanfuncs:
+            f(ndat)
+            assert_equal(ndat, _ndat)
+
+    def test_keepdims(self):
+        mat = np.eye(3)
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for axis in [None, 0, 1]:
+                tgt = rf(mat, axis=axis, keepdims=True)
+                res = nf(mat, axis=axis, keepdims=True)
+                assert_(res.ndim == tgt.ndim)
+
+    def test_out(self):
+        mat = np.eye(3)
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            resout = np.zeros(3)
+            tgt = rf(mat, axis=1)
+            res = nf(mat, axis=1, out=resout)
+            assert_almost_equal(res, resout)
+            assert_almost_equal(res, tgt)
+
+    def test_dtype_from_input(self):
+        codes = 'efdgFDG'
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for c in codes:
+                mat = np.eye(3, dtype=c)
+                tgt = rf(mat, axis=1).dtype.type
+                res = nf(mat, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                tgt = rf(mat, axis=None).dtype.type
+                res = nf(mat, axis=None).dtype.type
+                assert_(res is tgt)
+
+    def test_result_values(self):
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            tgt = [rf(d) for d in _rdat]
+            res = nf(_ndat, axis=1)
+            assert_almost_equal(res, tgt)
+
+    def test_allnans(self):
+        mat = np.array([np.nan]*9).reshape(3, 3)
+        for f in self.nanfuncs:
+            for axis in [None, 0, 1]:
+                assert_(np.isnan(f(mat, axis=axis)).all())
+
+    def test_masked(self):
+        mat = np.ma.fix_invalid(_ndat)
+        msk = mat._mask.copy()
+        for f in [nanmin]:
+            res = f(mat, axis=1)
+            tgt = f(_ndat, axis=1)
+            assert_equal(res, tgt)
+            assert_equal(mat._mask, msk)
+            assert_(not np.isinf(mat).any())
+
+
+class TestNanFunctions_ArgminArgmax(TestCase):
+
+    nanfuncs = [nanargmin, nanargmax]
+
+    def test_mutation(self):
+        # Check that passed array is not modified.
+        ndat = _ndat.copy()
+        for f in self.nanfuncs:
+            f(ndat)
+            assert_equal(ndat, _ndat)
+
+    def test_result_values(self):
+        for f, fcmp in zip(self.nanfuncs, [np.greater, np.less]):
+            for row in _ndat:
+                with warnings.catch_warnings():
+                    warnings.simplefilter('ignore')
+                    ind = f(row)
+                    val = row[ind]
+                    # comparing with NaN is tricky as the result
+                    # is always false except for NaN != NaN
+                    assert_(not np.isnan(val))
+                    assert_(not fcmp(val, row).any())
+                    assert_(not np.equal(val, row[:ind]).any())
+
+    def test_allnans(self):
+        mat = np.array([np.nan]*9).reshape(3, 3)
+        tgt = np.iinfo(np.intp).min
+        for f in self.nanfuncs:
+            for axis in [None, 0, 1]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    res = f(mat, axis=axis)
+                    assert_((res == tgt).all())
+                    assert_(len(w) == 1)
+                    assert_(issubclass(w[0].category, NanWarning))
+
+    def test_empty(self):
+        mat = np.zeros((0,3))
+        for f in self.nanfuncs:
+            for axis in [0, None]:
+                assert_raises(ValueError, f, mat, axis=axis)
+            for axis in [1]:
+                res = f(mat, axis=axis)
+                assert_equal(res, np.zeros(0))
+
+
+class TestNanFunctions_IntTypes(TestCase):
+
+    int_types = (
+            np.int8, np.int16, np.int32, np.int64, np.uint8,
+            np.uint16, np.uint32, np.uint64)
+
+    def setUp(self, *args, **kwargs):
+        self.mat = np.array([127, 39,  93,  87, 46])
+
+    def integer_arrays(self):
+        for dtype in self.int_types:
+            yield self.mat.astype(dtype)
+
+    def test_nanmin(self):
+        min_value = min(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(nanmin(mat), min_value)
+
+    def test_nanmax(self):
+        max_value = max(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(nanmax(mat), max_value)
+
+    def test_nanargmin(self):
+        min_arg = np.argmin(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(nanargmin(mat), min_arg)
+
+    def test_nanargmax(self):
+        max_arg = np.argmax(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(nanargmax(mat), max_arg)
+
+
+class TestNanFunctions_Sum(TestCase):
+
+    def test_mutation(self):
+        # Check that passed array is not modified.
+        ndat = _ndat.copy()
+        nansum(ndat)
+        assert_equal(ndat, _ndat)
+
+    def test_keepdims(self):
+        mat = np.eye(3)
+        for axis in [None, 0, 1]:
+            tgt = np.sum(mat, axis=axis, keepdims=True)
+            res = nansum(mat, axis=axis, keepdims=True)
+            assert_(res.ndim == tgt.ndim)
+
+    def test_out(self):
+        mat = np.eye(3)
+        resout = np.zeros(3)
+        tgt = np.sum(mat, axis=1)
+        res = nansum(mat, axis=1, out=resout)
+        assert_almost_equal(res, resout)
+        assert_almost_equal(res, tgt)
+
+    def test_dtype_from_dtype(self):
+        mat = np.eye(3)
+        codes = 'efdgFDG'
+        for c in codes:
+            tgt = np.sum(mat, dtype=np.dtype(c), axis=1).dtype.type
+            res = nansum(mat, dtype=np.dtype(c), axis=1).dtype.type
+            assert_(res is tgt)
+            # scalar case
+            tgt = np.sum(mat, dtype=np.dtype(c), axis=None).dtype.type
+            res = nansum(mat, dtype=np.dtype(c), axis=None).dtype.type
+            assert_(res is tgt)
+
+    def test_dtype_from_char(self):
+        mat = np.eye(3)
+        codes = 'efdgFDG'
+        for c in codes:
+            tgt = np.sum(mat, dtype=c, axis=1).dtype.type
+            res = nansum(mat, dtype=c, axis=1).dtype.type
+            assert_(res is tgt)
+            # scalar case
+            tgt = np.sum(mat, dtype=c, axis=None).dtype.type
+            res = nansum(mat, dtype=c, axis=None).dtype.type
+            assert_(res is tgt)
+
+    def test_dtype_from_input(self):
+        codes = 'efdgFDG'
+        for c in codes:
+            mat = np.eye(3, dtype=c)
+            tgt = np.sum(mat, axis=1).dtype.type
+            res = nansum(mat, axis=1).dtype.type
+            assert_(res is tgt)
+            # scalar case
+            tgt = np.sum(mat, axis=None).dtype.type
+            res = nansum(mat, axis=None).dtype.type
+            assert_(res is tgt)
+
+    def test_result_values(self):
+            tgt = [np.sum(d) for d in _rdat]
+            res = nansum(_ndat, axis=1)
+            assert_almost_equal(res, tgt)
+
+    def test_allnans(self):
+        # Check for FutureWarning and later change of return from
+        # NaN to zero.
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            res = nansum([np.nan]*3, axis=None)
+            if np.__version__[:3] < '1.9':
+                assert_(np.isnan(res), 'result is not NaN')
+                assert_(len(w) == 1, 'no warning raised')
+                assert_(issubclass(w[0].category, FutureWarning))
+            else:
+                assert_(res == 0, 'result is not 0')
+                assert_(len(w) == 0, 'warning raised')
+
+    def test_empty(self):
+        mat = np.zeros((0,3))
+        if np.__version__[:3] < '1.9':
+            tgt = [np.nan]*3
+            res = nansum(mat, axis=0)
+            assert_equal(res, tgt)
+            tgt = []
+            res = nansum(mat, axis=1)
+            assert_equal(res, tgt)
+            tgt = np.nan
+            res = nansum(mat, axis=None)
+            assert_equal(res, tgt)
+        else:
+            tgt = [0]*3
+            res = nansum(mat, axis=0)
+            assert_equal(res, tgt)
+            tgt = []
+            res = nansum(mat, axis=1)
+            assert_equal(res, tgt)
+            tgt = 0
+            res = nansum(mat, axis=None)
+            assert_equal(res, tgt)
+
+
+class TestNanFunctions_MeanVarStd(TestCase):
+
+    nanfuncs = [nanmean, nanvar, nanstd]
+    stdfuncs = [np.mean, np.var, np.std]
+
+    def test_mutation(self):
+        # Check that passed array is not modified.
+        ndat = _ndat.copy()
+        for f in self.nanfuncs:
+            f(ndat)
+            assert_equal(ndat, _ndat)
+
+    def test_dtype_error(self):
+        for f in self.nanfuncs:
+            for dtype in [np.bool_, np.int_, np.object]:
+                assert_raises( TypeError, f, _ndat, axis=1, dtype=np.int)
+
+    def test_out_dtype_error(self):
+        for f in self.nanfuncs:
+            for dtype in [np.bool_, np.int_, np.object]:
+                out = np.empty(_ndat.shape[0], dtype=dtype)
+                assert_raises( TypeError, f, _ndat, axis=1, out=out)
+
+    def test_keepdims(self):
+        mat = np.eye(3)
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for axis in [None, 0, 1]:
+                tgt = rf(mat, axis=axis, keepdims=True)
+                res = nf(mat, axis=axis, keepdims=True)
+                assert_(res.ndim == tgt.ndim)
+
+    def test_out(self):
+        mat = np.eye(3)
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            resout = np.zeros(3)
+            tgt = rf(mat, axis=1)
+            res = nf(mat, axis=1, out=resout)
+            assert_almost_equal(res, resout)
+            assert_almost_equal(res, tgt)
+
+    def test_dtype_from_dtype(self):
+        mat = np.eye(3)
+        codes = 'efdgFDG'
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for c in codes:
+                tgt = rf(mat, dtype=np.dtype(c), axis=1).dtype.type
+                res = nf(mat, dtype=np.dtype(c), axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                tgt = rf(mat, dtype=np.dtype(c), axis=None).dtype.type
+                res = nf(mat, dtype=np.dtype(c), axis=None).dtype.type
+                assert_(res is tgt)
+
+    def test_dtype_from_char(self):
+        mat = np.eye(3)
+        codes = 'efdgFDG'
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for c in codes:
+                tgt = rf(mat, dtype=c, axis=1).dtype.type
+                res = nf(mat, dtype=c, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                tgt = rf(mat, dtype=c, axis=None).dtype.type
+                res = nf(mat, dtype=c, axis=None).dtype.type
+                assert_(res is tgt)
+
+    def test_dtype_from_input(self):
+        codes = 'efdgFDG'
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            for c in codes:
+                mat = np.eye(3, dtype=c)
+                tgt = rf(mat, axis=1).dtype.type
+                res = nf(mat, axis=1).dtype.type
+                assert_(res is tgt)
+                # scalar case
+                tgt = rf(mat, axis=None).dtype.type
+                res = nf(mat, axis=None).dtype.type
+                assert_(res is tgt)
+
+    def test_ddof(self):
+        nanfuncs = [nanvar, nanstd]
+        stdfuncs = [np.var, np.std]
+        for nf, rf in zip(nanfuncs, stdfuncs):
+            for ddof in [0, 1]:
+                tgt = [rf(d, ddof=ddof) for d in _rdat]
+                res = nf(_ndat, axis=1, ddof=ddof)
+                assert_almost_equal(res, tgt)
+
+    def test_ddof_too_big(self):
+        nanfuncs = [nanvar, nanstd]
+        stdfuncs = [np.var, np.std]
+        dsize = [len(d) for d in _rdat]
+        for nf, rf in zip(nanfuncs, stdfuncs):
+            for ddof in range(5):
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    tgt = [ddof >= d for d in dsize]
+                    res = nf(_ndat, axis=1, ddof=ddof)
+                    assert_equal(np.isnan(res), tgt)
+                    if any(tgt):
+                        assert_(len(w) == 1)
+                        assert_(issubclass(w[0].category, NanWarning))
+                    else:
+                        assert_(len(w) == 0)
+
+    def test_result_values(self):
+        for nf, rf in zip(self.nanfuncs, self.stdfuncs):
+            tgt = [rf(d) for d in _rdat]
+            res = nf(_ndat, axis=1)
+            assert_almost_equal(res, tgt)
+
+    def test_allnans(self):
+        mat = np.array([np.nan]*9).reshape(3, 3)
+        for f in self.nanfuncs:
+            for axis in [None, 0, 1]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_(np.isnan(f(mat, axis=axis)).all())
+                    assert_(len(w) == 1)
+                    assert_(issubclass(w[0].category, NanWarning))
+
+    def test_empty(self):
+        mat = np.zeros((0,3))
+        for f in self.nanfuncs:
+            for axis in [0, None]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_(np.isnan(f(mat, axis=axis)).all())
+                    assert_(len(w) == 1)
+                    assert_(issubclass(w[0].category, NanWarning))
+            for axis in [1]:
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_equal(f(mat, axis=axis), np.zeros([]))
+                    assert_(len(w) == 0)
+
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/testing/nosetester.py b/numpy/testing/nosetester.py
index 212afe79b..79ae79765 100644
--- a/numpy/testing/nosetester.py
+++ b/numpy/testing/nosetester.py
@@ -377,6 +377,7 @@ class NoseTester(object):
             warnings.filterwarnings("ignore", message="numpy.dtype size changed")
             warnings.filterwarnings("ignore", message="numpy.ufunc size changed")
             warnings.filterwarnings("ignore", category=ModuleDeprecationWarning)
+            warnings.filterwarnings("ignore", category=FutureWarning)
 
             from .noseclasses import NumpyTestProgram