diff --git a/src/histogram/strategies.rs b/src/histogram/strategies.rs
index d25eb2c3..920936e2 100644
--- a/src/histogram/strategies.rs
+++ b/src/histogram/strategies.rs
@@ -181,15 +181,15 @@ impl<T> BinsBuildingStrategy for Sqrt<T>
 {
     type Elem = T;
 
-    /// **Panics** if the array is constant or if `a.len()==0` and division by 0 panics for `T`.
+    /// **Panics** if the array is constant or if `a.len()==0`.
     fn from_array<S>(a: &ArrayBase<S, Ix1>) -> Self
     where
         S: Data<Elem=Self::Elem>
     {
         let n_elems = a.len();
         let n_bins = (n_elems as f64).sqrt().round() as usize;
-        let min = a.min().clone();
-        let max = a.max().clone();
+        let min = a.min().unwrap().clone();
+        let max = a.max().unwrap().clone();
         let bin_width = compute_bin_width(min.clone(), max.clone(), n_bins);
         let builder = EquiSpaced::new(bin_width, min, max);
         Self { builder }
@@ -220,15 +220,15 @@ impl<T> BinsBuildingStrategy for Rice<T>
 {
     type Elem = T;
 
-    /// **Panics** if the array is constant or if `a.len()==0` and division by 0 panics for `T`.
+    /// **Panics** if the array is constant or if `a.len()==0`.
     fn from_array<S>(a: &ArrayBase<S, Ix1>) -> Self
     where
         S: Data<Elem=Self::Elem>
     {
         let n_elems = a.len();
         let n_bins = (2. * (n_elems as f64).powf(1./3.)).round() as usize;
-        let min = a.min().clone();
-        let max = a.max().clone();
+        let min = a.min().unwrap().clone();
+        let max = a.max().unwrap().clone();
         let bin_width = compute_bin_width(min.clone(), max.clone(), n_bins);
         let builder = EquiSpaced::new(bin_width, min, max);
         Self { builder }
@@ -259,15 +259,15 @@ impl<T> BinsBuildingStrategy for Sturges<T>
 {
     type Elem = T;
 
-    /// **Panics** if the array is constant or if `a.len()==0` and division by 0 panics for `T`.
+    /// **Panics** if the array is constant or if `a.len()==0`.
     fn from_array<S>(a: &ArrayBase<S, Ix1>) -> Self
     where
         S: Data<Elem=Self::Elem>
     {
         let n_elems = a.len();
         let n_bins = (n_elems as f64).log2().round() as usize + 1;
-        let min = a.min().clone();
-        let max = a.max().clone();
+        let min = a.min().unwrap().clone();
+        let max = a.max().unwrap().clone();
         let bin_width = compute_bin_width(min.clone(), max.clone(), n_bins);
         let builder = EquiSpaced::new(bin_width, min, max);
         Self { builder }
@@ -298,7 +298,7 @@ impl<T> BinsBuildingStrategy for FreedmanDiaconis<T>
 {
     type Elem = T;
 
-    /// **Panics** if `IQR==0` or if `a.len()==0` and division by 0 panics for `T`.
+    /// **Panics** if `IQR==0` or if `a.len()==0`.
     fn from_array<S>(a: &ArrayBase<S, Ix1>) -> Self
     where
         S: Data<Elem=Self::Elem>
@@ -306,13 +306,13 @@ impl<T> BinsBuildingStrategy for FreedmanDiaconis<T>
         let n_points = a.len();
 
         let mut a_copy = a.to_owned();
-        let first_quartile = a_copy.quantile_mut::<Nearest>(0.25);
-        let third_quartile = a_copy.quantile_mut::<Nearest>(0.75);
+        let first_quartile = a_copy.quantile_mut::<Nearest>(0.25).unwrap();
+        let third_quartile = a_copy.quantile_mut::<Nearest>(0.75).unwrap();
         let iqr = third_quartile - first_quartile;
 
         let bin_width = FreedmanDiaconis::compute_bin_width(n_points, iqr);
-        let min = a_copy.min().clone();
-        let max = a_copy.max().clone();
+        let min = a_copy.min().unwrap().clone();
+        let max = a_copy.max().unwrap().clone();
         let builder = EquiSpaced::new(bin_width, min, max);
         Self { builder }
     }
@@ -349,8 +349,7 @@ impl<T> BinsBuildingStrategy for Auto<T>
 {
     type Elem = T;
 
-    /// **Panics** if `IQR==0`, the array is constant or if
-    /// `a.len()==0` and division by 0 panics for `T`.
+    /// **Panics** if `IQR==0`, the array is constant, or `a.len()==0`.
     fn from_array<S>(a: &ArrayBase<S, Ix1>) -> Self
     where
         S: Data<Elem=Self::Elem>
@@ -403,7 +402,7 @@ impl<T> Auto<T>
 ///
 /// `bin_width = (max - min)/n`
 ///
-/// **Panics** if division by 0 panics for `T`.
+/// **Panics** if `n_bins == 0` and division by 0 panics for `T`.
 fn compute_bin_width<T>(min: T, max: T, n_bins: usize) -> T
 where
     T: Ord + Clone + FromPrimitive + NumOps + Zero,
diff --git a/src/quantile.rs b/src/quantile.rs
index 3e179539..ff7c310f 100644
--- a/src/quantile.rs
+++ b/src/quantile.rs
@@ -175,13 +175,6 @@ where
     S: Data<Elem = A>,
     D: Dimension,
 {
-    /// Finds the elementwise minimum of the array.
-    ///
-    /// **Panics** if the array is empty.
-    fn min(&self) -> &A
-    where
-        A: Ord;
-
     /// Finds the elementwise minimum of the array.
     ///
     /// Returns `None` if any of the pairwise orderings tested by the function
@@ -189,7 +182,7 @@ where
     /// floating-point NaN values in the array.)
     ///
     /// Additionally, returns `None` if the array is empty.
-    fn min_partialord(&self) -> Option<&A>
+    fn min(&self) -> Option<&A>
     where
         A: PartialOrd;
 
@@ -203,13 +196,6 @@ where
         A: MaybeNan,
         A::NotNan: Ord;
 
-    /// Finds the elementwise maximum of the array.
-    ///
-    /// **Panics** if the array is empty.
-    fn max(&self) -> &A
-    where
-        A: Ord;
-
     /// Finds the elementwise maximum of the array.
     ///
     /// Returns `None` if any of the pairwise orderings tested by the function
@@ -217,7 +203,7 @@ where
     /// floating-point NaN values in the array.)
     ///
     /// Additionally, returns `None` if the array is empty.
-    fn max_partialord(&self) -> Option<&A>
+    fn max(&self) -> Option<&A>
     where
         A: PartialOrd;
 
@@ -259,8 +245,8 @@ where
     /// - worst case: O(`m`^2);
     /// where `m` is the number of elements in the array.
     ///
-    /// **Panics** if `axis` is out of bounds or if `q` is not between
-    /// `0.` and `1.` (inclusive).
+    /// **Panics** if `axis` is out of bounds, if the axis has length 0, or if
+    /// `q` is not between `0.` and `1.` (inclusive).
     fn quantile_axis_mut<I>(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
     where
         D: RemoveAxis,
@@ -285,22 +271,11 @@ where
     S: Data<Elem = A>,
     D: Dimension,
 {
-    fn min(&self) -> &A
-    where
-        A: Ord,
-    {
-        let first = self
-            .iter()
-            .next()
-            .expect("Attempted to find min of empty array.");
-        self.fold(first, |acc, elem| if elem < acc { elem } else { acc })
-    }
-
-    fn min_partialord(&self) -> Option<&A>
+    fn min(&self) -> Option<&A>
     where
         A: PartialOrd,
     {
-        let first = self.iter().next()?;
+        let first = self.first()?;
         self.fold(Some(first), |acc, elem| match elem.partial_cmp(acc?)? {
             cmp::Ordering::Less => Some(elem),
             _ => acc,
@@ -312,7 +287,7 @@ where
         A: MaybeNan,
         A::NotNan: Ord,
     {
-        let first = self.iter().next().and_then(|v| v.try_as_not_nan());
+        let first = self.first().and_then(|v| v.try_as_not_nan());
         A::from_not_nan_ref_opt(self.fold_skipnan(first, |acc, elem| {
             Some(match acc {
                 Some(acc) => acc.min(elem),
@@ -321,22 +296,11 @@ where
         }))
     }
 
-    fn max(&self) -> &A
-    where
-        A: Ord,
-    {
-        let first = self
-            .iter()
-            .next()
-            .expect("Attempted to find max of empty array.");
-        self.fold(first, |acc, elem| if elem > acc { elem } else { acc })
-    }
-
-    fn max_partialord(&self) -> Option<&A>
+    fn max(&self) -> Option<&A>
     where
         A: PartialOrd,
     {
-        let first = self.iter().next()?;
+        let first = self.first()?;
         self.fold(Some(first), |acc, elem| match elem.partial_cmp(acc?)? {
             cmp::Ordering::Greater => Some(elem),
             _ => acc,
@@ -348,7 +312,7 @@ where
         A: MaybeNan,
         A::NotNan: Ord,
     {
-        let first = self.iter().next().and_then(|v| v.try_as_not_nan());
+        let first = self.first().and_then(|v| v.try_as_not_nan());
         A::from_not_nan_ref_opt(self.fold_skipnan(first, |acc, elem| {
             Some(match acc {
                 Some(acc) => acc.max(elem),
@@ -442,8 +406,10 @@ pub trait QuantileExt1d<A, S>
     /// - worst case: O(`m`^2);
     /// where `m` is the number of elements in the array.
     ///
+    /// Returns `None` if the array is empty.
+    ///
     /// **Panics** if `q` is not between `0.` and `1.` (inclusive).
-    fn quantile_mut<I>(&mut self, q: f64) -> A
+    fn quantile_mut<I>(&mut self, q: f64) -> Option<A>
     where
         A: Ord + Clone,
         S: DataMut,
@@ -454,13 +420,17 @@ impl<A, S> QuantileExt1d<A, S> for ArrayBase<S, Ix1>
     where
         S: Data<Elem = A>,
 {
-    fn quantile_mut<I>(&mut self, q: f64) -> A
+    fn quantile_mut<I>(&mut self, q: f64) -> Option<A>
     where
         A: Ord + Clone,
         S: DataMut,
         I: Interpolate<A>,
     {
-        self.quantile_axis_mut::<I>(Axis(0), q).into_scalar()
+        if self.is_empty() {
+            None
+        } else {
+            Some(self.quantile_axis_mut::<I>(Axis(0), q).into_scalar())
+        }
     }
 }
 
diff --git a/tests/quantile.rs b/tests/quantile.rs
index 7581e5fe..faaffd6f 100644
--- a/tests/quantile.rs
+++ b/tests/quantile.rs
@@ -11,16 +11,13 @@ use ndarray_stats::{
 #[test]
 fn test_min() {
     let a = array![[1, 5, 3], [2, 0, 6]];
-    assert_eq!(a.min(), &0);
-}
+    assert_eq!(a.min(), Some(&0));
 
-#[test]
-fn test_min_partialord() {
     let a = array![[1., 5., 3.], [2., 0., 6.]];
-    assert_eq!(a.min_partialord(), Some(&0.));
+    assert_eq!(a.min(), Some(&0.));
 
     let a = array![[1., 5., 3.], [2., ::std::f64::NAN, 6.]];
-    assert_eq!(a.min_partialord(), None);
+    assert_eq!(a.min(), None);
 }
 
 #[test]
@@ -41,16 +38,13 @@ fn test_min_skipnan_all_nan() {
 #[test]
 fn test_max() {
     let a = array![[1, 5, 7], [2, 0, 6]];
-    assert_eq!(a.max(), &7);
-}
+    assert_eq!(a.max(), Some(&7));
 
-#[test]
-fn test_max_partialord() {
     let a = array![[1., 5., 7.], [2., 0., 6.]];
-    assert_eq!(a.max_partialord(), Some(&7.));
+    assert_eq!(a.max(), Some(&7.));
 
     let a = array![[1., 5., 7.], [2., ::std::f64::NAN, 6.]];
-    assert_eq!(a.max_partialord(), None);
+    assert_eq!(a.max(), None);
 }
 
 #[test]
@@ -75,6 +69,20 @@ fn test_quantile_axis_mut_with_odd_axis_length() {
     assert!(p == a.subview(Axis(0), 1));
 }
 
+#[test]
+#[should_panic]
+fn test_quantile_axis_mut_with_zero_axis_length() {
+    let mut a = Array2::<i32>::zeros((5, 0));
+    a.quantile_axis_mut::<Lower>(Axis(1), 0.5);
+}
+
+#[test]
+fn test_quantile_axis_mut_with_empty_array() {
+    let mut a = Array2::<i32>::zeros((5, 0));
+    let p = a.quantile_axis_mut::<Lower>(Axis(0), 0.5);
+    assert_eq!(p.shape(), &[0]);
+}
+
 #[test]
 fn test_quantile_axis_mut_with_even_axis_length() {
     let mut b = arr2(&[[1, 3, 2, 10], [2, 4, 3, 11], [3, 5, 6, 12], [4, 6, 7, 13]]);