List of usage examples for org.apache.commons.math3.linear RealMatrix subtract
RealMatrix subtract(RealMatrix m) throws MatrixDimensionMismatchException;
From source file:com.itemanalysis.psychometrics.mixture.MvNormalComponentDistribution.java
/** * * @param x a matrix of dimension 1 x k, where k is the number of variables * @return// w w w . ja va2 s. com */ public double density(RealMatrix x) throws SingularMatrixException { double prob = 0.0; RealMatrix xTran = x.transpose(); int d = xTran.getRowDimension(); double det = new LUDecomposition(sigma).getDeterminant(); double nconst = 1.0 / Math.sqrt(det * Math.pow(2.0 * Math.PI, d)); RealMatrix Sinv = new LUDecomposition(sigma).getSolver().getInverse(); RealMatrix delta = xTran.subtract(mu); RealMatrix dist = (delta.transpose().multiply(Sinv).multiply(delta)); prob = nconst * Math.exp(-0.5 * dist.getEntry(0, 0)); return prob; }
From source file:edu.cudenver.bios.power.glmm.GLMMTest.java
/** * Calculate the sum of squares hypothesis matrix (the H matrix) * @param params matrices input by user//w w w. ja v a2 s.c o m * @return H matrix */ public RealMatrix getHypothesisSumOfSquares() { // thetaHat = C * Beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); // thetaHat - thetaNull. Multiple by negative one to do subtraction RealMatrix thetaDiff = thetaHat.subtract(thetaNull); // calculate the hypothesis sum of squares: (thetaHat - thetaNull)'[C(X'X)-1C'](thetaHat - thetaNull) RealMatrix hss = thetaDiff.transpose().multiply(M.multiply(thetaDiff)); return hss; }
From source file:com.bolatu.gezkoncsvlogger.GyroOrientation.RotationKalmanFilter.java
/** * Correct the current state estimate with an actual measurement. * * @param z// ww w.j a v a 2 s .c o m * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance).multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 // instead of calculating the inverse of S we can rearrange the formula, // and then solve the linear equation A x X = B with A = S', X = K' and // B = (H * P)' // K(k) * S = P(k)- * H' // S' * K(k)' = H * P(k)-' RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver() .solve(measurementMatrix.multiply(errorCovariance.transpose())).transpose(); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
From source file:com.joptimizer.algebra.CholeskySparseFactorizationTest.java
public void testSimple1_norescaling() throws Exception { log.debug("testSimple1_norescaling"); double[][] A = new double[][] { { 4, 0, 2, 2 }, { 0, 4, 2, -2 }, { 2, 2, 6, 0 }, { 2, -2, 0, 6 } }; //expected L//from w w w . ja v a 2 s. c o m double[][] EL = new double[][] { { 2, 0, 0, 0 }, { 0, 2, 0, 0 }, { 1, 1, 2, 0 }, { 1, -1, 0, 2 } }; SparseDoubleMatrix2D ASparse = new SparseDoubleMatrix2D(A); CholeskySparseFactorization cs = new CholeskySparseFactorization(ASparse); cs.factorize(); DoubleMatrix2D L = cs.getL(); DoubleMatrix2D LT = cs.getLT(); log.debug("L : " + ArrayUtils.toString(L.toArray())); log.debug("LT: " + ArrayUtils.toString(LT.toArray())); RealMatrix ELMatrix = MatrixUtils.createRealMatrix(EL); RealMatrix LMatrix = MatrixUtils.createRealMatrix(L.toArray()); RealMatrix LTMatrix = MatrixUtils.createRealMatrix(LT.toArray()); assertEquals((ELMatrix.subtract(LMatrix).getNorm()), 0.); assertEquals((ELMatrix.subtract(LTMatrix.transpose()).getNorm()), 0.); //A.x = b double[] b = new double[] { 1, 2, 3, 4 }; double[] x = cs.solve(F1.make(b)).toArray(); //check the norm ||A.x-b|| double norm = new Array2DRowRealMatrix(A).operate(new ArrayRealVector(x)).subtract(new ArrayRealVector(b)) .getNorm(); log.debug("norm: " + norm); assertEquals(0., norm, Utils.getDoubleMachineEpsilon()); //check the scaled residual double residual = Utils.calculateScaledResidual(A, x, b); log.debug("residual: " + residual); assertEquals(0., residual, Utils.getDoubleMachineEpsilon()); }
From source file:com.joptimizer.algebra.CholeskySparseFactorizationTest.java
public void testSimple1_rescaling() throws Exception { log.debug("testSimple1_rescaling"); double[][] A = new double[][] { { 4, 0, 2, 2 }, { 0, 4, 2, -2 }, { 2, 2, 6, 0 }, { 2, -2, 0, 6 } }; //expected L// w w w. j a va2s. co m double[][] EL = new double[][] { { 2, 0, 0, 0 }, { 0, 2, 0, 0 }, { 1, 1, 2, 0 }, { 1, -1, 0, 2 } }; SparseDoubleMatrix2D ASparse = new SparseDoubleMatrix2D(A); CholeskySparseFactorization cs = new CholeskySparseFactorization(ASparse, new Matrix1NormRescaler()); cs.factorize(); DoubleMatrix2D L = cs.getL(); DoubleMatrix2D LT = cs.getLT(); log.debug("L : " + ArrayUtils.toString(L.toArray())); log.debug("LT: " + ArrayUtils.toString(LT.toArray())); RealMatrix ELMatrix = MatrixUtils.createRealMatrix(EL); RealMatrix LMatrix = MatrixUtils.createRealMatrix(L.toArray()); RealMatrix LTMatrix = MatrixUtils.createRealMatrix(LT.toArray()); assertEquals((ELMatrix.subtract(LMatrix).getNorm()), 0., Utils.getDoubleMachineEpsilon()); assertEquals((ELMatrix.subtract(LTMatrix.transpose()).getNorm()), 0., Utils.getDoubleMachineEpsilon()); //A.x = b double[] b = new double[] { 1, 2, 3, 4 }; double[] x = cs.solve(F1.make(b)).toArray(); //check the norm ||A.x-b|| double norm = new Array2DRowRealMatrix(A).operate(new ArrayRealVector(x)).subtract(new ArrayRealVector(b)) .getNorm(); log.debug("norm: " + norm); assertEquals(0., norm, 1.e-12); //check the scaled residual double residual = Utils.calculateScaledResidual(A, x, b); log.debug("residual: " + residual); assertEquals(0., residual, Utils.getDoubleMachineEpsilon()); }
From source file:com.analog.lyric.dimple.test.solvers.gibbs.TestGibbsRolledUp.java
/** * Adapted from MATLAB algoRolledUpGraphs/testGibbsTableFactor.m *///from w ww . j a va 2 s . c om @Test public void testGibbsTableFactor() { final int N = 100; final int bufferSize = 1; // Create graph Bit xi = name("xi", new Bit()); Bit xo = name("xo", new Bit()); FactorGraph sg = new FactorGraph(xi, xo); IFactorTable table = FactorTable.create(DiscreteDomain.bit(), DiscreteDomain.bit()); table.setWeightsDense(new double[] { 0, 1, 1, 0 }); sg.addFactor(table, xi, xo); FactorGraph fg = new FactorGraph(); BitStream x = new BitStream("x"); fg.addRepeatedFactorWithBufferSize(sg, bufferSize, x, x.getSlice(1)); // Generate data final double[][] input = new double[N][]; double val = 1.0;//testRand.nextBoolean() ? 1 : 0; input[0] = new double[] { val, 1 - val }; for (int i = 1; i < N; ++i) { double p = testRand.nextDouble(); val = p > table.getWeightForIndices((int) input[i - 1][0], 0) ? 1 : 0; input[i] = new double[] { val, 1 - val }; } // Solve using sum-product SumProductSolverGraph sfg1 = requireNonNull(fg.setSolverFactory(new SumProductSolver())); x.setDataSource(new DoubleArrayDataSource(input)); DoubleArrayDataSink sink1 = new DoubleArrayDataSink(); x.setDataSink(sink1); sfg1.solve(); // Solve again using Gibbs GibbsSolverGraph sfg2 = requireNonNull(fg.setSolverFactory(new GibbsSolver())); x.setDataSource(new DoubleArrayDataSource(input)); DoubleArrayDataSink sink2 = new DoubleArrayDataSink(); x.setDataSink(sink2); sfg2.solve(); RealMatrix b1 = wrapRealMatrix(sink1.getArray()); RealMatrix b2 = wrapRealMatrix(sink2.getArray()); // RealMatrixFormat fmt = new RealMatrixFormat("[","]","","", "; ", ","); // System.out.println(fmt.format(wrapRealMatrix(input))); // System.out.println(fmt.format(b1)); // System.out.println(fmt.format(b2)); assertEquals(0.0, b1.subtract(b2).getNorm(), 1e-20); }
From source file:edu.cudenver.bios.power.glmm.GLMMTestUnivariateRepeatedMeasures.java
/** * Ensure that the within subject contrast is orthonormal for all * UNIREP tests//from w w w. ja v a2 s . c o m */ protected void createOrthonormalU() { RealMatrix UtU = U.transpose().multiply(U); double upperLeft = UtU.getEntry(0, 0); if (upperLeft != 0) UtU = UtU.scalarMultiply(1 / upperLeft); RealMatrix diffFromIdentity = UtU.subtract( org.apache.commons.math3.linear.MatrixUtils.createRealIdentityMatrix(UtU.getRowDimension())); // get maximum absolute value in U'U double maxValue = Double.NEGATIVE_INFINITY; for (int r = 0; r < diffFromIdentity.getRowDimension(); r++) { for (int c = 0; c < diffFromIdentity.getColumnDimension(); c++) { double entryVal = Math.abs(diffFromIdentity.getEntry(r, c)); if (entryVal > maxValue) maxValue = entryVal; } } if (maxValue > Precision.SAFE_MIN) { // U'U matrix deviates from identity, so create a U matrix that is orthonormal // TODO: thus UNIREP tests may use a different U matrix than HLT/PBT/WLR tests??? // TODO: displayed matrix results are incorrect now? U = new GramSchmidtOrthonormalization(U).getQ(); debug("U replaced by orthonormal", U); } }
From source file:com.itemanalysis.psychometrics.factoranalysis.ObliminCriteria.java
public void computeValues(RealMatrix L) { final int k = L.getColumnDimension(); final int p = L.getRowDimension(); RealMatrix I = new IdentityMatrix(p); RealMatrix L2 = MatrixUtils.multiplyElements(L, L); RealMatrix N = new Array2DRowRealMatrix(k, k); N.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override// ww w . jav a2s. co m public double visit(int row, int column, double value) { if (row == column) return 0.0; return 1.0; } }); RealMatrix C = new Array2DRowRealMatrix(p, p); C.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return gam / (double) p; } }); RealMatrix X = I.subtract(C).multiply(L2).multiply(N); double sum = MatrixUtils.sumMatrix(MatrixUtils.multiplyElements(L2, X)); functionValue = sum / 4.0; gradient = MatrixUtils.multiplyElements(L, X); }
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. */// w w w . j a va 2 s . com private void initialize(Test test, RealMatrix FEssence, RealMatrix FtFinverse, int perGroupN, FixedRandomMatrix CFixedRand, RealMatrix U, RealMatrix thetaNull, RealMatrix beta, RealMatrix sigmaError, RealMatrix sigmaG, boolean exact) throws PowerException { debug("entering initialize"); // reset member variables this.T1 = null; this.FT1 = null; this.S = null; this.mzSq = null; this.H0 = 0; this.sStar = 0; // cache inputs this.test = test; this.FEssence = FEssence; this.FtFinverse = FtFinverse; this.perGroupN = perGroupN; this.CFixedRand = CFixedRand; this.U = U; this.thetaNull = thetaNull; this.beta = beta; this.sigmaError = sigmaError; this.sigmaG = sigmaG; // calculate intermediate matrices // RealMatrix FEssence = params.getDesignEssence().getFullDesignMatrixFixed(); // TODO: do we ever get here with values that can cause integer overflow, // and if so, does it matter? this.N = (double) FEssence.getRowDimension() * perGroupN; this.exact = exact; try { // TODO: need to calculate H0, need to adjust H1 for Unirep // get design matrix for fixed parameters only qF = FEssence.getColumnDimension(); // a = CFixedRand.getCombinedMatrix().getRowDimension(); // get fixed contrasts RealMatrix Cfixed = CFixedRand.getFixedMatrix(); RealMatrix CGaussian = CFixedRand.getRandomMatrix(); // build intermediate terms h1, S if (FtFinverse == null) { FtFinverse = new LUDecomposition(FEssence.transpose().multiply(FEssence)).getSolver().getInverse(); debug("FEssence", FEssence); debug("FtFinverse = (FEssence transpose * FEssence) inverse", FtFinverse); } else { debug("FtFinverse", FtFinverse); } RealMatrix PPt = Cfixed.multiply(FtFinverse.scalarMultiply(1 / (double) perGroupN)) .multiply(Cfixed.transpose()); debug("Cfixed", Cfixed); debug("n = " + perGroupN); debug("PPt = Cfixed * FtF inverse * (1/n) * Cfixed transpose", PPt); T1 = forceSymmetric(new LUDecomposition(PPt).getSolver().getInverse()); debug("T1 = PPt inverse", T1); FT1 = new CholeskyDecomposition(T1).getL(); debug("FT1 = Cholesky decomposition (L) of T1", FT1); // calculate theta difference // RealMatrix thetaNull = params.getTheta(); RealMatrix C = CFixedRand.getCombinedMatrix(); // RealMatrix beta = params.getScaledBeta(); // RealMatrix U = params.getWithinSubjectContrast(); // thetaHat = C * beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); debug("C", C); debug("beta", beta); debug("U", U); debug("thetaHat = C * beta * U", thetaHat); // thetaDiff = thetaHat - thetaNull RealMatrix thetaDiff = thetaHat.subtract(thetaNull); debug("thetaNull", thetaNull); debug("thetaDiff = thetaHat - thetaNull", thetaDiff); // TODO: specific to HLT or UNIREP RealMatrix sigmaStarInverse = getSigmaStarInverse(U, sigmaError, test); debug("sigmaStarInverse", sigmaStarInverse); RealMatrix H1matrix = thetaDiff.transpose().multiply(T1).multiply(thetaDiff).multiply(sigmaStarInverse); debug("H1matrix = thetaDiff transpose * T1 * thetaDiff * sigmaStarInverse", H1matrix); H1 = H1matrix.getTrace(); debug("H1 = " + H1); // Matrix which represents the non-centrality parameter as a linear combination of chi-squared r.v.'s. S = FT1.transpose().multiply(thetaDiff).multiply(sigmaStarInverse).multiply(thetaDiff.transpose()) .multiply(FT1).scalarMultiply(1 / H1); debug("S = FT1 transpose * thetaDiff * sigmaStar inverse * thetaDiff transpose * FT1 * (1/H1)", S); // We use the S matrix to generate the F-critical, numerical df's, and denominator df's // for a central F distribution. The resulting F distribution is used as an approximation // for the distribution of the non-centrality parameter. // See formulas 18-21 and A8,A10 from Glueck & Muller (2003) for details. EigenDecomposition sEigenDecomp = new EigenDecomposition(S); sEigenValues = sEigenDecomp.getRealEigenvalues(); // calculate H0 if (sEigenValues.length > 0) H0 = H1 * (1 - sEigenValues[0]); if (H0 <= 0) H0 = 0; // count the # of positive eigen values for (double value : sEigenValues) { if (value > 0) sStar++; } // TODO: throw error if sStar is <= 0 // TODO: NO: throw error if sStar != sEigenValues.length instead??? double stddevG = Math.sqrt(sigmaG.getEntry(0, 0)); RealMatrix svec = sEigenDecomp.getVT(); mzSq = svec.multiply(FT1.transpose()).multiply(CGaussian).scalarMultiply(1 / stddevG); for (int i = 0; i < mzSq.getRowDimension(); i++) { for (int j = 0; j < mzSq.getColumnDimension(); j++) { double entry = mzSq.getEntry(i, j); mzSq.setEntry(i, j, entry * entry); // TODO: is there an apache function to do this? } } debug("exiting initialize normally"); } catch (RuntimeException e) { LOGGER.warn("exiting initialize abnormally", e); throw new PowerException(e.getMessage(), PowerErrorEnum.INVALID_DISTRIBUTION_NONCENTRALITY_PARAMETER); } }
From source file:com.itemanalysis.psychometrics.factoranalysis.GPArotation.java
/** * Conducts orthogonal rotation of factor loadings. * * @param A matrix of orthogonal factor loadings * @return a matrix of rotated factor loadings. * @throws ConvergenceException//from w w w. j av a2 s. co m */ private RotationResults GPForth(RealMatrix A, boolean normalize, int maxIter, double eps) throws ConvergenceException { int ncol = A.getColumnDimension(); if (normalize) { //elementwise division by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value / W.getEntry(row, column); } }); } RealMatrix Tmat = new IdentityMatrix(ncol); double alpha = 1; RealMatrix L = A.multiply(Tmat); gpFunction.computeValues(L); double f = gpFunction.getValue(); RealMatrix VgQ = gpFunction.getGradient(); RealMatrix G = A.transpose().multiply(VgQ); double VgQtF = gpFunction.getValue(); RealMatrix VgQt = gpFunction.getGradient(); RealMatrix Tmatt = null; int iter = 0; double s = eps + 0.5; double s2 = 0; int innnerIter = 11; while (iter < maxIter) { RealMatrix M = Tmat.transpose().multiply(G); RealMatrix S = (M.add(M.transpose())); S = S.scalarMultiply(0.5); RealMatrix Gp = G.subtract(Tmat.multiply(S)); s = Math.sqrt((Gp.transpose().multiply(Gp)).getTrace()); s2 = Math.pow(s, 2); if (s < eps) break; alpha *= 2.0; for (int j = 0; j < innnerIter; j++) { Gp = Gp.scalarMultiply(alpha); RealMatrix X = (Tmat.subtract(Gp)); SingularValueDecomposition SVD = new SingularValueDecomposition(X); Tmatt = SVD.getU().multiply(SVD.getV().transpose()); L = A.multiply(Tmatt); gpFunction.computeValues(L); VgQt = gpFunction.getGradient(); VgQtF = gpFunction.getValue(); if (VgQtF < f - 0.5 * s2 * alpha) { break; } alpha /= 2.0; } Tmat = Tmatt; f = VgQtF; G = A.transpose().multiply(VgQt); iter++; } boolean convergence = s < eps; if (!convergence) { throw new ConvergenceException(); } if (normalize) { //elementwise multiplication by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value * W.getEntry(row, column); } }); } RealMatrix Phi = Tmat.transpose().multiply(Tmat); RotationResults result = new RotationResults(gpFunction.getValue(), L, Phi, Tmat, rotationMethod); return result; }