List of usage examples for org.apache.commons.math3.linear EigenDecomposition getVT
public RealMatrix getVT()
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. *///from ww w . j a v a 2 s. c o m 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:org.ojalgo.benchmark.lab.library.ACM.java
@Override public DecompositionOperation<RealMatrix, RealMatrix> getOperationEvD(final int dim) { final RealMatrix[] ret = this.makeArray(3); return (matrix) -> { final EigenDecomposition svd = new EigenDecomposition(matrix); ret[0] = svd.getV();/*from ww w. ja v a 2 s . c o m*/ ret[1] = svd.getD(); ret[2] = svd.getVT(); return ret; }; }
From source file:playground.sergioo.facilitiesGenerator2012.WorkFacilitiesGeneration.java
private static Set<PointPerson> getPCATransformation(Collection<PointPerson> points) { RealMatrix pointsM = new Array2DRowRealMatrix(points.iterator().next().getDimension(), points.size()); int k = 0;/*from ww w. j a v a 2 s .c o m*/ for (PointND<Double> point : points) { for (int f = 0; f < point.getDimension(); f++) pointsM.setEntry(f, k, point.getElement(f)); k++; } RealMatrix means = new Array2DRowRealMatrix(pointsM.getRowDimension(), 1); for (int r = 0; r < means.getRowDimension(); r++) { double mean = 0; for (int c = 0; c < pointsM.getColumnDimension(); c++) mean += pointsM.getEntry(r, c) / pointsM.getColumnDimension(); means.setEntry(r, 0, mean); } RealMatrix deviations = new Array2DRowRealMatrix(pointsM.getRowDimension(), pointsM.getColumnDimension()); for (int r = 0; r < deviations.getRowDimension(); r++) for (int c = 0; c < deviations.getColumnDimension(); c++) deviations.setEntry(r, c, pointsM.getEntry(r, c) - means.getEntry(r, 0)); RealMatrix covariance = deviations.multiply(deviations.transpose()) .scalarMultiply(1 / (double) pointsM.getColumnDimension()); EigenDecomposition eigenDecomposition = new EigenDecomposition(covariance, 0); RealMatrix eigenVectorsT = eigenDecomposition.getVT(); RealVector eigenValues = new ArrayRealVector(eigenDecomposition.getD().getRowDimension()); for (int r = 0; r < eigenDecomposition.getD().getRowDimension(); r++) eigenValues.setEntry(r, eigenDecomposition.getD().getEntry(r, r)); for (int i = 0; i < eigenValues.getDimension(); i++) { for (int j = i + 1; j < eigenValues.getDimension(); j++) if (eigenValues.getEntry(i) < eigenValues.getEntry(j)) { double tempValue = eigenValues.getEntry(i); eigenValues.setEntry(i, eigenValues.getEntry(j)); eigenValues.setEntry(j, tempValue); RealVector tempVector = eigenVectorsT.getRowVector(i); eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(j)); eigenVectorsT.setRowVector(j, tempVector); } eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(i).mapMultiply(Math.sqrt(1 / eigenValues.getEntry(i)))); } RealVector standardDeviations = new ArrayRealVector(pointsM.getRowDimension()); for (int r = 0; r < covariance.getRowDimension(); r++) standardDeviations.setEntry(r, Math.sqrt(covariance.getEntry(r, r))); double zValue = standardDeviations.dotProduct(new ArrayRealVector(pointsM.getRowDimension(), 1)); RealMatrix zScore = deviations.scalarMultiply(1 / zValue); pointsM = eigenVectorsT.multiply(zScore); Set<PointPerson> pointsC = new HashSet<PointPerson>(); k = 0; for (PointPerson point : points) { PointPerson pointC = new PointPerson(point.getId(), point.getOccupation(), new Double[] { pointsM.getEntry(0, k), pointsM.getEntry(1, k) }, point.getPlaceType()); pointC.setWeight(point.getWeight()); pointsC.add(pointC); k++; } return pointsC; }