org.apache.commons.math.linear.RealMatrix Java Examples

@Test
public void testLeastSquares2()
    throws FunctionEvaluationException {

    final RealMatrix factors =
        new Array2DRowRealMatrix(new double[][] {
                { 1.0, 0.0 },
                { 0.0, 1.0 }
            }, false);
    LeastSquaresConverter ls = new LeastSquaresConverter(new MultivariateVectorialFunction() {
            public double[] value(double[] variables) {
                return factors.operate(variables);
            }
        }, new double[] { 2.0, -3.0 }, new double[] { 10.0, 0.1 });
    NelderMead optimizer = new NelderMead();
    optimizer.setConvergenceChecker(new SimpleScalarValueChecker(-1.0, 1.0e-6));
    optimizer.setMaxEvaluations(200);
    RealPointValuePair optimum =
        optimizer.optimize(ls, GoalType.MINIMIZE, new double[] { 10.0, 10.0 });
    assertEquals( 2.0, optimum.getPointRef()[0], 5.0e-5);
    assertEquals(-3.0, optimum.getPointRef()[1], 8.0e-4);
    assertTrue(optimizer.getEvaluations() > 60);
    assertTrue(optimizer.getEvaluations() < 80);
    assertTrue(optimum.getValue() < 1.0e-6);
}
 

/**
 * Verifies that setting X, Y and covariance separately has the same effect as newSample(X,Y,cov).
 */
@Test
public void testNewSample2() throws Exception {
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    double[][] covariance = MatrixUtils.createRealIdentityMatrix(4).scalarMultiply(2).getData();
    GLSMultipleLinearRegression regression = new GLSMultipleLinearRegression();
    regression.newSampleData(y, x, covariance);
    RealMatrix combinedX = regression.X.copy();
    RealVector combinedY = regression.Y.copy();
    RealMatrix combinedCovInv = regression.getOmegaInverse();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    assertEquals(combinedX, regression.X);
    assertEquals(combinedY, regression.Y);
    assertEquals(combinedCovInv, regression.getOmegaInverse());
}
 

public void testOverdetermined() {
    final Random r    = new Random(5559252868205245l);
    int          p    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
    int          q    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
    RealMatrix   a    = createTestMatrix(r, p, q);
    RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);

    // build a perturbed system: A.X + noise = B
    RealMatrix b = a.multiply(xRef);
    final double noise = 0.001;
    b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(int row, int column, double value) {
            return value * (1.0 + noise * (2 * r.nextDouble() - 1));
        }
    });

    // despite perturbation, the least square solution should be pretty good
    RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
    assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);

}
 

public void testOverdetermined() {
    final Random r    = new Random(5559252868205245l);
    int          p    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
    int          q    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
    RealMatrix   a    = createTestMatrix(r, p, q);
    RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);

    // build a perturbed system: A.X + noise = B
    RealMatrix b = a.multiply(xRef);
    final double noise = 0.001;
    b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(int row, int column, double value) {
            return value * (1.0 + noise * (2 * r.nextDouble() - 1));
        }
    });

    // despite perturbation, the least square solution should be pretty good
    RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
    assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);

}
 

/**
 * Verify that direct t-tests using standard error estimates are consistent
 * with reported p-values
 */
public void testStdErrorConsistency() throws Exception {
    TDistribution tDistribution = new TDistributionImpl(45);
    RealMatrix matrix = createRealMatrix(swissData, 47, 5);
    PearsonsCorrelation corrInstance = new PearsonsCorrelation(matrix);
    RealMatrix rValues = corrInstance.getCorrelationMatrix();
    RealMatrix pValues = corrInstance.getCorrelationPValues();
    RealMatrix stdErrors = corrInstance.getCorrelationStandardErrors();
    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < i; j++) {
            double t = Math.abs(rValues.getEntry(i, j)) / stdErrors.getEntry(i, j);
            double p = 2 * (1 - tDistribution.cumulativeProbability(t));
            assertEquals(p, pValues.getEntry(i, j), 10E-15);
        }
    }
}
 

/**
 * <p>Compute the "hat" matrix.
 * </p>
 * <p>The hat matrix is defined in terms of the design matrix X
 *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
 * </p>
 * <p>The implementation here uses the QR decomposition to compute the
 * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
 * p-dimensional identity matrix augmented by 0's.  This computational
 * formula is from "The Hat Matrix in Regression and ANOVA",
 * David C. Hoaglin and Roy E. Welsch,
 * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
 *
 * @return the hat matrix
 */
public RealMatrix calculateHat() {
    // Create augmented identity matrix
    RealMatrix Q = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = Q.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j =0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }

    // Compute and return Hat matrix
    return Q.multiply(augI).multiply(Q.transpose());
}
 

/**
 * Verifies operation on indefinite matrix
 */
public void testZeroDivide() {
    RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] {
            { 0.0, 1.0, -1.0 }, 
            { 1.0, 1.0, 0.0 }, 
            { -1.0,0.0, 1.0 }        
    });
    EigenDecomposition ed = new EigenDecompositionImpl(indefinite, MathUtils.SAFE_MIN);
    checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12);
    double isqrt3 = 1/Math.sqrt(3.0);
    checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12);
    double isqrt2 = 1/Math.sqrt(2.0);
    checkEigenVector((new double[] {0.0,-isqrt2,-isqrt2}), ed, 1E-12);
    double isqrt6 = 1/Math.sqrt(6.0);
    checkEigenVector((new double[] {2*isqrt6,-isqrt6,isqrt6}), ed, 1E-12);
}
 

public void testMeanAndCorrelation() throws DimensionMismatchException {

        VectorialMean meanStat = new VectorialMean(mean.length);
        VectorialCovariance covStat = new VectorialCovariance(mean.length, true);
        for (int i = 0; i < 10000; ++i) {
            double[] v = generator.nextVector();
            meanStat.increment(v);
            covStat.increment(v);
        }

        double[] estimatedMean = meanStat.getResult();
        double scale;
        RealMatrix estimatedCorrelation = covStat.getResult();
        for (int i = 0; i < estimatedMean.length; ++i) {
            assertEquals(mean[i], estimatedMean[i], 0.07);
            for (int j = 0; j < i; ++j) {
                scale = standardDeviation[i] * standardDeviation[j];
                assertEquals(0, estimatedCorrelation.getEntry(i, j) / scale, 0.03);
            }
            scale = standardDeviation[i] * standardDeviation[i];
            assertEquals(1, estimatedCorrelation.getEntry(i, i) / scale, 0.025);
        }

    }
 

/** test that R is upper triangular */
public void testRUpperTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData3x4);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData4x3);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    Random r = new Random(643895747384642l);
    int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
    int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
    matrix = createTestMatrix(r, p, q);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = createTestMatrix(r, p, q);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

}
 

public void testMath226()
    throws DimensionMismatchException, NotPositiveDefiniteMatrixException {
    double[] mean = { 1, 1, 10, 1 };
    double[][] cov = {
            { 1, 3, 2, 6 },
            { 3, 13, 16, 2 },
            { 2, 16, 38, -1 },
            { 6, 2, -1, 197 }
    };
    RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
    JDKRandomGenerator jg = new JDKRandomGenerator();
    jg.setSeed(5322145245211l);
    NormalizedRandomGenerator rg = new GaussianRandomGenerator(jg);
    CorrelatedRandomVectorGenerator sg =
        new CorrelatedRandomVectorGenerator(mean, covRM, 0.00001, rg);

    for (int i = 0; i < 10; i++) {
        double[] generated = sg.nextVector();
        assertTrue(Math.abs(generated[0] - 1) > 0.1);
    }

}
 

/**
 * Returns a matrix of p-values associated with the (two-sided) null
 * hypothesis that the corresponding correlation coefficient is zero.
 * <p><code>getCorrelationPValues().getEntry(i,j)</code> is the probability
 * that a random variable distributed as <code>t<sub>n-2</sub></code> takes
 * a value with absolute value greater than or equal to <br>
 * <code>|r|((n - 2) / (1 - r<sup>2</sup>))<sup>1/2</sup></code></p>
 * <p>The values in the matrix are sometimes referred to as the
 * <i>significance</i> of the corresponding correlation coefficients.</p>
 *
 * @return matrix of p-values
 * @throws MathException if an error occurs estimating probabilities
 */
public RealMatrix getCorrelationPValues() throws MathException {
    TDistribution tDistribution = new TDistributionImpl(nObs - 2);
    int nVars = correlationMatrix.getColumnDimension();
    double[][] out = new double[nVars][nVars];
    for (int i = 0; i < nVars; i++) {
        for (int j = 0; j < nVars; j++) {
            if (i == j) {
                out[i][j] = 0d;
            } else {
                double r = correlationMatrix.getEntry(i, j);
                double t = Math.abs(r * Math.sqrt((nObs - 2)/(1 - r * r)));
                out[i][j] = 2 * (1 - tDistribution.cumulativeProbability(t));
            }
        }
    }
    return new BlockRealMatrix(out);
}
 

/**
 * Get the covariance matrix of the optimized parameters.
 *
 * @return the covariance matrix.
 * @throws org.apache.commons.math.linear.SingularMatrixException
 * if the covariance matrix cannot be computed (singular problem).
 * @throws org.apache.commons.math.exception.MathUserException if the
 * jacobian function throws one.
 */
public double[][] getCovariances() {
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    // compute the covariances matrix
    RealMatrix inverse =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
    return inverse.getData();
}
 

private void checkBiDiagonal(RealMatrix m) {
    final int rows = m.getRowDimension();
    final int cols = m.getColumnDimension();
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            if (rows < cols) {
                if ((i < j) || (i > j + 1)) {
                    Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
                }
            } else {
                if ((i < j - 1) || (i > j)) {
                    Assert.assertEquals(0, m.getEntry(i, j), 1.0e-16);
                }
            }
        }
    }
}
 

private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
    RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
    m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
        @Override
        public double visit(int row, int column, double value)
            throws MatrixVisitorException {
            return 2.0 * r.nextDouble() - 1.0;
        }
    });
    return m;
}
 

private void checkdimensions(RealMatrix matrix) {
    final int m = matrix.getRowDimension();
    final int n = matrix.getColumnDimension();
    BiDiagonalTransformer transformer = new BiDiagonalTransformer(matrix);
    assertEquals(m, transformer.getU().getRowDimension());
    assertEquals(m, transformer.getU().getColumnDimension());
    assertEquals(m, transformer.getB().getRowDimension());
    assertEquals(n, transformer.getB().getColumnDimension());
    assertEquals(n, transformer.getV().getRowDimension());
    assertEquals(n, transformer.getV().getColumnDimension());

}
 

@Override
public void testConsistency() {
    RealMatrix matrix = createRealMatrix(longleyData, 16, 7);
    SpearmansCorrelation corrInstance = new SpearmansCorrelation(matrix);
    double[][] data = matrix.getData();
    double[] x = matrix.getColumn(0);
    double[] y = matrix.getColumn(1);
    assertEquals(new SpearmansCorrelation().correlation(x, y),
            corrInstance.getCorrelationMatrix().getEntry(0, 1), Double.MIN_VALUE);
    TestUtils.assertEquals("Correlation matrix", corrInstance.getCorrelationMatrix(),
            new SpearmansCorrelation().computeCorrelationMatrix(data), Double.MIN_VALUE);
}
 

private void checkdimensions(RealMatrix matrix) {
    final int m = matrix.getRowDimension();
    final int n = matrix.getColumnDimension();
    BiDiagonalTransformer transformer = new BiDiagonalTransformer(matrix);
    Assert.assertEquals(m, transformer.getU().getRowDimension());
    Assert.assertEquals(m, transformer.getU().getColumnDimension());
    Assert.assertEquals(m, transformer.getB().getRowDimension());
    Assert.assertEquals(n, transformer.getB().getColumnDimension());
    Assert.assertEquals(n, transformer.getV().getRowDimension());
    Assert.assertEquals(n, transformer.getV().getColumnDimension());

}
 

/**
 * Returns an array whose i<sup>th</sup> entry is the standard deviation of the
 * i<sup>th</sup> entries of the arrays that have been added using
 * {@link #addValue(double[])}
 *
 * @return the array of component standard deviations
 */
public double[] getStandardDeviation() {
    double[] stdDev = new double[k];
    if (getN() < 1) {
        Arrays.fill(stdDev, Double.NaN);
    } else if (getN() < 2) {
        Arrays.fill(stdDev, 0.0);
    } else {
        RealMatrix matrix = covarianceImpl.getResult();
        for (int i = 0; i < k; ++i) {
            stdDev[i] = Math.sqrt(matrix.getEntry(i, i));
        }
    }
    return stdDev;
}
 

/**
 * Computes the correlation matrix for the columns of the
 * input matrix.
 * 
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    int nVars = matrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        for (int j = 0; j < i; j++) {
          double corr = correlation(matrix.getColumn(i), matrix.getColumn(j));
          outMatrix.setEntry(i, j, corr);
          outMatrix.setEntry(j, i, corr);
        }
        outMatrix.setEntry(i, i, 1d);
    }
    return outMatrix;
}
 

/**
 * Throws IllegalArgumentException of the matrix does not have at least
 * two columns and two rows
 *
 * @param matrix matrix to check for sufficiency
 */
private void checkSufficientData(final RealMatrix matrix) {
    int nRows = matrix.getRowDimension();
    int nCols = matrix.getColumnDimension();
    if (nRows < 2 || nCols < 2) {
        throw MathRuntimeException.createIllegalArgumentException(
                "insufficient data: only {0} rows and {1} columns.",
                nRows, nCols);
    }
}
 

/**
 * Replace NaNs in given matrix with given value
 */
public static void setNaNs(RealMatrix matrix, int val) {
    for (int i = 0; i < matrix.getRowDimension(); i++) {
        for (int j = 0; j < matrix.getColumnDimension(); j++) {
            if (Double.isNaN(matrix.getEntry(i, j))) {
                matrix.setEntry(i, j, val);
            }
        }
    }
}
 

public void testDimension3() {
    RealMatrix matrix =
        MatrixUtils.createRealMatrix(new double[][] {
                               {  39632.0, -4824.0, -16560.0 },
                               {  -4824.0,  8693.0,   7920.0 },
                               { -16560.0,  7920.0,  17300.0 }
                           });
    EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
    assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
    assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
    assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

/**
 * @param mat
 *            Input matrix.
 * @param n
 *            Number of row replicates.
 * @param m
 *            Number of column replicates.
 * @return Matrix which replicates the input matrix in both directions.
 */
private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
    int rd = mat.getRowDimension();
    int cd = mat.getColumnDimension();
    double[][] d = new double[n * rd][m * cd];
    for (int r = 0; r < n * rd; r++) {
        for (int c = 0; c < m * cd; c++) {
            d[r][c] = mat.getEntry(r % rd, c % cd);
        }
    }
    return new Array2DRowRealMatrix(d, false);
}
 

/**
 * Computes the correlation matrix for the columns of the
 * input matrix.
 *
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    int nVars = matrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        for (int j = 0; j < i; j++) {
          double corr = correlation(matrix.getColumn(i), matrix.getColumn(j));
          outMatrix.setEntry(i, j, corr);
          outMatrix.setEntry(j, i, corr);
        }
        outMatrix.setEntry(i, i, 1d);
    }
    return outMatrix;
}
 

public void checkAEqualUSVt(final RealMatrix matrix) {
    SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
    RealMatrix u = svd.getU();
    RealMatrix s = svd.getS();
    RealMatrix v = svd.getV();
    double norm = u.multiply(s).multiply(v.transpose()).subtract(matrix).getNorm();
    assertEquals(0, norm, normTolerance);

}
 

@Test
public void testSimplistic() {
    VectorialCovariance stat = new VectorialCovariance(2, true);
    stat.increment(new double[] {-1.0,  1.0});
    stat.increment(new double[] { 1.0, -1.0});
    RealMatrix c = stat.getResult();
    Assert.assertEquals( 2.0, c.getEntry(0, 0), 1.0e-12);
    Assert.assertEquals(-2.0, c.getEntry(1, 0), 1.0e-12);
    Assert.assertEquals( 2.0, c.getEntry(1, 1), 1.0e-12);
}
 

/**
 * Matrix with eigenvalues {8, -1, -1}
 */
public void testRepeatedEigenvalue() {
    RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
            {3,  2,  4},
            {2,  0,  2},
            {4,  2,  3}
    });
    EigenDecomposition ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
    checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
    checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
}
 

/**
 * Throws IllegalArgumentException of the matrix does not have at least
 * two columns and two rows
 *
 * @param matrix matrix to check for sufficiency
 */
private void checkSufficientData(final RealMatrix matrix) {
    int nRows = matrix.getRowDimension();
    int nCols = matrix.getColumnDimension();
    if (nRows < 2 || nCols < 2) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.INSUFFICIENT_ROWS_AND_COLUMNS,
                nRows, nCols);
    }
}