diff --git a/src/histogram/histograms.rs b/src/histogram/histograms.rs
index 825aadb7..0a0529cd 100644
--- a/src/histogram/histograms.rs
+++ b/src/histogram/histograms.rs
@@ -78,6 +78,20 @@ impl<A: Ord> Histogram<A> {
     pub fn grid(&self) -> &Grid<A> {
         &self.grid
     }
+
+    /// Returns the cumulative distribution function of a histogram. 
+    /// Equivalent to the numpy histogram function cumsum
+    pub fn cumulative_sum(&self) -> ArrayD<usize> {
+        let mut cdf = self.counts.clone();
+        for i in 0..self.ndim() {
+            for j in 1..cdf.shape()[i] { 
+                let temp = cdf.index_axis(Axis(i), j - 1).to_owned();
+                let mut ax = cdf.index_axis_mut(Axis(i), j);
+                ax += &temp;
+            }
+        }
+        cdf
+    }
 }
 
 /// Extension trait for `ArrayBase` providing methods to compute histograms.
@@ -163,3 +177,39 @@ impl<A, S> HistogramExt<A, S> for ArrayBase<S, Ix2>
         histogram
     }
 }
+
+#[cfg(test)]
+mod auto_tests {
+    use super::*;
+    use histogram::histograms::HistogramExt;
+    use histogram::bins::{Bins, Edges};
+
+    #[test]
+    fn histogram_cdf_1d() {
+        let data = arr2(&[[1], [2], [3], [4], [1], [1], [4], [5], [8], [8], [8]]);
+        let grid = Grid::from(vec![
+                              Bins::new(
+                                  Edges::from(vec![0,1,2,3,4,5,6,7,8,9]))]);
+
+        let histogram = data.histogram(grid);
+                       //0, 1  2  3  4  5  6  7  8
+        let cdf = arr1(&[0, 3, 4, 5, 7, 8, 8, 8, 11]).into_dyn();
+        assert_eq!(histogram.cumulative_sum(), cdf);
+    }
+    
+    #[test]
+    fn histogram_cdf_2d() {
+        let data = arr2(&[[0, 2], [4, 4], [1, 1], [3, 3], [0, 2]]);
+        let grid = Grid::from(vec![
+                              Bins::new(Edges::from(vec![0, 1, 2, 3, 4])),
+                              Bins::new(Edges::from(vec![0, 1, 2, 3, 4]))]);
+
+        let histogram = data.histogram(grid);
+
+        let cdf = arr2(&[[0, 0, 2, 2],
+                         [0, 1, 3, 3],
+                         [0, 1, 3, 3],
+                         [0, 1, 3, 4]]).into_dyn();
+        assert_eq!(histogram.cumulative_sum(), cdf);
+    }
+}