diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs
index fe3e197..41ec415 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -168,16 +168,10 @@ pub trait BaseVector<T: RealNumber>: Clone + Debug {
 
     /// Computes the arithmetic mean.
     fn mean(&self) -> T {
-        let n = self.len();
-        let mut mean = T::zero();
-
-        for i in 0..n {
-            mean += self.get(i);
-        }
-        mean / T::from_usize(n).unwrap()
+        self.sum() / T::from_usize(self.len()).unwrap()
     }
-    /// Computes the standard deviation.
-    fn std(&self) -> T {
+    /// Computes variance.
+    fn var(&self) -> T {
         let n = self.len();
 
         let mut mu = T::zero();
@@ -189,7 +183,11 @@ pub trait BaseVector<T: RealNumber>: Clone + Debug {
             sum += xi * xi;
         }
         mu /= div;
-        (sum / div - mu * mu).sqrt()
+        sum / div - mu * mu
+    }
+    /// Computes the standard deviation.
+    fn std(&self) -> T {
+        self.var().sqrt()
     }
 }
 
@@ -592,4 +590,11 @@ mod tests {
 
         assert!((m.std() - 0.81f64).abs() < 1e-2);
     }
+
+    #[test]
+    fn var() {
+        let m = vec![1., 2., 3., 4.];
+
+        assert!((m.var() - 1.25f64).abs() < std::f64::EPSILON);
+    }
 }
diff --git a/src/linalg/stats.rs b/src/linalg/stats.rs
index ecb7ceb..fc339e0 100644
--- a/src/linalg/stats.rs
+++ b/src/linalg/stats.rs
@@ -35,8 +35,8 @@ pub trait MatrixStats<T: RealNumber>: BaseMatrix<T> {
         x
     }
 
-    /// Computes the standard deviation along the specified axis.
-    fn std(&self, axis: u8) -> Vec<T> {
+    /// Computes variance along the specified axis.
+    fn var(&self, axis: u8) -> Vec<T> {
         let (n, m) = match axis {
             0 => {
                 let (n, m) = self.shape();
@@ -61,7 +61,24 @@ pub trait MatrixStats<T: RealNumber>: BaseMatrix<T> {
                 sum += a * a;
             }
             mu /= div;
-            x[i] = (sum / div - mu * mu).sqrt();
+            x[i] = sum / div - mu * mu;
+        }
+
+        x
+    }
+
+    /// Computes the standard deviation along the specified axis.
+    fn std(&self, axis: u8) -> Vec<T> {
+
+        let mut x = self.var(axis);
+
+        let n = match axis {
+            0 => self.shape().1,
+            _ => self.shape().0,
+        };
+
+        for i in 0..n {            
+            x[i] = x[i].sqrt();
         }
 
         x
@@ -122,6 +139,19 @@ mod tests {
         assert!(m.std(1).approximate_eq(&expected_1, 1e-2));
     }
 
+    #[test]
+    fn var() {
+        let m = DenseMatrix::from_2d_array(&[
+            &[1., 2., 3., 4.],
+            &[5., 6., 7., 8.]
+        ]);
+        let expected_0 = vec![4., 4., 4., 4.];
+        let expected_1 = vec![1.25, 1.25];
+
+        assert!(m.var(0).approximate_eq(&expected_0, std::f64::EPSILON));
+        assert!(m.var(1).approximate_eq(&expected_1, std::f64::EPSILON));
+    }   
+
     #[test]
     fn scale() {
         let mut m = DenseMatrix::from_2d_array(&[&[1., 2., 3.], &[4., 5., 6.]]);