PearsonsCorrelation.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  17. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */
  21. package org.hipparchus.stat.correlation;

  22. import org.hipparchus.distribution.continuous.TDistribution;
  23. import org.hipparchus.exception.LocalizedCoreFormats;
  24. import org.hipparchus.exception.MathIllegalArgumentException;
  25. import org.hipparchus.linear.BlockRealMatrix;
  26. import org.hipparchus.linear.RealMatrix;
  27. import org.hipparchus.stat.LocalizedStatFormats;
  28. import org.hipparchus.stat.regression.SimpleRegression;
  29. import org.hipparchus.util.FastMath;
  30. import org.hipparchus.util.MathArrays;
  31. import org.hipparchus.util.MathUtils;

  32. /**
  33.  * Computes Pearson's product-moment correlation coefficients for pairs of arrays
  34.  * or columns of a matrix.
  35.  * <p>
  36.  * The constructors that take <code>RealMatrix</code> or
  37.  * <code>double[][]</code> arguments generate correlation matrices.  The
  38.  * columns of the input matrices are assumed to represent variable values.
  39.  * Correlations are given by the formula:
  40.  * <p>
  41.  * <code>cor(X, Y) = &Sigma;[(x<sub>i</sub> - E(X))(y<sub>i</sub> - E(Y))] / [(n - 1)s(X)s(Y)]</code>
  42.  * <p>
  43.  * where <code>E(X)</code> is the mean of <code>X</code>, <code>E(Y)</code>
  44.  * is the mean of the <code>Y</code> values and s(X), s(Y) are standard deviations.
  45.  * <p>
  46.  * To compute the correlation coefficient for a single pair of arrays, use {@link #PearsonsCorrelation()}
  47.  * to construct an instance with no data and then {@link #correlation(double[], double[])}.
  48.  * Correlation matrices can also be computed directly from an instance with no data using
  49.  * {@link #computeCorrelationMatrix(double[][])}. In order to use {@link #getCorrelationMatrix()},
  50.  * {@link #getCorrelationPValues()},  or {@link #getCorrelationStandardErrors()}; however, one of the
  51.  * constructors supplying data or a covariance matrix must be used to create the instance.
  52.  */
  53. public class PearsonsCorrelation {

  54.     /** correlation matrix */
  55.     private final RealMatrix correlationMatrix;

  56.     /** number of observations */
  57.     private final int nObs;

  58.     /**
  59.      * Create a PearsonsCorrelation instance without data.
  60.      */
  61.     public PearsonsCorrelation() {
  62.         super();
  63.         correlationMatrix = null;
  64.         nObs = 0;
  65.     }

  66.     /**
  67.      * Create a PearsonsCorrelation from a rectangular array
  68.      * whose columns represent values of variables to be correlated.
  69.      *
  70.      * Throws MathIllegalArgumentException if the input array does not have at least
  71.      * two columns and two rows.  Pairwise correlations are set to NaN if one
  72.      * of the correlates has zero variance.
  73.      *
  74.      * @param data rectangular array with columns representing variables
  75.      * @throws MathIllegalArgumentException if the input data array is not
  76.      * rectangular with at least two rows and two columns.
  77.      * @see #correlation(double[], double[])
  78.      */
  79.     public PearsonsCorrelation(double[][] data) {
  80.         this(new BlockRealMatrix(data));
  81.     }

  82.     /**
  83.      * Create a PearsonsCorrelation from a RealMatrix whose columns
  84.      * represent variables to be correlated.
  85.      *
  86.      * Throws MathIllegalArgumentException if the matrix does not have at least
  87.      * two columns and two rows.  Pairwise correlations are set to NaN if one
  88.      * of the correlates has zero variance.
  89.      *
  90.      * @param matrix matrix with columns representing variables to correlate
  91.      * @throws MathIllegalArgumentException if the matrix does not contain sufficient data
  92.      * @see #correlation(double[], double[])
  93.      */
  94.     public PearsonsCorrelation(RealMatrix matrix) {
  95.         nObs = matrix.getRowDimension();
  96.         correlationMatrix = computeCorrelationMatrix(matrix);
  97.     }

  98.     /**
  99.      * Create a PearsonsCorrelation from a {@link Covariance}.  The correlation
  100.      * matrix is computed by scaling the Covariance's covariance matrix.
  101.      * The Covariance instance must have been created from a data matrix with
  102.      * columns representing variable values.
  103.      *
  104.      * @param covariance Covariance instance
  105.      */
  106.     public PearsonsCorrelation(Covariance covariance) {
  107.         RealMatrix covarianceMatrix = covariance.getCovarianceMatrix();
  108.         MathUtils.checkNotNull(covarianceMatrix, LocalizedStatFormats.COVARIANCE_MATRIX);
  109.         nObs = covariance.getN();
  110.         correlationMatrix = covarianceToCorrelation(covarianceMatrix);
  111.     }

  112.     /**
  113.      * Create a PearsonsCorrelation from a covariance matrix. The correlation
  114.      * matrix is computed by scaling the covariance matrix.
  115.      *
  116.      * @param covarianceMatrix covariance matrix
  117.      * @param numberOfObservations the number of observations in the dataset used to compute
  118.      * the covariance matrix
  119.      */
  120.     public PearsonsCorrelation(RealMatrix covarianceMatrix, int numberOfObservations) {
  121.         nObs = numberOfObservations;
  122.         correlationMatrix = covarianceToCorrelation(covarianceMatrix);
  123.     }

  124.     /**
  125.      * Returns the correlation matrix.
  126.      *
  127.      * <p>This method will return null if the argumentless constructor was used
  128.      * to create this instance, even if {@link #computeCorrelationMatrix(double[][])}
  129.      * has been called before it is activated.</p>
  130.      *
  131.      * @return correlation matrix
  132.      */
  133.     public RealMatrix getCorrelationMatrix() {
  134.         return correlationMatrix;
  135.     }

  136.     /**
  137.      * Returns a matrix of standard errors associated with the estimates
  138.      * in the correlation matrix.<br>
  139.      * <code>getCorrelationStandardErrors().getEntry(i,j)</code> is the standard
  140.      * error associated with <code>getCorrelationMatrix.getEntry(i,j)</code>
  141.      *
  142.      * <p>The formula used to compute the standard error is <br>
  143.      * <code>SE<sub>r</sub> = ((1 - r<sup>2</sup>) / (n - 2))<sup>1/2</sup></code>
  144.      * where <code>r</code> is the estimated correlation coefficient and
  145.      * <code>n</code> is the number of observations in the source dataset.</p>
  146.      *
  147.      * <p>To use this method, one of the constructors that supply an input
  148.      * matrix must have been used to create this instance.</p>
  149.      *
  150.      * @return matrix of correlation standard errors
  151.      * @throws NullPointerException if this instance was created with no data
  152.      */
  153.     public RealMatrix getCorrelationStandardErrors() {
  154.         int nVars = correlationMatrix.getColumnDimension();
  155.         double[][] out = new double[nVars][nVars];
  156.         for (int i = 0; i < nVars; i++) {
  157.             for (int j = 0; j < nVars; j++) {
  158.                 double r = correlationMatrix.getEntry(i, j);
  159.                 out[i][j] = FastMath.sqrt((1 - r * r) /(nObs - 2));
  160.             }
  161.         }
  162.         return new BlockRealMatrix(out);
  163.     }

  164.     /**
  165.      * Returns a matrix of p-values associated with the (two-sided) null
  166.      * hypothesis that the corresponding correlation coefficient is zero.
  167.      *
  168.      * <p><code>getCorrelationPValues().getEntry(i,j)</code> is the probability
  169.      * that a random variable distributed as <code>t<sub>n-2</sub></code> takes
  170.      * a value with absolute value greater than or equal to <br>
  171.      * <code>|r|((n - 2) / (1 - r<sup>2</sup>))<sup>1/2</sup></code></p>
  172.      *
  173.      * <p>The values in the matrix are sometimes referred to as the
  174.      * <i>significance</i> of the corresponding correlation coefficients.</p>
  175.      *
  176.      * <p>To use this method, one of the constructors that supply an input
  177.      * matrix must have been used to create this instance.</p>
  178.      *
  179.      * @return matrix of p-values
  180.      * @throws org.hipparchus.exception.MathIllegalStateException
  181.      * if an error occurs estimating probabilities
  182.      * @throws NullPointerException if this instance was created with no data
  183.      */
  184.     public RealMatrix getCorrelationPValues() {
  185.         TDistribution tDistribution = new TDistribution(nObs - 2);
  186.         int nVars = correlationMatrix.getColumnDimension();
  187.         double[][] out = new double[nVars][nVars];
  188.         for (int i = 0; i < nVars; i++) {
  189.             for (int j = 0; j < nVars; j++) {
  190.                 if (i == j) {
  191.                     out[i][j] = 0d;
  192.                 } else {
  193.                     double r = correlationMatrix.getEntry(i, j);
  194.                     double t = FastMath.abs(r * FastMath.sqrt((nObs - 2)/(1 - r * r)));
  195.                     out[i][j] = 2 * tDistribution.cumulativeProbability(-t);
  196.                 }
  197.             }
  198.         }
  199.         return new BlockRealMatrix(out);
  200.     }


  201.     /**
  202.      * Computes the correlation matrix for the columns of the
  203.      * input matrix, using {@link #correlation(double[], double[])}.
  204.      *
  205.      * Throws MathIllegalArgumentException if the matrix does not have at least
  206.      * two columns and two rows.  Pairwise correlations are set to NaN if one
  207.      * of the correlates has zero variance.
  208.      *
  209.      * @param matrix matrix with columns representing variables to correlate
  210.      * @return correlation matrix
  211.      * @throws MathIllegalArgumentException if the matrix does not contain sufficient data
  212.      * @see #correlation(double[], double[])
  213.      */
  214.     public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
  215.         checkSufficientData(matrix);
  216.         int nVars = matrix.getColumnDimension();
  217.         RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
  218.         for (int i = 0; i < nVars; i++) {
  219.             for (int j = 0; j < i; j++) {
  220.               double corr = correlation(matrix.getColumn(i), matrix.getColumn(j));
  221.               outMatrix.setEntry(i, j, corr);
  222.               outMatrix.setEntry(j, i, corr);
  223.             }
  224.             outMatrix.setEntry(i, i, 1d);
  225.         }
  226.         return outMatrix;
  227.     }

  228.     /**
  229.      * Computes the correlation matrix for the columns of the
  230.      * input rectangular array.  The columns of the array represent values
  231.      * of variables to be correlated.
  232.      *
  233.      * Throws MathIllegalArgumentException if the matrix does not have at least
  234.      * two columns and two rows or if the array is not rectangular. Pairwise
  235.      * correlations are set to NaN if one of the correlates has zero variance.
  236.      *
  237.      * @param data matrix with columns representing variables to correlate
  238.      * @return correlation matrix
  239.      * @throws MathIllegalArgumentException if the array does not contain sufficient data
  240.      * @see #correlation(double[], double[])
  241.      */
  242.     public RealMatrix computeCorrelationMatrix(double[][] data) {
  243.        return computeCorrelationMatrix(new BlockRealMatrix(data));
  244.     }

  245.     /**
  246.      * Computes the Pearson's product-moment correlation coefficient between two arrays.
  247.      *
  248.      * <p>Throws MathIllegalArgumentException if the arrays do not have the same length
  249.      * or their common length is less than 2.  Returns {@code NaN} if either of the arrays
  250.      * has zero variance (i.e., if one of the arrays does not contain at least two distinct
  251.      * values).</p>
  252.      *
  253.      * @param xArray first data array
  254.      * @param yArray second data array
  255.      * @return Returns Pearson's correlation coefficient for the two arrays
  256.      * @throws MathIllegalArgumentException if the arrays lengths do not match
  257.      * @throws MathIllegalArgumentException if there is insufficient data
  258.      */
  259.     public double correlation(final double[] xArray, final double[] yArray) {
  260.         MathArrays.checkEqualLength(xArray, yArray);
  261.         if (xArray.length < 2) {
  262.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DIMENSION,
  263.                                                    xArray.length, 2);
  264.         }

  265.         SimpleRegression regression = new SimpleRegression();
  266.         for(int i = 0; i < xArray.length; i++) {
  267.             regression.addData(xArray[i], yArray[i]);
  268.         }
  269.         return regression.getR();
  270.     }

  271.     /**
  272.      * Derives a correlation matrix from a covariance matrix.
  273.      *
  274.      * <p>Uses the formula <br>
  275.      * <code>r(X,Y) = cov(X,Y)/s(X)s(Y)</code> where
  276.      * <code>r(&middot;,&middot;)</code> is the correlation coefficient and
  277.      * <code>s(&middot;)</code> means standard deviation.</p>
  278.      *
  279.      * @param covarianceMatrix the covariance matrix
  280.      * @return correlation matrix
  281.      */
  282.     public RealMatrix covarianceToCorrelation(RealMatrix covarianceMatrix) {
  283.         int nVars = covarianceMatrix.getColumnDimension();
  284.         RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
  285.         for (int i = 0; i < nVars; i++) {
  286.             double sigma = FastMath.sqrt(covarianceMatrix.getEntry(i, i));
  287.             outMatrix.setEntry(i, i, 1d);
  288.             for (int j = 0; j < i; j++) {
  289.                 double entry = covarianceMatrix.getEntry(i, j) /
  290.                        (sigma * FastMath.sqrt(covarianceMatrix.getEntry(j, j)));
  291.                 outMatrix.setEntry(i, j, entry);
  292.                 outMatrix.setEntry(j, i, entry);
  293.             }
  294.         }
  295.         return outMatrix;
  296.     }

  297.     /**
  298.      * Throws MathIllegalArgumentException if the matrix does not have at least
  299.      * two columns and two rows.
  300.      *
  301.      * @param matrix matrix to check for sufficiency
  302.      * @throws MathIllegalArgumentException if there is insufficient data
  303.      */
  304.     private void checkSufficientData(final RealMatrix matrix) {
  305.         int nRows = matrix.getRowDimension();
  306.         int nCols = matrix.getColumnDimension();
  307.         if (nRows < 2 || nCols < 2) {
  308.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_ROWS_AND_COLUMNS,
  309.                                                    nRows, nCols);
  310.         }
  311.     }
  312. }