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

The following examples show how to use org.apache.commons.math.linear.RealMatrix. These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
/**
 * 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());
}
 
Example 2
Source Project: astor   Source File: QRSolverTest.java    License: GNU General Public License v2.0 6 votes vote down vote up
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);

}
 
Example 3
Source Project: astor   Source File: NelderMeadTest.java    License: GNU General Public License v2.0 6 votes vote down vote up
@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);
}
 
Example 4
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);
                }
            }
        }
    }
}
 
Example 5
Source Project: astor   Source File: PearsonsCorrelation.java    License: GNU General Public License v2.0 6 votes vote down vote up
/**
 * 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);
}
 
Example 6
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);
    }

}
 
Example 7
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);
        }

    }
 
Example 8
/** 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());

}
 
Example 9
/**
 * 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);
}
 
Example 10
/**
 * <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());
}
 
Example 11
/**
 * 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);
        }
    }
}
 
Example 12
Source Project: astor   Source File: QRSolverTest.java    License: GNU General Public License v2.0 6 votes vote down vote up
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);

}
 
Example 13
/**
 * 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();
}
 
Example 14
Source Project: astor   Source File: PearsonsCorrelation.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * 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;
}
 
Example 15
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());

}
 
Example 16
/**
 * Test Longley dataset against R.
 */
public void testLongly() throws Exception {
    RealMatrix matrix = createRealMatrix(longleyData, 16, 7);
    PearsonsCorrelation corrInstance = new PearsonsCorrelation(matrix);
    RealMatrix correlationMatrix = corrInstance.getCorrelationMatrix();
    double[] rData = new double[] {
            1.000000000000000, 0.9708985250610560, 0.9835516111796693, 0.5024980838759942,
            0.4573073999764817, 0.960390571594376, 0.9713294591921188,
            0.970898525061056, 1.0000000000000000, 0.9915891780247822, 0.6206333925590966,
            0.4647441876006747, 0.979163432977498, 0.9911491900672053,
            0.983551611179669, 0.9915891780247822, 1.0000000000000000, 0.6042609398895580,
            0.4464367918926265, 0.991090069458478, 0.9952734837647849,
            0.502498083875994, 0.6206333925590966, 0.6042609398895580, 1.0000000000000000,
            -0.1774206295018783, 0.686551516365312, 0.6682566045621746,
            0.457307399976482, 0.4647441876006747, 0.4464367918926265, -0.1774206295018783,
            1.0000000000000000, 0.364416267189032, 0.4172451498349454,
            0.960390571594376, 0.9791634329774981, 0.9910900694584777, 0.6865515163653120,
            0.3644162671890320, 1.000000000000000, 0.9939528462329257,
            0.971329459192119, 0.9911491900672053, 0.9952734837647849, 0.6682566045621746,
            0.4172451498349454, 0.993952846232926, 1.0000000000000000
    };
    TestUtils.assertEquals("correlation matrix", createRealMatrix(rData, 7, 7), correlationMatrix, 10E-15);

    double[] rPvalues = new double[] {
            4.38904690369668e-10,
            8.36353208910623e-12, 7.8159700933611e-14,
            0.0472894097790304, 0.01030636128354301, 0.01316878049026582,
            0.0749178049642416, 0.06971758330341182, 0.0830166169296545, 0.510948586323452,
            3.693245043123738e-09, 4.327782576751815e-11, 1.167954621905665e-13, 0.00331028281967516, 0.1652293725106684,
            3.95834476307755e-10, 1.114663916723657e-13, 1.332267629550188e-15, 0.00466039138541463, 0.1078477071581498, 7.771561172376096e-15
    };
    RealMatrix rPMatrix = createLowerTriangularRealMatrix(rPvalues, 7);
    fillUpper(rPMatrix, 0d);
    TestUtils.assertEquals("correlation p values", rPMatrix, corrInstance.getCorrelationPValues(), 10E-15);
}
 
Example 17
/** Test based on a dimension 4 Hadamard matrix. */
public void testHadamard() {
    RealMatrix matrix = new Array2DRowRealMatrix(new double[][] {
            {15.0 / 2.0,  5.0 / 2.0,  9.0 / 2.0,  3.0 / 2.0 },
            { 5.0 / 2.0, 15.0 / 2.0,  3.0 / 2.0,  9.0 / 2.0 },
            { 9.0 / 2.0,  3.0 / 2.0, 15.0 / 2.0,  5.0 / 2.0 },
            { 3.0 / 2.0,  9.0 / 2.0,  5.0 / 2.0, 15.0 / 2.0 }
    }, false);
    SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
    assertEquals(16.0, svd.getSingularValues()[0], 1.0e-14);
    assertEquals( 8.0, svd.getSingularValues()[1], 1.0e-14);
    assertEquals( 4.0, svd.getSingularValues()[2], 1.0e-14);
    assertEquals( 2.0, svd.getSingularValues()[3], 1.0e-14);

    RealMatrix fullCovariance = new Array2DRowRealMatrix(new double[][] {
            {  85.0 / 1024, -51.0 / 1024, -75.0 / 1024,  45.0 / 1024 },
            { -51.0 / 1024,  85.0 / 1024,  45.0 / 1024, -75.0 / 1024 },
            { -75.0 / 1024,  45.0 / 1024,  85.0 / 1024, -51.0 / 1024 },
            {  45.0 / 1024, -75.0 / 1024, -51.0 / 1024,  85.0 / 1024 }
    }, false);
    assertEquals(0.0,
                 fullCovariance.subtract(svd.getCovariance(0.0)).getNorm(),
                 1.0e-14);

    RealMatrix halfCovariance = new Array2DRowRealMatrix(new double[][] {
            {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
            {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 },
            {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
            {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 }
    }, false);
    assertEquals(0.0,
                 halfCovariance.subtract(svd.getCovariance(6.0)).getNorm(),
                 1.0e-14);

}
 
Example 18
/** test that LT is transpose of L */
@Test
public void testLTTransposed() throws MathException {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    CholeskyDecomposition llt = new CholeskyDecompositionImpl(matrix);
    RealMatrix l  = llt.getL();
    RealMatrix lt = llt.getLT();
    double norm = l.subtract(lt.transpose()).getNorm();
    assertEquals(0, norm, 1.0e-15);
}
 
Example 19
Source Project: astor   Source File: AbstractEstimator.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int n = problem.getMeasurements().length;
    final int m = problem.getUnboundParameters().length;
    final int max  = m * n;
    double[][] jTj = new double[m][m];
    for (int i = 0; i < m; ++i) {
        for (int j = i; j < m; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += m) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem");
    }

}
 
Example 20
Source Project: astor   Source File: CMAESOptimizer.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * @param m
 *            Input matrix.
 * @return Row matrix representing the sums of the rows.
 */
private static RealMatrix sumRows(final RealMatrix m) {
    double[][] d = new double[1][m.getColumnDimension()];
    for (int c = 0; c < m.getColumnDimension(); c++) {
        double sum = 0;
        for (int r = 0; r < m.getRowDimension(); r++) {
            sum += m.getEntry(r, c);
        }
        d[0][c] = sum;
    }
    return new Array2DRowRealMatrix(d, false);
}
 
Example 21
Source Project: astor   Source File: PearsonsCorrelation.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * 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);
    }
}
 
Example 22
@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);
}
 
Example 23
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());

}
 
Example 24
/** test matrices values */
public void testMatricesValues1() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 0.0, 0.0 },
            { 0.5, 1.0, 0.0 },
            { 0.5, 0.2, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 2.0,  5.0, 3.0 },
            { 0.0, -2.5, 6.5 },
            { 0.0,  0.0, 0.2 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 1, 2, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 
Example 25
/**
 * Verifies that calculateYVariance and calculateResidualVariance return consistent
 * values with direct variance computation from Y, residuals, respectively.
 */
protected void checkVarianceConsistency(OLSMultipleLinearRegression model) throws Exception {
    // Check Y variance consistency
    TestUtils.assertEquals(StatUtils.variance(model.Y.getData()), model.calculateYVariance(), 0);
    
    // Check residual variance consistency
    double[] residuals = model.calculateResiduals().getData();
    RealMatrix X = model.X;
    TestUtils.assertEquals(
            StatUtils.variance(model.calculateResiduals().getData()) * (residuals.length - 1),
            model.calculateErrorVariance() * (X.getRowDimension() - X.getColumnDimension()), 1E-20);
    
}
 
Example 26
/** test that U is upper triangular */
public void testUUpperTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix u = new LUDecompositionImpl(matrix).getU();
    for (int i = 0; i < u.getRowDimension(); i++) {
        for (int j = 0; j < i; j++) {
            assertEquals(u.getEntry(i, j), 0, entryTolerance);
        }
    }
}
 
Example 27
/**
 * Verifies that newSampleData methods consistently insert unitary columns
 * in design matrix.  Confirms the fix for MATH-411.
 */
@Test
public void testNewSample() throws Exception {
    double[] design = new double[] {
      1, 19, 22, 33,
      2, 20, 30, 40,
      3, 25, 35, 45,
      4, 27, 37, 47
    };
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    AbstractMultipleLinearRegression regression = createRegression();
    regression.newSampleData(design, 4, 3);
    RealMatrix flatX = regression.X.copy();
    RealVector flatY = regression.Y.copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.X);
    Assert.assertEquals(flatY, regression.Y);
    
    // No intercept
    regression.setNoIntercept(true);
    regression.newSampleData(design, 4, 3);
    flatX = regression.X.copy();
    flatY = regression.Y.copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.X);
    Assert.assertEquals(flatY, regression.Y);
}
 
Example 28
Source Project: Juicebox   Source File: MatrixTools.java    License: MIT License 5 votes vote down vote up
/**
 * From Matrix M, extract out M[r1:r2,c1:c2]
 * r2, c2 not inclusive (~python numpy, not like Matlab)
 *
 * @return extracted matrix region M[r1:r2,c1:c2]
 */
public static RealMatrix extractLocalMatrixRegion(RealMatrix matrix, int r1, int r2, int c1, int c2) {

    int numRows = r2 - r1;
    int numColumns = c2 - c1;
    RealMatrix extractedRegion = cleanArray2DMatrix(numRows, numColumns);

    for (int i = 0; i < numRows; i++) {
        for (int j = 0; j < numColumns; j++) {
            extractedRegion.setEntry(i, j, matrix.getEntry(r1 + i, c1 + j));
        }
    }
    return extractedRegion;
}
 
Example 29
Source Project: astor   Source File: Covariance.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * Throws IllegalArgumentException of the matrix does not have at least
 * two columns and two rows
 * @param matrix matrix to check
 */
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);
    }
}
 
Example 30
Source Project: astor   Source File: AbstractEstimator.java    License: GNU General Public License v2.0 5 votes vote down vote up
/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {
 
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int rows = problem.getMeasurements().length;
    final int cols = problem.getUnboundParameters().length;
    final int max  = cols * rows;
    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 < max; k += cols) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem");
    }

}