diff --git a/src/benchmark/measurement.rs b/src/benchmark/measurement.rs
index 718ca88..da35e9b 100644
--- a/src/benchmark/measurement.rs
+++ b/src/benchmark/measurement.rs
@@ -54,13 +54,14 @@ impl Measurements {
         mean(self.wall_clock_times())
     }
 
-    /// The standard deviation of wall clock times, if at least two runs were measured.
+    /// The standard deviation of all wall clock times. Not available if only one run has been performed
     pub fn stddev(&self) -> Option<Time> {
-        if self.len() < 2 {
+        let times: Vec<_> = self.wall_clock_times().collect(); // TODO: Avoid collecting
+
+        if times.len() < 2 {
             None
         } else {
-            let times: Vec<_> = self.wall_clock_times().collect();
-            Some(standard_deviation(&times, self.time_wall_clock_mean()))
+            Some(standard_deviation(times))
         }
     }
 
diff --git a/src/quantity/statistics.rs b/src/quantity/statistics.rs
index 8f26c6e..9434b41 100644
--- a/src/quantity/statistics.rs
+++ b/src/quantity/statistics.rs
@@ -1,8 +1,12 @@
-//! Statistics for non-empty samples of finite measurements.
+use std::ops::Add;
+use std::ops::AddAssign;
+use std::ops::Div;
 
-use std::ops::{Add, AddAssign, Div};
+use uom::num_traits;
+use uom::si::f64::Ratio;
+use uom::si::ratio::ratio;
 
-use crate::quantity::{byte, ratio, second, Information, Ratio, Time, Zero};
+use super::{byte, second, Information, Time};
 
 /// A min function that assumes no NaNs and at least one element
 pub fn min<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
@@ -20,14 +24,14 @@ pub fn max<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
         .expect("'max' requires at least one element")
 }
 
-/// Arithmetic mean, summing measurements in their original order.
+/// A mean function that assumes at least one element
 pub fn mean<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
 where
-    Q: AddAssign + Zero + Div<Ratio, Output = P>,
+    Q: AddAssign + num_traits::Zero + Div<Ratio, Output = P>,
     P: Into<Q>,
 {
     let mut sum = Q::zero();
-    let mut count = 0usize;
+    let mut count = 0;
     for value in values {
         sum += value;
         count += 1;
@@ -37,25 +41,20 @@ where
     (sum / count).into()
 }
 
-/// Median without modifying the input sample.
 pub fn median<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
 where
     Q: Copy + PartialOrd + Add<Output = Q> + Div<Ratio, Output = P>,
     P: Into<Q>,
 {
     let mut values = values.into_iter().collect::<Vec<_>>();
-    assert!(
-        !values.is_empty(),
-        "median requires at least one measurement"
-    );
-    values.sort_unstable_by(|a, b| a.partial_cmp(b).expect("No NaN measurements"));
+    values.sort_by(|a, b| a.partial_cmp(b).expect("No NaN values"));
 
     let len = values.len();
     if len % 2 == 0 {
         let mid = len / 2;
         let a = &values[mid - 1];
         let b = &values[mid];
-        ((*b + *a) / Ratio::new::<ratio>(2.)).into()
+        ((*a + *b) / Ratio::new::<ratio>(2.)).into()
     } else {
         values[len / 2]
     }
@@ -88,22 +87,24 @@ impl UnsafeRawValue for Information {
     }
 }
 
-/// Sample standard deviation, using an already computed mean.
-pub fn standard_deviation<Q: UnsafeRawValue>(values: &[Q], mean: Q) -> Q {
-    let mean = mean.unsafe_raw_value();
-    assert!(
-        values.len() > 1,
-        "standard deviation requires at least two measurements"
-    );
-    let sum = values
-        .iter()
-        .map(|value| (value.unsafe_raw_value() - mean) * (value.unsafe_raw_value() - mean))
-        .fold(0.0, |sum, deviation| sum + deviation);
-    assert!(sum >= 0.0, "invalid sum of squared deviations");
+fn standard_deviation_f64(values: &[f64]) -> f64 {
+    let mean_value = mean(values.iter().copied());
 
-    // Preserve the operation order of statistical 1.0: dividing by n - 1
-    // before taking the square root also preserves exported floating-point values.
-    Q::unsafe_from_raw_value((sum / (values.len() as f64 - 1.0)).sqrt())
+    let mut squared_deviations = 0.;
+    let mut n = 0;
+    for value in values {
+        let deviation = value - mean_value;
+        squared_deviations += deviation * deviation;
+        n += 1;
+    }
+
+    (1. / ((n - 1) as f64) * squared_deviations).sqrt()
+}
+
+pub fn standard_deviation<Q: UnsafeRawValue>(values: impl IntoIterator<Item = Q> + Clone) -> Q {
+    let values: Vec<_> = values.into_iter().map(|q| q.unsafe_raw_value()).collect();
+    let result = standard_deviation_f64(&values);
+    Q::unsafe_from_raw_value(result)
 }
 
 /// Compute modified Z-scores for a given sample. A (unmodified) Z-score is defined by
@@ -219,7 +220,7 @@ mod tests {
             Time::new::<millisecond>(200.0),
             Time::new::<microsecond>(300_000.0),
         ];
-        let result = standard_deviation(&values, mean(values));
+        let result = standard_deviation(values);
         assert_relative_eq!(result.get::<millisecond>(), 100.0);
     }
 }