feat: + ridge regression

src/linalg/mod.rs

@@ -48,6 +48,7 @@ pub mod nalgebra_bindings;
 pub mod ndarray_bindings;
 /// QR factorization that factors a matrix into a product of an orthogonal matrix and an upper triangular matrix.
 pub mod qr;
+pub mod stats;
 /// Singular value decomposition.
 pub mod svd;
 
@@ -60,6 +61,7 @@ use cholesky::CholeskyDecomposableMatrix;
 use evd::EVDDecomposableMatrix;
 use lu::LUDecomposableMatrix;
 use qr::QRDecomposableMatrix;
+use stats::MatrixStats;
 use svd::SVDDecomposableMatrix;
 
 /// Column or row vector
@@ -163,6 +165,32 @@ pub trait BaseVector<T: RealNumber>: Clone + Debug {
     ///assert_eq!(a.unique(), vec![-7., -6., -2., 1., 2., 3., 4.]);
     /// ```
     fn unique(&self) -> Vec<T>;
+
+    /// Compute 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()
+    }
+    /// Compute the standard deviation.
+    fn std(&self) -> T {
+        let n = self.len();
+
+        let mut mu = T::zero();
+        let mut sum = T::zero();
+        let div = T::from_usize(n).unwrap();
+        for i in 0..n {
+            let xi = self.get(i);
+            mu += xi;
+            sum += xi * xi;
+        }
+        mu /= div;
+        (sum / div - mu * mu).sqrt()
+    }
 }
 
 /// Generic matrix type.
@@ -510,6 +538,7 @@ pub trait Matrix<T: RealNumber>:
     + QRDecomposableMatrix<T>
     + LUDecomposableMatrix<T>
     + CholeskyDecomposableMatrix<T>
+    + MatrixStats<T>
     + PartialEq
     + Display
 {
@@ -545,3 +574,22 @@ impl<'a, T: RealNumber, M: BaseMatrix<T>> Iterator for RowIter<'a, T, M> {
         res
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::linalg::BaseVector;
+
+    #[test]
+    fn mean() {
+        let m = vec![1., 2., 3.];
+
+        assert_eq!(m.mean(), 2.0);
+    }
+
+    #[test]
+    fn std() {
+        let m = vec![1., 2., 3.];
+
+        assert!((m.std() - 0.81f64).abs() < 1e-2);
+    }
+}

src/linalg/naive/dense_matrix.rs

@@ -12,6 +12,7 @@ use crate::linalg::cholesky::CholeskyDecomposableMatrix;
 use crate::linalg::evd::EVDDecomposableMatrix;
 use crate::linalg::lu::LUDecomposableMatrix;
 use crate::linalg::qr::QRDecomposableMatrix;
+use crate::linalg::stats::MatrixStats;
 use crate::linalg::svd::SVDDecomposableMatrix;
 use crate::linalg::Matrix;
 pub use crate::linalg::{BaseMatrix, BaseVector};
@@ -445,6 +446,8 @@ impl<T: RealNumber> LUDecomposableMatrix<T> for DenseMatrix<T> {}
 
 impl<T: RealNumber> CholeskyDecomposableMatrix<T> for DenseMatrix<T> {}
 
+impl<T: RealNumber> MatrixStats<T> for DenseMatrix<T> {}
+
 impl<T: RealNumber> Matrix<T> for DenseMatrix<T> {}
 
 impl<T: RealNumber> PartialEq for DenseMatrix<T> {

src/linalg/nalgebra_bindings.rs

@@ -46,6 +46,7 @@ use crate::linalg::cholesky::CholeskyDecomposableMatrix;
 use crate::linalg::evd::EVDDecomposableMatrix;
 use crate::linalg::lu::LUDecomposableMatrix;
 use crate::linalg::qr::QRDecomposableMatrix;
+use crate::linalg::stats::MatrixStats;
 use crate::linalg::svd::SVDDecomposableMatrix;
 use crate::linalg::Matrix as SmartCoreMatrix;
 use crate::linalg::{BaseMatrix, BaseVector};
@@ -550,6 +551,11 @@ impl<T: RealNumber + Scalar + AddAssign + SubAssign + MulAssign + DivAssign + Su
 {
 }
 
+impl<T: RealNumber + Scalar + AddAssign + SubAssign + MulAssign + DivAssign + Sum + 'static>
+    MatrixStats<T> for Matrix<T, Dynamic, Dynamic, VecStorage<T, Dynamic, Dynamic>>
+{
+}
+
 impl<T: RealNumber + Scalar + AddAssign + SubAssign + MulAssign + DivAssign + Sum + 'static>
     SmartCoreMatrix<T> for Matrix<T, Dynamic, Dynamic, VecStorage<T, Dynamic, Dynamic>>
 {

src/linalg/ndarray_bindings.rs

@@ -53,6 +53,7 @@ use crate::linalg::cholesky::CholeskyDecomposableMatrix;
 use crate::linalg::evd::EVDDecomposableMatrix;
 use crate::linalg::lu::LUDecomposableMatrix;
 use crate::linalg::qr::QRDecomposableMatrix;
+use crate::linalg::stats::MatrixStats;
 use crate::linalg::svd::SVDDecomposableMatrix;
 use crate::linalg::Matrix;
 use crate::linalg::{BaseMatrix, BaseVector};
@@ -500,6 +501,11 @@ impl<T: RealNumber + ScalarOperand + AddAssign + SubAssign + MulAssign + DivAssi
 {
 }
 
+impl<T: RealNumber + ScalarOperand + AddAssign + SubAssign + MulAssign + DivAssign + Sum>
+    MatrixStats<T> for ArrayBase<OwnedRepr<T>, Ix2>
+{
+}
+
 impl<T: RealNumber + ScalarOperand + AddAssign + SubAssign + MulAssign + DivAssign + Sum> Matrix<T>
     for ArrayBase<OwnedRepr<T>, Ix2>
 {

src/linalg/stats.rs

@@ -0,0 +1,139 @@
+//! # Various Statistical Methods
+//!
+//!
+
+use crate::linalg::BaseMatrix;
+use crate::math::num::RealNumber;
+
+/// Defines baseline implementations for various statistical functions
+pub trait MatrixStats<T: RealNumber>: BaseMatrix<T> {
+    /// Compute the arithmetic mean along the specified axis.
+    fn mean(&self, axis: u8) -> Vec<T> {
+        let (n, m) = match axis {
+            0 => {
+                let (n, m) = self.shape();
+                (m, n)
+            }
+            _ => self.shape(),
+        };
+
+        let mut x: Vec<T> = vec![T::zero(); n];
+
+        let div = T::from_usize(m).unwrap();
+
+        for i in 0..n {
+            for j in 0..m {
+                x[i] += match axis {
+                    0 => self.get(j, i),
+                    _ => self.get(i, j),
+                };
+            }
+            x[i] /= div;
+        }
+
+        x
+    }
+
+    /// Compute the standard deviation along the specified axis.
+    fn std(&self, axis: u8) -> Vec<T> {
+        let (n, m) = match axis {
+            0 => {
+                let (n, m) = self.shape();
+                (m, n)
+            }
+            _ => self.shape(),
+        };
+
+        let mut x: Vec<T> = vec![T::zero(); n];
+
+        let div = T::from_usize(m).unwrap();
+
+        for i in 0..n {
+            let mut mu = T::zero();
+            let mut sum = T::zero();
+            for j in 0..m {
+                let a = match axis {
+                    0 => self.get(j, i),
+                    _ => self.get(i, j),
+                };
+                mu += a;
+                sum += a * a;
+            }
+            mu /= div;
+            x[i] = (sum / div - mu * mu).sqrt();
+        }
+
+        x
+    }
+
+    /// standardize values by removing the mean and scaling to unit variance
+    fn scale_mut(&mut self, mean: &Vec<T>, std: &Vec<T>, axis: u8) {
+        let (n, m) = match axis {
+            0 => {
+                let (n, m) = self.shape();
+                (m, n)
+            }
+            _ => self.shape(),
+        };
+
+        for i in 0..n {
+            for j in 0..m {
+                match axis {
+                    0 => self.set(j, i, (self.get(j, i) - mean[i]) / std[i]),
+                    _ => self.set(i, j, (self.get(i, j) - mean[i]) / std[i]),
+                }
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::linalg::naive::dense_matrix::DenseMatrix;
+    use crate::linalg::BaseVector;
+
+    #[test]
+    fn mean() {
+        let m = DenseMatrix::from_2d_array(&[
+            &[1., 2., 3., 1., 2.],
+            &[4., 5., 6., 3., 4.],
+            &[7., 8., 9., 5., 6.],
+        ]);
+        let expected_0 = vec![4., 5., 6., 3., 4.];
+        let expected_1 = vec![1.8, 4.4, 7.];
+
+        assert_eq!(m.mean(0), expected_0);
+        assert_eq!(m.mean(1), expected_1);
+    }
+
+    #[test]
+    fn std() {
+        let m = DenseMatrix::from_2d_array(&[
+            &[1., 2., 3., 1., 2.],
+            &[4., 5., 6., 3., 4.],
+            &[7., 8., 9., 5., 6.],
+        ]);
+        let expected_0 = vec![2.44, 2.44, 2.44, 1.63, 1.63];
+        let expected_1 = vec![0.74, 1.01, 1.41];
+
+        assert!(m.std(0).approximate_eq(&expected_0, 1e-2));
+        assert!(m.std(1).approximate_eq(&expected_1, 1e-2));
+    }
+
+    #[test]
+    fn scale() {
+        let mut m = DenseMatrix::from_2d_array(&[&[1., 2., 3.], &[4., 5., 6.]]);
+        let expected_0 = DenseMatrix::from_2d_array(&[&[-1., -1., -1.], &[1., 1., 1.]]);
+        let expected_1 = DenseMatrix::from_2d_array(&[&[-1.22, 0.0, 1.22], &[-1.22, 0.0, 1.22]]);
+
+        {
+            let mut m = m.clone();
+            m.scale_mut(&m.mean(0), &m.std(0), 0);
+            assert!(m.approximate_eq(&expected_0, std::f32::EPSILON));
+        }
+
+        m.scale_mut(&m.mean(1), &m.std(1), 1);
+        assert!(m.approximate_eq(&expected_1, 1e-2));
+    }
+}