List of usage examples for org.apache.commons.math3.linear RealMatrix transpose
RealMatrix transpose();
From source file:org.knime.al.util.noveltydetection.knfst.OneClassKNFST.java
private NoveltyScores score(final RealMatrix kernelMatrix) { // projected test samples: final RealMatrix projectionVectors = kernelMatrix.transpose().multiply(m_projection); // differences to the target value: final RealMatrix diff = projectionVectors.subtract(MatrixFunctions .ones(kernelMatrix.getColumnDimension(), 1).scalarMultiply(m_targetPoints.getEntry(0, 0))); // distances to the target value: final RealVector scoresVector = MatrixFunctions .sqrt(MatrixFunctions.rowSums(MatrixFunctions.multiplyElementWise(diff, diff))); return new NoveltyScores(scoresVector.toArray(), projectionVectors); }
From source file:org.moeaframework.util.RotationMatrixBuilderTest.java
/** * Asserts that the matrix is a rotation matrix; that is, the matrix is * orthogonal and has a determinant of {@code 1}. * //from w w w . jav a2 s. c o m * @param rm the matrix to test */ public void testRotationMatrix(RealMatrix rm) { LUDecomposition lu = new LUDecomposition(rm); Assert.assertEquals(1.0, lu.getDeterminant(), Settings.EPS); TestUtils.assertEquals(rm.transpose(), lu.getSolver().getInverse(), new AbsoluteError(0.05)); }
From source file:org.ojalgo.benchmark.contestant.ACM.java
@Override public BenchmarkContestant<RealMatrix>.TransposedMultiplier getTransposedMultiplier() { return new TransposedMultiplier() { @Override// ww w .j av a2s. c om public RealMatrix apply(final RealMatrix left, final RealMatrix right) { return left.multiply(right.transpose()); } }; }
From source file:org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { // get the position in body frame final Transform fromBodyFrame = bodyFrame.getTransformTo(frame, date); final Transform toBodyFrame = fromBodyFrame.getInverse(); final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D()); // compute gradient and Hessian final GradientHessian gh = gradientHessian(date, positionBody); // gradient of the non-central part of the gravity field final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray(); // Hessian of the non-central part of the gravity field final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false); final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix()); final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot); // distribute all partial derivatives in a compact acceleration vector final int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); final double[] derivatives = new double[1 + parameters]; final DerivativeStructure[] accDer = new DerivativeStructure[3]; for (int i = 0; i < 3; ++i) { // first element is value of acceleration (i.e. gradient of field) derivatives[0] = gInertial[i];/*from w ww . ja v a 2 s .co m*/ // next three elements are one row of the Jacobian of acceleration (i.e. Hessian of field) derivatives[1] = hInertial.getEntry(i, 0); derivatives[2] = hInertial.getEntry(i, 1); derivatives[3] = hInertial.getEntry(i, 2); // next elements (three or four depending on mass being used or not) are left as 0 accDer[i] = new DerivativeStructure(parameters, order, derivatives); } return new FieldVector3D<DerivativeStructure>(accDer); }
From source file:org.pmad.gmm.MyMixMNDEM.java
/** * Fit a mixture model to the data supplied to the constructor. * * The quality of the fit depends on the concavity of the data provided to * the constructor and the initial mixture provided to this function. If the * data has many local optima, multiple runs of the fitting function with * different initial mixtures may be required to find the optimal solution. * If a SingularMatrixException is encountered, it is possible that another * initialization would work.// ww w . ja va2 s. co m * * @param initialMixture Model containing initial values of weights and * multivariate normals * @param maxIterations Maximum iterations allowed for fit * @param threshold Convergence threshold computed as difference in * logLikelihoods between successive iterations * @throws SingularMatrixException if any component's covariance matrix is * singular during fitting * @throws NotStrictlyPositiveException if numComponents is less than one * or threshold is less than Double.MIN_VALUE * @throws DimensionMismatchException if initialMixture mean vector and data * number of columns are not equal */ public void fit(final MyMixMND initialMixture, final int maxIterations, final double threshold) throws SingularMatrixException, NotStrictlyPositiveException, DimensionMismatchException { if (maxIterations < 1) { throw new NotStrictlyPositiveException(maxIterations); } if (threshold < Double.MIN_VALUE) { throw new NotStrictlyPositiveException(threshold); } final int n = data.length; // Number of data columns. Jagged data already rejected in constructor, // so we can assume the lengths of each row are equal. final int numCols = data[0].length; final int k = initialMixture.getComponents().size(); final int numMeanColumns = initialMixture.getComponents().get(0).getSecond().getMeans().length; if (numMeanColumns != numCols) { throw new DimensionMismatchException(numMeanColumns, numCols); } int numIterations = 0; double previousLogLikelihood = 0d; logLikelihood = Double.NEGATIVE_INFINITY; // Initialize model to fit to initial mixture. fittedModel = new MyMixMND(initialMixture.getComponents()); while (numIterations++ <= maxIterations && FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { System.out.println(numIterations); previousLogLikelihood = logLikelihood; double sumLogLikelihood = 0d; // Mixture components final List<Pair<Double, MyMND>> components = fittedModel.getComponents(); // Weight and distribution of each component final double[] weights = new double[k]; final MyMND[] mvns = new MyMND[k]; for (int j = 0; j < k; j++) { weights[j] = components.get(j).getFirst(); mvns[j] = components.get(j).getSecond(); } // E-step: compute the data dependent parameters of the expectation // function. // The percentage of row's total density between a row and a // component final double[][] gamma = new double[n][k]; // Sum of gamma for each component final double[] gammaSums = new double[k]; // for (double d : weights) { // System.out.print(d+" "); // } // System.out.println(); // Sum of gamma times its row for each each component final double[][] gammaDataProdSums = new double[k][numCols]; for (int i = 0; i < n; i++) { final double rowDensity = fittedModel.density(data[i]); sumLogLikelihood += FastMath.log(rowDensity); for (int j = 0; j < k; j++) { gamma[i][j] = weights[j] * mvns[j].density(data[i]); if (rowDensity == 0) { if (gamma[i][j] == 0) { gamma[i][j] = weights[j]; } else { System.err.println("ele"); } } else { gamma[i][j] /= rowDensity; } gammaSums[j] += gamma[i][j]; for (int col = 0; col < numCols; col++) { gammaDataProdSums[j][col] += gamma[i][j] * data[i][col]; } } } if (Utils.hasNaN(gamma)) { System.out.println("gamma has NaN"); } logLikelihood = sumLogLikelihood / n; // M-step: compute the new parameters based on the expectation // function. final double[] newWeights = new double[k]; final double[][] newMeans = new double[k][numCols]; for (int j = 0; j < k; j++) { newWeights[j] = gammaSums[j] / n; for (int col = 0; col < numCols; col++) { newMeans[j][col] = gammaDataProdSums[j][col] / gammaSums[j]; } } // Compute new covariance matrices final RealMatrix[] newCovMats = new RealMatrix[k]; for (int j = 0; j < k; j++) { newCovMats[j] = new Array2DRowRealMatrix(numCols, numCols); } for (int i = 0; i < n; i++) { for (int j = 0; j < k; j++) { final RealMatrix vec = new Array2DRowRealMatrix(MathArrays.ebeSubtract(data[i], newMeans[j])); final RealMatrix dataCov = vec.multiply(vec.transpose()).scalarMultiply(gamma[i][j]); newCovMats[j] = newCovMats[j].add(dataCov); } } // Converting to arrays for use by fitted model final double[][][] newCovMatArrays = new double[k][numCols][numCols]; for (int j = 0; j < k; j++) { newCovMats[j] = newCovMats[j].scalarMultiply(1d / gammaSums[j]); newCovMatArrays[j] = newCovMats[j].getData(); // if (Utils.hasNaN(newCovMatArrays[j])) { // System.out.println("covmet nan"); // } } // if (Utils.hasNaN(newMeans)) { // System.out.println("neams nan"); // } // Update current model fittedModel = new MyMixMND(newWeights, newMeans, newCovMatArrays); } if (FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { // Did not converge before the maximum number of iterations // throw new ConvergenceException(); System.out.println("Did not converge"); } }
From source file:org.rhwlab.dispim.nucleus.NucleusData.java
public RealMatrix quadraticSurfaceMatrix() { int D = 4;/*from w w w. j a v a2 s.c om*/ int Dm1 = 3; RealMatrix T = new Array2DRowRealMatrix(D, D); for (int r = 0; r < D; ++r) { for (int c = 0; c < D; ++c) { if (r == c) { T.setEntry(r, c, 1.0); } else { T.setEntry(r, c, 0.0); } } } T.setEntry(Dm1, 0, -xC); T.setEntry(Dm1, 1, -yC); T.setEntry(Dm1, 2, -zC); RealMatrix TT = T.transpose(); RealMatrix C = new Array2DRowRealMatrix(D, D); for (int r = 0; r < adjustedA.getRowDimension(); ++r) { for (int c = 0; c < adjustedA.getColumnDimension(); ++c) { C.setEntry(r, c, .8 * adjustedA.getEntry(r, c)); // make the effective radii bigger to get separation of the nuclei } } for (int d = 0; d < Dm1; ++d) { C.setEntry(Dm1, d, 0.0); C.setEntry(d, Dm1, 0.0); } C.setEntry(Dm1, Dm1, -1.0); return T.multiply(C.multiply(TT)); }
From source file:org.ssascaling.model.timeseries.learner.apache.OLSMultipleLinearRegression.java
/** * <p>Compute the "hat" matrix./* w ww . j av a 2 s . com*/ * </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()); }
From source file:org.ssascaling.model.timeseries.learner.apache.OLSMultipleLinearRegression.java
/** * <p>Calculates the variance-covariance matrix of the regression parameters. * </p>//from w w w. j a va 2 s .c o m * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup> * </p> * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup> * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of * R included, where p = the length of the beta vector.</p> * * @return The beta variance-covariance matrix */ @Override protected RealMatrix calculateBetaVariance() { int p = getX().getColumnDimension(); RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1, 0, p - 1); RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse(); return Rinv.multiply(Rinv.transpose()); }
From source file:org.ssascaling.model.timeseries.learner.apache.QRDecomposition.java
/** * Calculates the QR-decomposition of the given matrix. * * @param matrix The matrix to decompose. * @param threshold Singularity threshold. *///from w ww .j a v a 2 s. c o m public QRDecomposition(RealMatrix matrix, double threshold) { this.threshold = threshold; final int m = matrix.getRowDimension(); final int n = matrix.getColumnDimension(); qrt = matrix.transpose().getData(); rDiag = new double[FastMath.min(m, n)]; cachedQ = null; cachedQT = null; cachedR = null; cachedH = null; /* * The QR decomposition of a matrix A is calculated using Householder * reflectors by repeating the following operations to each minor * A(minor,minor) of A: */ for (int minor = 0; minor < FastMath.min(m, n); minor++) { final double[] qrtMinor = qrt[minor]; /* * Let x be the first column of the minor, and a^2 = |x|^2. * x will be in the positions qr[minor][minor] through qr[m][minor]. * The first column of the transformed minor will be (a,0,0,..)' * The sign of a is chosen to be opposite to the sign of the first * component of x. Let's find a: */ double xNormSqr = 0; for (int row = minor; row < m; row++) { final double c = qrtMinor[row]; xNormSqr += c * c; } final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); rDiag[minor] = a; if (a != 0.0) { /* * Calculate the normalized reflection vector v and transform * the first column. We know the norm of v beforehand: v = x-ae * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> = * a^2+a^2-2a<x,e> = 2a*(a - <x,e>). * Here <x, e> is now qr[minor][minor]. * v = x-ae is stored in the column at qr: */ qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor]) /* * Transform the rest of the columns of the minor: * They will be transformed by the matrix H = I-2vv'/|v|^2. * If x is a column vector of the minor, then * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v. * Therefore the transformation is easily calculated by * subtracting the column vector (2<x,v>/|v|^2)v from x. * * Let 2<x,v>/|v|^2 = alpha. From above we have * |v|^2 = -2a*(qr[minor][minor]), so * alpha = -<x,v>/(a*qr[minor][minor]) */ for (int col = minor + 1; col < n; col++) { final double[] qrtCol = qrt[col]; double alpha = 0; for (int row = minor; row < m; row++) { alpha -= qrtCol[row] * qrtMinor[row]; } alpha /= a * qrtMinor[minor]; // Subtract the column vector alpha*v from x. for (int row = minor; row < m; row++) { qrtCol[row] -= alpha * qrtMinor[row]; } } } } }
From source file:org.wso2.extension.siddhi.execution.var.models.parametric.ParametricVaRCalculator.java
/** * @return the var of the portfolio Calculate var for the given portfolio using updated covariances and means */// w w w . j ava 2s . c om @Override public Double processData(Portfolio portfolio, Event event) { Asset asset = getAssetPool().get(event.getSymbol()); if (asset.getNumberOfReturnValues() > 1) { RealMatrix varCovarMatrix = new Array2DRowRealMatrix(getVarCovarMatrix(portfolio)); RealMatrix weightageMatrix = new Array2DRowRealMatrix(getWeightageMatrix(portfolio)); RealMatrix meanMatrix = new Array2DRowRealMatrix(getMeanMatrix(portfolio)); RealMatrix portfolioVarianceMatrix = weightageMatrix.multiply(varCovarMatrix) .multiply(weightageMatrix.transpose()); RealMatrix portfolioMeanMatrix = weightageMatrix.multiply(meanMatrix.transpose()); double portfolioVariance = portfolioVarianceMatrix.getData()[0][0]; double portfolioMean = portfolioMeanMatrix.getData()[0][0]; if (portfolioVariance == 0) { // a normal distribution cannot be defined when sd = 0 return null; } double portfolioStandardDeviation = Math.sqrt(portfolioVariance); NormalDistribution normalDistribution = new NormalDistribution(portfolioMean, portfolioStandardDeviation); double zValue = normalDistribution.inverseCumulativeProbability(1 - getConfidenceInterval()); double var = zValue * portfolio.getTotalPortfolioValue(); return var; } return null; }