Merge pull request #3866 from charris/refactor-1.9-nanfunctions

Refactor 1.9 nanfunctions

2 files changed, 343 insertions, 210 deletions

diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py
index 44cc2b163..5766084ab 100644
--- a/numpy/lib/nanfunctions.py
+++ b/numpy/lib/nanfunctions.py
@@ -13,23 +13,18 @@ Functions
 - `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'
+    'nanvar', 'nanstd'
     ]
 
-class NanWarning(RuntimeWarning): pass
-
 
 def _replace_nan(a, val):
     """
@@ -103,7 +98,8 @@ 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.
+    division is in place. Note that this is only called with `a` an inexact
+    type.
 
     Parameters
     ----------
@@ -140,37 +136,38 @@ def _divide_by_count(a, b, out=None):
 
 def nanmin(a, axis=None, out=None, keepdims=False):
     """
-    Return the minimum of an array or minimum along an axis, ignoring any
-    NaNs.
+    Return minimum of an array or minimum along an axis, ignoring any NaNs.
+    When all-NaN slices are encountered a ``RuntimeWarning`` is raised and
+    Nan is returned for that slice.
 
     Parameters
     ----------
     a : array_like
-        Array containing numbers whose minimum is desired. If `a` is not
-        an array, a conversion is attempted.
+        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.
+        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`.
+        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.
+        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
     --------
@@ -193,8 +190,8 @@ def nanmin(a, axis=None, out=None, keepdims=False):
     -----
     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.
+    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.
 
@@ -216,32 +213,51 @@ def nanmin(a, axis=None, out=None, keepdims=False):
     -inf
 
     """
-    return np.fmin.reduce(a, axis=axis, out=out, keepdims=keepdims)
+    if not isinstance(a, np.ndarray) or type(a) is np.ndarray:
+        # Fast, but not safe for subclasses of ndarray
+        res = np.fmin.reduce(a, axis=axis, out=out, keepdims=keepdims)
+        if np.isnan(res).any():
+            warnings.warn("All-NaN axis encountered", RuntimeWarning)
+    else:
+        # Slow, but safe for subclasses of ndarray
+        a, mask = _replace_nan(a, +np.inf)
+        res = np.amin(a, axis=axis, out=out, keepdims=keepdims)
+        if mask is None:
+            return res
+
+        # Check for all-NaN axis
+        mask = np.all(mask, axis=axis, keepdims=keepdims)
+        if np.any(mask):
+            res = _copyto(res, mask, np.nan)
+            warnings.warn("All-NaN axis encountered", RuntimeWarning)
+    return res
 
 
 def nanmax(a, axis=None, out=None, keepdims=False):
     """
-    Return the maximum of an array or maximum along an axis, ignoring any NaNs.
+    Return the maximum of an array or maximum along an axis, ignoring any
+    NaNs.  When all-NaN slices are encountered a ``RuntimeWarning`` is
+    raised and NaN is returned for that slice.
 
     Parameters
     ----------
     a : array_like
-        Array containing numbers whose maximum is desired. If `a` is not
-        an array, a conversion is attempted.
+        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.
+        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`.
+        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
 
@@ -273,8 +289,8 @@ def nanmax(a, axis=None, out=None, keepdims=False):
     -----
     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.
+    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.
 
@@ -296,15 +312,31 @@ def nanmax(a, axis=None, out=None, keepdims=False):
     inf
 
     """
-    return np.fmax.reduce(a, axis=axis, out=out, keepdims=keepdims)
+    if not isinstance(a, np.ndarray) or type(a) is np.ndarray:
+        # Fast, but not safe for subclasses of ndarray
+        res = np.fmax.reduce(a, axis=axis, out=out, keepdims=keepdims)
+        if np.isnan(res).any():
+            warnings.warn("All-NaN slice encountered", RuntimeWarning)
+    else:
+        # Slow, but safe for subclasses of ndarray
+        a, mask = _replace_nan(a, -np.inf)
+        res = np.amax(a, axis=axis, out=out, keepdims=keepdims)
+        if mask is None:
+            return res
+
+        # Check for all-NaN axis
+        mask = np.all(mask, axis=axis, keepdims=keepdims)
+        if np.any(mask):
+            res = _copyto(res, mask, np.nan)
+            warnings.warn("All-NaN axis encountered", RuntimeWarning)
+    return res
 
 
 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.
+    NaNs. For all-NaN slices ``ValueError`` is raised. Warning: the results
+    cannot be trusted if a slice contains only NaNs and Infs.
 
     Parameters
     ----------
@@ -336,23 +368,19 @@ def nanargmin(a, axis=None):
 
     """
     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
+    res = np.argmin(a, axis=axis)
+    if mask is not None:
+        mask = np.all(mask, axis=axis)
+        if np.any(mask):
+            raise ValueError("All-NaN slice encountered")
+    return res
 
 
 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.
+    NaNs. For all-NaN slices ``ValueError`` is raised. Warning: the
+    results cannot be trusted if a slice contains only NaNs and -Infs.
 
 
     Parameters
@@ -385,24 +413,21 @@ def nanargmax(a, axis=None):
 
     """
     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
+    res = np.argmax(a, axis=axis)
+    if mask is not None:
+        mask = np.all(mask, axis=axis)
+        if np.any(mask):
+            raise ValueError("All-NaN slice encountered")
+    return res
 
 
 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.
+    Return the sum of array elements over a given axis treating Not a
+    Numbers (NaNs) as zero.
 
-    In Numpy versions <= 1.8 Nan is returned for slices that
-    are all-NaN or empty. In later versions zero is returned.
+    In Numpy versions <= 1.8 Nan is returned for slices that are all-NaN or
+    empty. In later versions zero is returned.
 
     Parameters
     ----------
@@ -410,19 +435,23 @@ def nansum(a, axis=None, dtype=None, out=None, keepdims=0):
         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.
+        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.
+        The type of the returned array and of the accumulator in which the
+        elements are summed.  By default, the dtype of `a` is used.  An
+        exception is when `a` has an integer type with less precision than
+        the platform (u)intp. In that case, the default will be either
+        (u)int32 or (u)int64 depending on whether the platform is 32 or 64
+        bits. For inexact inputs, dtype 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.
+        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
@@ -444,17 +473,13 @@ def nansum(a, axis=None, dtype=None, out=None, keepdims=0):
 
     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.
+    If both positive and negative infinity are present, the sum will be Not
+    A Number (NaN).
 
+    Numpy integer arithmetic is modular. If the size of a sum exceeds the
+    size of an integer accumulator, its value will wrap around and the
+    result will be incorrect. Specifying ``dtype=double`` can alleviate
+    that problem.
 
     Examples
     --------
@@ -478,7 +503,7 @@ def nansum(a, axis=None, dtype=None, out=None, keepdims=0):
 
     """
     a, mask = _replace_nan(a, 0)
-    return a.sum(axis, dtype, out, keepdims)
+    return np.sum(a, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
 
 
 def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
@@ -489,7 +514,7 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
     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.
+    For all-NaN slices, NaN is returned and a `RuntimeWarning` is raised.
 
     .. versionadded:: 1.8.0
 
@@ -503,17 +528,17 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
         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.
+        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.
+        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`.
+        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
     -------
@@ -533,11 +558,11 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
     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.
+    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
     --------
@@ -552,7 +577,7 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
     """
     arr, mask = _replace_nan(a, 0)
     if mask is None:
-        return np.mean(arr, axis, dtype=dtype, out=out, keepdims=keepdims)
+        return np.mean(arr, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
 
     if dtype is not None:
         dtype = np.dtype(dtype)
@@ -564,28 +589,28 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False):
     # 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)
+        cnt = np.sum(~mask, axis=axis, dtype=np.intp, keepdims=keepdims)
+        tot = np.sum(arr, axis=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)
+        warnings.warn("Mean of empty slice", RuntimeWarning)
         # 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):
+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
+    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.
+    For all-NaN slices or slices with zero degrees of freedom, NaN is
+    returned and a `RuntimeWarning` is raised.
 
     .. versionadded:: 1.8.0
 
@@ -638,9 +663,9 @@ def nanvar(a, axis=None, dtype=None, out=None, ddof=0,
     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.
+    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.
@@ -664,7 +689,8 @@ def nanvar(a, axis=None, dtype=None, out=None, ddof=0,
     """
     arr, mask = _replace_nan(a, 0)
     if mask is None:
-        return np.var(arr, axis, dtype=dtype, out=out, keepdims=keepdims)
+        return np.var(arr, axis=axis, dtype=dtype, out=out, ddof=ddof,
+                      keepdims=keepdims)
 
     if dtype is not None:
         dtype = np.dtype(dtype)
@@ -677,30 +703,29 @@ def nanvar(a, axis=None, dtype=None, out=None, ddof=0,
         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)
+        cnt = np.sum(~mask, axis=axis, dtype=np.intp, keepdims=True)
+        avg = np.sum(arr, axis=axis, dtype=dtype, keepdims=True)
+        avg = _divide_by_count(avg, cnt)
 
         # Compute squared deviation from mean.
-        x = arr - avg
-        np.copyto(x, 0, where=mask)
+        arr -= avg
+        arr = _copyto(arr, 0, mask)
         if issubclass(arr.dtype.type, np.complexfloating):
-            sqr = np.multiply(x, x.conj(), out=x).real
+            sqr = np.multiply(arr, arr.conj(), out=arr).real
         else:
-            sqr = np.multiply(x, x, out=x)
-
-        # adjust cnt.
-        if not keepdims:
-            cnt = cnt.squeeze(axis)
-        cnt -= ddof
+            sqr = np.multiply(arr, arr, out=arr)
 
         # Compute variance.
-        var = np.add.reduce(sqr, axis, dtype=dtype, out=out, keepdims=keepdims)
-        var = _divide_by_count(var, cnt)
+        var = np.sum(sqr, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
+        if var.ndim < cnt.ndim:
+            # Subclasses of ndarray may ignore keepdims, so check here.
+            cnt = cnt.squeeze(axis)
+        dof = cnt - ddof
+        var = _divide_by_count(var, dof)
 
-    isbad = (cnt <= 0)
-    if isbad.any():
-        warnings.warn("Degrees of freedom <= 0 for slice.", NanWarning)
+    isbad = (dof <= 0)
+    if np.any(isbad):
+        warnings.warn("Degrees of freedom <= 0 for slice.", RuntimeWarning)
         # NaN, inf, or negative numbers are all possible bad
         # values, so explicitly replace them with NaN.
         var = _copyto(var, np.nan, isbad)
@@ -712,11 +737,13 @@ 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.
+    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.
+    For all-NaN slices or slices with zero degrees of freedom, NaN is
+    returned and a `RuntimeWarning` is raised.
 
     .. versionadded:: 1.8.0
 
@@ -729,12 +756,12 @@ def nanstd(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
         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.
+        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.
+        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
@@ -761,26 +788,26 @@ def nanstd(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     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))``.
+    deviations from the mean: ``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
+    ``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.
+    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.
+    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
     --------
@@ -793,7 +820,8 @@ def nanstd(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     array([ 0.,  0.5])
 
     """
-    var = nanvar(a, axis, dtype, out, ddof, keepdims)
+    var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
+                 keepdims=keepdims)
     if isinstance(var, np.ndarray):
         std = np.sqrt(var, out=var)
     else:
diff --git a/numpy/lib/tests/test_nanfunctions.py b/numpy/lib/tests/test_nanfunctions.py
index 70e7865db..af01a7167 100644
--- a/numpy/lib/tests/test_nanfunctions.py
+++ b/numpy/lib/tests/test_nanfunctions.py
@@ -7,31 +7,25 @@ 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]]
-)
+# Test data
+_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.1610, 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])
-]
+# 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]
+    nanfuncs = [np.nanmin, np.nanmax]
     stdfuncs = [np.min, np.max]
 
     def test_mutation(self):
@@ -81,22 +75,50 @@ class TestNanFunctions_MinMax(TestCase):
         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())
+                with warnings.catch_warnings(record=True) as w:
+                    warnings.simplefilter('always')
+                    assert_(np.isnan(f(mat, axis=axis)).all())
+                    assert_(len(w) == 1, 'no warning raised')
+                    assert_(issubclass(w[0].category, RuntimeWarning))
+            # Check scalars
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter('always')
+                assert_(np.isnan(f(np.nan)))
+                assert_(len(w) == 1, 'no warning raised')
+                assert_(issubclass(w[0].category, RuntimeWarning))
 
     def test_masked(self):
         mat = np.ma.fix_invalid(_ndat)
         msk = mat._mask.copy()
-        for f in [nanmin]:
+        for f in [np.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())
 
+    def test_scalar(self):
+        for f in self.nanfuncs:
+            assert_(f(0.) == 0.)
+
+    def test_matrices(self):
+        # Check that it works and that type and
+        # shape are preserved
+        mat = np.matrix(np.eye(3))
+        for f in self.nanfuncs:
+            res = f(mat, axis=0)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (1, 3))
+            res = f(mat, axis=1)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (3, 1))
+            res = f(mat)
+            assert_(np.isscalar(res))
+
 
 class TestNanFunctions_ArgminArgmax(TestCase):
 
-    nanfuncs = [nanargmin, nanargmax]
+    nanfuncs = [np.nanargmin, np.nanargmax]
 
     def test_mutation(self):
         # Check that passed array is not modified.
@@ -120,15 +142,10 @@ class TestNanFunctions_ArgminArgmax(TestCase):
 
     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))
+                assert_raises(ValueError, f, mat, axis=axis)
+            assert_raises(ValueError, f, np.nan)
 
     def test_empty(self):
         mat = np.zeros((0, 3))
@@ -139,39 +156,83 @@ class TestNanFunctions_ArgminArgmax(TestCase):
                 res = f(mat, axis=axis)
                 assert_equal(res, np.zeros(0))
 
+    def test_scalar(self):
+        for f in self.nanfuncs:
+            assert_(f(0.) == 0.)
+
+    def test_matrices(self):
+        # Check that it works and that type and
+        # shape are preserved
+        mat = np.matrix(np.eye(3))
+        for f in self.nanfuncs:
+            res = f(mat, axis=0)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (1, 3))
+            res = f(mat, axis=1)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (3, 1))
+            res = f(mat)
+            assert_(np.isscalar(res))
+
 
 class TestNanFunctions_IntTypes(TestCase):
 
-    int_types = (
-        np.int8, np.int16, np.int32, np.int64, np.uint8,
-        np.uint16, np.uint32, np.uint64)
+    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])
+    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)
+        tgt = np.min(self.mat)
         for mat in self.integer_arrays():
-            assert_equal(nanmin(mat), min_value)
+            assert_equal(np.nanmin(mat), tgt)
 
     def test_nanmax(self):
-        max_value = max(self.mat)
+        tgt = np.max(self.mat)
         for mat in self.integer_arrays():
-            assert_equal(nanmax(mat), max_value)
+            assert_equal(np.nanmax(mat), tgt)
 
     def test_nanargmin(self):
-        min_arg = np.argmin(self.mat)
+        tgt = np.argmin(self.mat)
         for mat in self.integer_arrays():
-            assert_equal(nanargmin(mat), min_arg)
+            assert_equal(np.nanargmin(mat), tgt)
 
     def test_nanargmax(self):
-        max_arg = np.argmax(self.mat)
+        tgt = np.argmax(self.mat)
         for mat in self.integer_arrays():
-            assert_equal(nanargmax(mat), max_arg)
+            assert_equal(np.nanargmax(mat), tgt)
+
+    def test_nansum(self):
+        tgt = np.sum(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(np.nansum(mat), tgt)
+
+    def test_nanmean(self):
+        tgt = np.mean(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(np.nanmean(mat), tgt)
+
+    def test_nanvar(self):
+        tgt = np.var(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(np.nanvar(mat), tgt)
+
+        tgt = np.var(mat, ddof=1)
+        for mat in self.integer_arrays():
+            assert_equal(np.nanvar(mat, ddof=1), tgt)
+
+    def test_nanstd(self):
+        tgt = np.std(self.mat)
+        for mat in self.integer_arrays():
+            assert_equal(np.nanstd(mat), tgt)
+
+        tgt = np.std(self.mat, ddof=1)
+        for mat in self.integer_arrays():
+            assert_equal(np.nanstd(mat, ddof=1), tgt)
 
 
 class TestNanFunctions_Sum(TestCase):
@@ -179,21 +240,21 @@ class TestNanFunctions_Sum(TestCase):
     def test_mutation(self):
         # Check that passed array is not modified.
         ndat = _ndat.copy()
-        nansum(ndat)
+        np.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)
+            res = np.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)
+        res = np.nansum(mat, axis=1, out=resout)
         assert_almost_equal(res, resout)
         assert_almost_equal(res, tgt)
 
@@ -202,11 +263,11 @@ class TestNanFunctions_Sum(TestCase):
         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
+            res = np.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
+            res = np.nansum(mat, dtype=np.dtype(c), axis=None).dtype.type
             assert_(res is tgt)
 
     def test_dtype_from_char(self):
@@ -214,11 +275,11 @@ class TestNanFunctions_Sum(TestCase):
         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
+            res = np.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
+            res = np.nansum(mat, dtype=c, axis=None).dtype.type
             assert_(res is tgt)
 
     def test_dtype_from_input(self):
@@ -226,43 +287,65 @@ class TestNanFunctions_Sum(TestCase):
         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
+            res = np.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
+            res = np.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)
+            res = np.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.
+        # Check for FutureWarning
         with warnings.catch_warnings(record=True) as w:
             warnings.simplefilter('always')
-            res = nansum([np.nan]*3, axis=None)
+            res = np.nansum([np.nan]*3, axis=None)
             assert_(res == 0, 'result is not 0')
             assert_(len(w) == 0, 'warning raised')
+            # Check scalar
+            res = np.nansum(np.nan)
+            assert_(res == 0, 'result is not 0')
+            assert_(len(w) == 0, 'warning raised')
+            # Check there is no warning for not all-nan
+            np.nansum([0]*3, axis=None)
+            assert_(len(w) == 0, 'unwanted warning raised')
 
     def test_empty(self):
         mat = np.zeros((0, 3))
         tgt = [0]*3
-        res = nansum(mat, axis=0)
+        res = np.nansum(mat, axis=0)
         assert_equal(res, tgt)
         tgt = []
-        res = nansum(mat, axis=1)
+        res = np.nansum(mat, axis=1)
         assert_equal(res, tgt)
         tgt = 0
-        res = nansum(mat, axis=None)
+        res = np.nansum(mat, axis=None)
         assert_equal(res, tgt)
 
+    def test_scalar(self):
+        assert_(np.nansum(0.) == 0.)
+
+    def test_matrices(self):
+        # Check that it works and that type and
+        # shape are preserved
+        mat = np.matrix(np.eye(3))
+        res = np.nansum(mat, axis=0)
+        assert_(isinstance(res, np.matrix))
+        assert_(res.shape == (1, 3))
+        res = np.nansum(mat, axis=1)
+        assert_(isinstance(res, np.matrix))
+        assert_(res.shape == (3, 1))
+        res = np.nansum(mat)
+        assert_(np.isscalar(res))
+
 
 class TestNanFunctions_MeanVarStd(TestCase):
 
-    nanfuncs = [nanmean, nanvar, nanstd]
+    nanfuncs = [np.nanmean, np.nanvar, np.nanstd]
     stdfuncs = [np.mean, np.var, np.std]
 
     def test_mutation(self):
@@ -275,13 +358,13 @@ class TestNanFunctions_MeanVarStd(TestCase):
     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)
+                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)
+                assert_raises( TypeError, f, _ndat, axis=1, out=out)
 
     def test_keepdims(self):
         mat = np.eye(3)
@@ -333,14 +416,14 @@ class TestNanFunctions_MeanVarStd(TestCase):
                 mat = np.eye(3, dtype=c)
                 tgt = rf(mat, axis=1).dtype.type
                 res = nf(mat, axis=1).dtype.type
-                assert_(res is tgt)
+                assert_(res is tgt, "res %s, tgt %s" % (res, 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]
+        nanfuncs = [np.nanvar, np.nanstd]
         stdfuncs = [np.var, np.std]
         for nf, rf in zip(nanfuncs, stdfuncs):
             for ddof in [0, 1]:
@@ -349,7 +432,7 @@ class TestNanFunctions_MeanVarStd(TestCase):
                 assert_almost_equal(res, tgt)
 
     def test_ddof_too_big(self):
-        nanfuncs = [nanvar, nanstd]
+        nanfuncs = [np.nanvar, np.nanstd]
         stdfuncs = [np.var, np.std]
         dsize = [len(d) for d in _rdat]
         for nf, rf in zip(nanfuncs, stdfuncs):
@@ -361,7 +444,7 @@ class TestNanFunctions_MeanVarStd(TestCase):
                     assert_equal(np.isnan(res), tgt)
                     if any(tgt):
                         assert_(len(w) == 1)
-                        assert_(issubclass(w[0].category, NanWarning))
+                        assert_(issubclass(w[0].category, RuntimeWarning))
                     else:
                         assert_(len(w) == 0)
 
@@ -379,7 +462,11 @@ class TestNanFunctions_MeanVarStd(TestCase):
                     warnings.simplefilter('always')
                     assert_(np.isnan(f(mat, axis=axis)).all())
                     assert_(len(w) == 1)
-                    assert_(issubclass(w[0].category, NanWarning))
+                    assert_(issubclass(w[0].category, RuntimeWarning))
+                    # Check scalar
+                    assert_(np.isnan(f(np.nan)))
+                    assert_(len(w) == 2)
+                    assert_(issubclass(w[0].category, RuntimeWarning))
 
     def test_empty(self):
         mat = np.zeros((0, 3))
@@ -389,13 +476,31 @@ class TestNanFunctions_MeanVarStd(TestCase):
                     warnings.simplefilter('always')
                     assert_(np.isnan(f(mat, axis=axis)).all())
                     assert_(len(w) == 1)
-                    assert_(issubclass(w[0].category, NanWarning))
+                    assert_(issubclass(w[0].category, RuntimeWarning))
             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)
 
+    def test_scalar(self):
+        for f in self.nanfuncs:
+            assert_(f(0.) == 0.)
+
+    def test_matrices(self):
+        # Check that it works and that type and
+        # shape are preserved
+        mat = np.matrix(np.eye(3))
+        for f in self.nanfuncs:
+            res = f(mat, axis=0)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (1, 3))
+            res = f(mat, axis=1)
+            assert_(isinstance(res, np.matrix))
+            assert_(res.shape == (3, 1))
+            res = f(mat)
+            assert_(np.isscalar(res))
+
 
 if __name__ == "__main__":
     run_module_suite()