List of usage examples for org.apache.commons.math3.linear RealMatrix add
RealMatrix add(RealMatrix m) throws MatrixDimensionMismatchException;
From source file:hivemall.utils.math.StatsUtils.java
/** * @param mu1 mean vector of the first normal distribution * @param sigma1 covariance matrix of the first normal distribution * @param mu2 mean vector of the second normal distribution * @param sigma2 covariance matrix of the second normal distribution * @return the Hellinger distance between two multivariate normal distributions * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples *//*from www.j a v a2 s . c o m*/ public static double hellingerDistance(@Nonnull final RealVector mu1, @Nonnull final RealMatrix sigma1, @Nonnull final RealVector mu2, @Nonnull final RealMatrix sigma2) { RealVector muSub = mu1.subtract(mu2); RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d); LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean); double denominator = Math.sqrt(LUsigmaMean.getDeterminant()); if (denominator == 0.d) { return 1.d; // avoid divide by zero } RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0 double sigma1Det = MatrixUtils.det(sigma1); double sigma2Det = MatrixUtils.det(sigma2); double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d) * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub)); return 1.d - numerator / denominator; }
From source file:ch.zhaw.iamp.rct.weights.Weights.java
private static RealMatrix addNoise(final RealMatrix A) { RealMatrix buffer = A.copy(); RealMatrix noise = MatrixUtils.createRealMatrix(A.getRowDimension(), A.getColumnDimension()); for (int i = 0; i < A.getRowDimension(); ++i) { for (int j = 0; j < A.getColumnDimension(); ++j) { double noiseSign = Math.random() > 0.5 ? 1 : -1; double noiseAmplitude = 1; noise.setEntry(i, j, noiseSign * Math.random() * noiseAmplitude); }/*w ww . j a v a 2 s .c o m*/ } buffer = buffer.add(noise); return buffer; }
From source file:com.yahoo.egads.utilities.SpectralMethods.java
public static RealMatrix mFilter(RealMatrix data, int windowSize, FilteringMethod method, double methodParameter) { int n = data.getRowDimension(); int m = data.getColumnDimension(); int k = n - windowSize + 1; int i = 0, ind = 0; double[] temp; double sum = 0; RealMatrix hankelMat = SpectralMethods.createHankelMatrix(data, windowSize); SingularValueDecomposition svd = new SingularValueDecomposition(hankelMat); double[] singularValues = svd.getSingularValues(); switch (method) { case VARIANCE: temp = new double[singularValues.length - 1]; for (i = 1; i < singularValues.length; ++i) { sum += (singularValues[i] * singularValues[i]); }// w w w. j av a 2 s . com for (i = 0; i < temp.length; ++i) { temp[i] = (singularValues[i + 1] * singularValues[i + 1]) / sum; } sum = 0; for (i = temp.length - 1; i >= 0; --i) { sum += temp[i]; if (sum >= 1 - methodParameter) { ind = i; break; } } break; case EXPLICIT: ind = (int) Math.max(Math.min(methodParameter - 1, singularValues.length - 1), 0); break; case K_GAP: final double[] eigenGaps = new double[singularValues.length - 1]; Integer[] index = new Integer[singularValues.length - 1]; for (i = 0; i < eigenGaps.length; ++i) { eigenGaps[i] = singularValues[i] - singularValues[i + 1]; index[i] = i; } Arrays.sort(index, new Comparator<Integer>() { @Override public int compare(Integer o1, Integer o2) { return Double.compare(eigenGaps[o1], eigenGaps[o2]); } }); int maxIndex = 0; for (i = index.length - (int) methodParameter; i < index.length; ++i) { if (index[i] > maxIndex) { maxIndex = index[i]; } } ind = Math.min(maxIndex, singularValues.length / 3); break; case SMOOTHNESS: double[] variances = new double[singularValues.length]; for (i = 1; i < singularValues.length; ++i) { variances[i] = (singularValues[i] * singularValues[i]); } variances[0] = variances[1]; double smoothness = SpectralMethods .computeSmoothness(Arrays.copyOfRange(variances, 1, variances.length)); if (methodParameter - smoothness < 0.01) { methodParameter += 0.01; } double invalidS = smoothness; int validIndex = 1, invalidIndex = singularValues.length; while (true) { if (invalidS >= methodParameter) { ind = invalidIndex - 1; break; } else if (invalidIndex - validIndex <= 1) { ind = validIndex - 1; break; } int ii = (validIndex + invalidIndex) / 2; double[] tempVariances = Arrays.copyOf(Arrays.copyOfRange(variances, 0, ii + 1), singularValues.length); double s = SpectralMethods.computeSmoothness(tempVariances); if (s >= methodParameter) { validIndex = ii; } else { invalidIndex = ii; invalidS = s; } } break; case EIGEN_RATIO: int startIndex = 0, endIndex = singularValues.length - 1; if (singularValues[endIndex] / singularValues[0] >= methodParameter) { ind = endIndex; } else { while (true) { int midIndex = (startIndex + endIndex) / 2; if (singularValues[midIndex] / singularValues[0] >= methodParameter) { if (singularValues[midIndex + 1] / singularValues[0] < methodParameter) { ind = midIndex; break; } else { startIndex = midIndex; } } else { endIndex = midIndex; } } } break; case GAP_RATIO: double[] gaps = new double[singularValues.length - 1]; for (i = 0; i < gaps.length; ++i) { gaps[i] = singularValues[i] - singularValues[i + 1]; } ind = 0; for (i = gaps.length - 1; i >= 0; --i) { if (gaps[i] / singularValues[0] >= methodParameter) { ind = i; break; } } break; default: ind = singularValues.length - 1; break; } ind = Math.max(0, Math.min(ind, singularValues.length - 1)); RealMatrix truncatedHankelMatrix = MatrixUtils.createRealMatrix(k, m * windowSize); RealMatrix mU = svd.getU(); RealMatrix mVT = svd.getVT(); for (i = 0; i <= ind; ++i) { truncatedHankelMatrix = truncatedHankelMatrix .add(mU.getColumnMatrix(i).multiply(mVT.getRowMatrix(i)).scalarMultiply(singularValues[i])); } return SpectralMethods.averageHankelMatrix(truncatedHankelMatrix, windowSize); }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 618" */// www . j av a 2 s . c om private RealMatrix buildS(double[] X) { RealMatrix S = this.G.scalarMultiply(-1); for (int i = 0; i < dim; i++) { S = S.add(this.Fi[i].scalarMultiply(-1 * X[i])); } return S; }
From source file:com.analog.lyric.dimple.test.solvers.sumproduct.TestSampledFactors.java
/** * Generates a random covariance matrix with given dimension. *//*w w w .j a va 2 s . c o m*/ RealMatrix randCovariance(int n) { RealMatrix A = MatrixUtils.createRealMatrix(n, n); // Randomize A.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return testRand.nextDouble(); } }); RealMatrix B = A.add(A.transpose()); // B is symmetric double minEig = Doubles.min(new EigenDecomposition(B).getRealEigenvalues()); double r = testRand.nextGaussian(); r *= r; r += Math.ulp(1.0); RealMatrix I = MatrixUtils.createRealIdentityMatrix(n); RealMatrix C = B.add(I.scalarMultiply(r - minEig)); return C; }
From source file:edu.cudenver.bios.power.glmm.GLMMTestWilksLambda.java
/** * Compute a Wilks Lamba statistic//from w w w .j a v a 2 s. c om * * @param H hypothesis sum of squares matrix * @param E error sum of squares matrix * @returns F statistic */ private double getWilksLambda(RealMatrix H, RealMatrix E, DistributionType type) throws IllegalArgumentException { if (!H.isSquare() || !E.isSquare() || H.getColumnDimension() != E.getRowDimension()) throw new IllegalArgumentException( "Failed to compute Wilks Lambda: hypothesis and error matrices must be square and same dimensions"); double a = C.getRowDimension(); double b = U.getColumnDimension(); double s = (a < b) ? a : b; double p = beta.getColumnDimension(); RealMatrix adjustedH = H; if (type != DistributionType.DATA_ANALYSIS_NULL && ((s == 1 && p > 1) || fMethod == FApproximation.RAO_TWO_MOMENT_OMEGA_MULT)) { adjustedH = H.scalarMultiply((totalN - rank) / totalN); } RealMatrix T = adjustedH.add(E); RealMatrix inverseT = new LUDecomposition(T).getSolver().getInverse(); RealMatrix EinverseT = E.multiply(inverseT); double lambda = new LUDecomposition(EinverseT).getDeterminant(); return lambda; }
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 a v a 2 s. c om */ 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; }
From source file:edu.cudenver.bios.power.glmm.GLMMTestPillaiBartlett.java
/** * Compute a Pillai Bartlett Trace statistic * //from w ww.j a va 2 s .c o m * @param H hypothesis sum of squares matrix * @param E error sum of squares matrix * @returns F statistic */ private double getPillaiBartlettTrace(RealMatrix H, RealMatrix E) { if (!H.isSquare() || !E.isSquare() || H.getColumnDimension() != E.getRowDimension()) throw new IllegalArgumentException( "Failed to compute Pillai-Bartlett Trace: hypothesis and error matrices must be square and same dimensions"); double a = C.getRowDimension(); double b = U.getColumnDimension(); double s = (a < b) ? a : b; double p = beta.getColumnDimension(); RealMatrix adjustedH = H; if ((s == 1 && p > 1) || fMethod == FApproximation.PILLAI_ONE_MOMENT_OMEGA_MULT || fMethod == FApproximation.MULLER_TWO_MOMENT_OMEGA_MULT) { adjustedH = H.scalarMultiply(((double) (totalN - rank) / (double) totalN)); } RealMatrix T = adjustedH.add(E); RealMatrix inverseT = new LUDecomposition(T).getSolver().getInverse(); RealMatrix HinverseT = adjustedH.multiply(inverseT); return HinverseT.getTrace(); }
From source file:edu.oregonstate.eecs.mcplan.ml.LinearDiscriminantAnalysis.java
/** * @param data The elements of 'data' will be modified. * @param label//from ww w. ja v a 2s .co m * @param Nclasses * @param shrinkage_intensity */ public LinearDiscriminantAnalysis(final ArrayList<double[]> data, final int[] label, final int Nclasses, final double shrinkage) { assert (data.size() == label.length); final int Ndata = data.size(); final int Ndim = data.get(0).length; // Partition data by class final ArrayList<ArrayList<double[]>> classes = new ArrayList<ArrayList<double[]>>(Nclasses); for (int i = 0; i < Nclasses; ++i) { classes.add(new ArrayList<double[]>()); } for (int i = 0; i < data.size(); ++i) { classes.get(label[i]).add(data.get(i)); } // Mean center the data final VectorMeanVarianceAccumulator mv = new VectorMeanVarianceAccumulator(Ndim); for (int i = 0; i < Ndata; ++i) { mv.add(data.get(i)); } mean = mv.mean(); // Subtract global mean for (final double[] x : data) { Fn.vminus_inplace(x, mean); } // Calculate class means and covariances final double[][] class_mean = new double[Nclasses][Ndim]; final RealMatrix[] class_cov = new RealMatrix[Nclasses]; for (int i = 0; i < Nclasses; ++i) { final ArrayList<double[]> Xc = classes.get(i); final VectorMeanVarianceAccumulator mv_i = new VectorMeanVarianceAccumulator(Ndim); final StorelessCovariance cov = new StorelessCovariance(Ndim); for (int j = 0; j < Xc.size(); ++j) { final double[] x = Xc.get(j); mv_i.add(x); cov.increment(x); } class_mean[i] = mv_i.mean(); class_cov[i] = cov.getCovarianceMatrix(); } // Between-class scatter. // Note that 'data' is mean-centered, so the global mean is 0. RealMatrix Sb_builder = new Array2DRowRealMatrix(Ndim, Ndim); for (int i = 0; i < Nclasses; ++i) { final RealVector mu_i = new ArrayRealVector(class_mean[i]); final RealMatrix xxt = mu_i.outerProduct(mu_i); Sb_builder = Sb_builder.add(xxt.scalarMultiply(classes.get(i).size() / ((double) Ndata - 1))); } this.Sb = Sb_builder; Sb_builder = null; // Within-class scatter with shrinkage estimate: // Sw = (1.0 - shrinkage) * \sum Sigma_i + shrinkage * I RealMatrix Sw_builder = new Array2DRowRealMatrix(Ndim, Ndim); for (int i = 0; i < Nclasses; ++i) { final RealMatrix Sigma_i = class_cov[i]; final RealMatrix scaled = Sigma_i.scalarMultiply((1.0 - shrinkage) * (classes.get(i).size() - 1)); Sw_builder = Sw_builder.add(scaled); } for (int i = 0; i < Ndim; ++i) { Sw_builder.setEntry(i, i, Sw_builder.getEntry(i, i) + shrinkage); } this.Sw = Sw_builder; Sw_builder = null; // Invert Sw System.out.println("[LDA] Sw inverse"); final RealMatrix Sw_inv = new LUDecomposition(Sw).getSolver().getInverse(); final RealMatrix F = Sw_inv.multiply(Sb); System.out.println("[LDA] Eigendecomposition"); eigenvalues = new double[Nclasses - 1]; eigenvectors = new ArrayList<RealVector>(Nclasses - 1); final EigenDecomposition evd = new EigenDecomposition(F); for (int j = 0; j < Nclasses - 1; ++j) { final double eigenvalue = evd.getRealEigenvalue(j); eigenvalues[j] = eigenvalue; // final double scale = 1.0 / Math.sqrt( eigenvalue ); // eigenvectors.add( evd.getEigenvector( j ).mapMultiply( scale ) ); eigenvectors.add(evd.getEigenvector(j)); } }
From source file:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java
private void init() { final int step = Math.max(1, n_ / k_); final double unif = 1.0 / k_; double acc = 0.0; final RandomPermutationIterator<double[]> r = new RandomPermutationIterator<double[]>(data_, rng_); final RandomPermutationIterator<double[]> rrepeat = new RandomPermutationIterator<double[]>(data_, r.permutation());//from w ww. j ava 2 s .c o m for (int i = 0; i < k_; ++i) { final RealVector mu = new ArrayRealVector(d_); for (int j = 0; j < step; ++j) { final double[] x = r.next(); final RealVector v = new ArrayRealVector(x); mu.combineToSelf(1.0, 1.0, v); } final double Zinv = 1.0 / step; mu.mapMultiplyToSelf(Zinv); RealMatrix Sigma = new Array2DRowRealMatrix(d_, d_); for (int j = 0; j < step; ++j) { final double[] x = rrepeat.next(); final RealVector v = new ArrayRealVector(x); v.combineToSelf(1.0, -1.0, mu); Sigma = Sigma.add(v.outerProduct(v)); } Sigma = Sigma.scalarMultiply(Zinv); pi_[i] = unif; acc += unif; mu_[i] = mu; Sigma_[i] = Sigma; //MatrixUtils.createRealIdentityMatrix( d_ ); } pi_[k_ - 1] += (1.0 - acc); // Round-off error }