001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 * image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause).
006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009
010package imagingbook.common.math;
011
012//import org.apache.commons.math3.linear.MatrixUtils;
013import org.apache.commons.math3.linear.EigenDecomposition;
014import org.apache.commons.math3.linear.MatrixUtils;
015import org.apache.commons.math3.linear.RealMatrix;
016import org.apache.commons.math3.stat.correlation.Covariance;
017
018/**
019 * This class defines static methods for statistical calculations.
020 * 
021 * @author WB
022 *
023 */
024public abstract class Statistics {
025        
026        private Statistics() {}
027
028        /**
029         * Calculates the mean vector for a sequence of sample vectors.
030         * @param samples a 2D array of m-dimensional vectors ({@code }double[n][m]})
031         * @return the mean vector for the sample data (of length m)
032         */
033        public static double[] meanVector(double[][] samples) {
034                final int n = samples.length;
035                final int m = samples[0].length;
036                double[] mean = new double[m];
037                for (int k = 0; k < n; k++) {
038                        for (int i = 0; i < m; i++) {
039                                mean[i] = mean[i] + samples[k][i];
040                        }
041                }
042                for (int i = 0; i < m; i++) {
043                        mean[i] = mean[i] / n;
044                }
045                return mean;
046        }
047
048        /**
049         * Calculates the covariance matrix for a sequence of sample vectors. Takes a sequence of n data samples, each of
050         * dimension m. The data element {@code samples[i][j]} refers to the j-th component of sample i. No statistical bias
051         * correction is applied. Uses {@link Covariance} from Apache Commons Math.
052         *
053         * @param samples a 2D array of m-dimensional vectors ({@code }double[n][m]})
054         * @return the covariance matrix for the sample data (of dimension m x m)
055         * @see Covariance
056         */
057        public static double[][] covarianceMatrix(double[][] samples) {
058                return covarianceMatrix(samples, false);
059        }
060
061        /**
062         * Calculates the covariance matrix for a sequence of sample vectors. Takes a sequence of n data samples, each of
063         * dimension m. The data element {@code samples[i][j]} refers to the j-th component of sample i. Statistical bias
064         * correction is optionally applied. Uses {@link Covariance} from Apache Commons Math.
065         *
066         * @param samples a 2D array of m-dimensional vectors (double[n][m]).
067         * @param biasCorrect if {@code true}, statistical bias correction is applied.
068         * @return the covariance matrix for the sample data (of dimension m x m).
069         * @see Covariance
070         */
071        public static double[][] covarianceMatrix(double[][] samples, boolean biasCorrect) {
072                Covariance cov = new Covariance(samples, biasCorrect);
073                return cov.getCovarianceMatrix().getData();
074        }
075
076        /**
077         * Conditions the supplied covariance matrix by enforcing positive eigenvalues.
078         *
079         * @param cov original covariance matrix
080         * @param minDiagVal the minimum positive value of diagonal elements
081         * @return conditioned covariance matrix
082         */
083        public static double[][] conditionCovarianceMatrix(double[][] cov, double minDiagVal) {
084                RealMatrix C = MatrixUtils.createRealMatrix(cov);
085                return conditionCovarianceMatrix(C, minDiagVal).getData();
086        }
087
088        /**
089         * Conditions the supplied covariance matrix by enforcing positive eigenvalues.
090         *
091         * @param cov original covariance matrix
092         * @param minDiagVal the minimum positive value of diagonal elements
093         * @return conditioned covariance matrix
094         */
095        public static RealMatrix conditionCovarianceMatrix(RealMatrix cov, double minDiagVal) {
096                if (!cov.isSquare()) {
097                        throw new IllegalArgumentException("covariance matrix must be square");
098                }
099                EigenDecomposition ed = new EigenDecomposition(cov);  // S  ->  V . D . V^T
100                RealMatrix V  = ed.getV();
101                RealMatrix D  = ed.getD();      // diagonal matrix of eigenvalues
102                RealMatrix VT = ed.getVT();
103                for (int i = 0; i < D.getRowDimension(); i++) {
104                        D.setEntry(i, i, Math.max(D.getEntry(i, i), minDiagVal));       // setting eigenvalues to zero is not enough!
105                }
106                return V.multiply(D).multiply(VT);
107        }
108        
109//      /** 
110//       * example from UTICS-C Appendix:
111//       * N = 4 samples
112//       * K = 3 dimensions
113//       * @param args ignored
114//       */
115//      public static void main(String[] args) {
116//              
117//              boolean BIAS_CORRECT = false;
118//              
119//              // example: n = 4 samples of dimension m = 3:
120//              // samples[i][j], i = column (sample index), j = row (dimension index).
121//              double[][] samples = { 
122//                              {75, 37, 12},   // i = 0
123//                              {41, 27, 20},   // i = 1
124//                              {93, 81, 11},   // i = 2
125//                              {12, 48, 52}    // i = 3
126//              };
127//              
128//              // covariance matrix Cov (3x3)
129//              double[][] cov = covarianceMatrix(samples, BIAS_CORRECT);
130//              System.out.println("cov = \n" + Matrix.toString(cov));
131//              
132//              System.out.println();
133//              
134//              double[][] icov = Matrix.inverse(cov);
135//              System.out.println("icov = \n" + Matrix.toString(icov));
136//              
137//              System.out.println();
138//              
139//              double trace = MatrixUtils.createRealMatrix(cov).getTrace();
140//              System.out.println("trace(cov) = " + trace);
141//              double Fnorm = MatrixUtils.createRealMatrix(cov).getFrobeniusNorm();
142//              System.out.println("Fnorm(cov) = " + Fnorm);
143//      }
144        
145/* Results (bias-corrected):
146cov = 
147{{1296.250, 442.583, -627.250}, 
148{442.583, 550.250, -70.917}, 
149{-627.250, -70.917, 370.917}}
150
151icov = 
152{{0.024, -0.014, 0.038}, 
153{-0.014, 0.011, -0.022}, 
154{0.038, -0.022, 0.063}}
155*/
156        
157/* verified with Mathematica
158X1 = {75, 37, 12}; X2 = {41, 27, 20}; X3 = {93, 81, 11}; X4 = {12, 48, 52};
159samples = {X1, X2, X3, X4}
160N[Covariance[samples]]
161-> {{1296.25, 442.583, -627.25}, {442.583, 550.25, -70.9167}, {-627.25, -70.9167, 370.917}}
162*/
163        
164/* Results (NON bias-corrected):
165cov = 
166{{972.188, 331.938, -470.438}, 
167{331.938, 412.688, -53.188}, 
168{-470.438, -53.188, 278.188}}
169
170icov = {{0.032, -0.019, 0.051}, 
171{-0.019, 0.014, -0.030}, 
172{0.051, -0.030, 0.083}}
173*/
174
175}