@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+import functools
 import warnings
 
 import numpy as np
@@ -49,40 +50,29 @@ def _move_tuple_axes_first(array, axis):
     return array_new
 
 
-def _nanmean(array, axis=None):
-    """Bottleneck nanmean function that handle tuple axis."""
-    if isinstance(axis, tuple):
-        array = _move_tuple_axes_first(array, axis=axis)
-        axis = 0
-
-    if isinstance(array, Quantity):
-        return array.__array_wrap__(bottleneck.nanmean(array, axis=axis))
-    else:
-        return bottleneck.nanmean(array, axis=axis)
+def _apply_bottleneck(function, array, axis=None, **kwargs):
+    """Wrap bottleneck function to handle tuple axis.
 
-
-def _nanmedian(array, axis=None):
-    """Bottleneck nanmedian function that handle tuple axis."""
+    Also takes care to ensure the output is of the expected type,
+    i.e., a quantity, numpy array, or numpy scalar.
+    """
     if isinstance(axis, tuple):
         array = _move_tuple_axes_first(array, axis=axis)
         axis = 0
 
+    result = function(array, axis=axis, **kwargs)
     if isinstance(array, Quantity):
-        return array.__array_wrap__(bottleneck.nanmedian(array, axis=axis))
+        return array.__array_wrap__(result)
+    elif isinstance(result, float):
+        # For compatibility with numpy, always return a numpy scalar.
+        return np.float64(result)
     else:
-        return bottleneck.nanmedian(array, axis=axis)
+        return result
 
 
-def _nanstd(array, axis=None, ddof=0):
-    """Bottleneck nanstd function that handle tuple axis."""
-    if isinstance(axis, tuple):
-        array = _move_tuple_axes_first(array, axis=axis)
-        axis = 0
-
-    if isinstance(array, Quantity):
-        return array.__array_wrap__(bottleneck.nanstd(array, axis=axis, ddof=ddof))
-    else:
-        return bottleneck.nanstd(array, axis=axis, ddof=ddof)
+_nanmean = functools.partial(_apply_bottleneck, bottleneck.nanmean)
+_nanmedian = functools.partial(_apply_bottleneck, bottleneck.nanmedian)
+_nanstd = functools.partial(_apply_bottleneck, bottleneck.nanstd)
 
 
 def _nanmadstd(array, axis=None):
@@ -316,15 +306,7 @@ def _compute_bounds(self, data, axis=None):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=RuntimeWarning)
             cen = self._cenfunc_parsed(data, axis=axis)
-            if np.isscalar(cen) and not isinstance(cen, np.float64):
-                # change bottleneck float scalar to np.float64
-                cen = np.float64(cen)
-
             std = self._stdfunc_parsed(data, axis=axis)
-            if np.isscalar(std) and not isinstance(std, np.float64):
-                # change bottleneck float scalar to np.float64
-                std = np.float64(std)
-
             self._min_value = cen - (std * self.sigma_lower)
             self._max_value = cen + (std * self.sigma_upper)