List of usage examples for org.apache.commons.math3.linear RealMatrix transpose
RealMatrix transpose();
From source file:com.datumbox.framework.machinelearning.featureselection.continuous.PCA.java
@Override protected void estimateModelParameters(Dataset originaldata) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = originaldata.size(); int d = originaldata.getColumnSize(); //convert data into matrix Map<Object, Integer> feature2ColumnId = modelParameters.getFeature2ColumnId(); MatrixDataset matrixDataset = MatrixDataset.newInstance(originaldata, false, feature2ColumnId); RealMatrix X = matrixDataset.getX(); //calculate means and subtract them from data double[] meanValues = new double[d]; for (Map.Entry<Object, Integer> entry : feature2ColumnId.entrySet()) { Object feature = entry.getKey(); Integer columnId = entry.getValue(); meanValues[columnId] = Descriptives .mean(originaldata.extractColumnValues(feature).toFlatDataCollection()); for (int row = 0; row < n; ++row) { X.addToEntry(row, columnId, -meanValues[columnId]); //inplace subtraction!!! }//www . ja v a 2 s . co m } modelParameters.setMean(meanValues); RealMatrix components; double[] eigenValues; /* if(d>n) { //turned off because of the algorithm could not be validated //Karhunen Lowe Transform to speed up calculations //nxn matrix RealMatrix covarianceNN = (X.multiply(X.transpose())).scalarMultiply(1.0/(n-1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceNN); eigenValues = decomposition.getRealEigenvalues(); RealMatrix eigenVectors = decomposition.getV(); double[] sqrtInverseEigenValues = new double[eigenValues.length]; for(int i=0;i<eigenValues.length;++i) { if(eigenValues[i]==0.0) { sqrtInverseEigenValues[i] = 0.0; } else { sqrtInverseEigenValues[i] = 1.0/Math.sqrt(eigenValues[i]); } } components = X.transpose().multiply(eigenVectors); //Components = X'*V*L^-0.5; To whiten them we multiply with L^0.5 which //cancels out the previous multiplication. So below we multiply by //L^-0.5 ONLY if we don't whiten. if(!knowledgeBase.getTrainingParameters().isWhitened()) { components = components.multiply(new DiagonalMatrix(sqrtInverseEigenValues)); } } else { //Normal PCA goes here } */ //dxd matrix RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0 / (n - 1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceDD); eigenValues = decomposition.getRealEigenvalues(); components = decomposition.getV(); //Whiten Components W = U*L^0.5; To whiten them we multiply with L^0.5. if (knowledgeBase.getTrainingParameters().isWhitened()) { double[] sqrtEigenValues = new double[eigenValues.length]; for (int i = 0; i < eigenValues.length; ++i) { sqrtEigenValues[i] = Math.sqrt(eigenValues[i]); } components = components.multiply(new DiagonalMatrix(sqrtEigenValues)); } //the eigenvalues and their components are sorted by descending order no need to resort them Integer maxDimensions = knowledgeBase.getTrainingParameters().getMaxDimensions(); Double varianceThreshold = knowledgeBase.getTrainingParameters().getVarianceThreshold(); if (varianceThreshold != null && varianceThreshold <= 1) { double sum = 0.0; double totalVariance = StatUtils.sum(eigenValues); int varCounter = 0; for (double l : eigenValues) { sum += l / totalVariance; ++varCounter; if (sum >= varianceThreshold) { break; } } if (maxDimensions == null || maxDimensions > varCounter) { maxDimensions = varCounter; } } if (maxDimensions != null && maxDimensions < d) { //keep only the maximum selected eigenvalues double[] newEigenValues = new double[maxDimensions]; System.arraycopy(eigenValues, 0, newEigenValues, 0, maxDimensions); eigenValues = newEigenValues; //keep only the maximum selected eigenvectors components = components.getSubMatrix(0, components.getRowDimension() - 1, 0, maxDimensions - 1); } modelParameters.setRows(components.getRowDimension()); modelParameters.setCols(components.getColumnDimension()); modelParameters.setEigenValues(eigenValues); modelParameters.setComponents(components.getData()); }
From source file:com.github.tteofili.looseen.yay.SGM.java
/** * perform weights learning from the training examples using (configurable) mini batch gradient descent algorithm * * @param samples the training examples//from w w w . j a v a 2 s .c o m * @return the final cost with the updated weights * @throws Exception if BGD fails to converge or any numerical error happens */ private double learnWeights(Sample... samples) throws Exception { int iterations = 0; double cost = Double.MAX_VALUE; int j = 0; // momentum RealMatrix vb = MatrixUtils.createRealMatrix(biases[0].getRowDimension(), biases[0].getColumnDimension()); RealMatrix vb2 = MatrixUtils.createRealMatrix(biases[1].getRowDimension(), biases[1].getColumnDimension()); RealMatrix vw = MatrixUtils.createRealMatrix(weights[0].getRowDimension(), weights[0].getColumnDimension()); RealMatrix vw2 = MatrixUtils.createRealMatrix(weights[1].getRowDimension(), weights[1].getColumnDimension()); long start = System.currentTimeMillis(); int c = 1; RealMatrix x = MatrixUtils.createRealMatrix(configuration.batchSize, samples[0].getInputs().length); RealMatrix y = MatrixUtils.createRealMatrix(configuration.batchSize, samples[0].getOutputs().length); while (true) { int i = 0; for (int k = j * configuration.batchSize; k < j * configuration.batchSize + configuration.batchSize; k++) { Sample sample = samples[k % samples.length]; x.setRow(i, sample.getInputs()); y.setRow(i, sample.getOutputs()); i++; } j++; long time = (System.currentTimeMillis() - start) / 1000; if (iterations % (1 + (configuration.maxIterations / 100)) == 0 && time > 60 * c) { c += 1; // System.out.println("cost: " + cost + ", accuracy: " + evaluate(this) + " after " + iterations + " iterations in " + (time / 60) + " minutes (" + ((double) iterations / time) + " ips)"); } RealMatrix w0t = weights[0].transpose(); RealMatrix w1t = weights[1].transpose(); RealMatrix hidden = rectifierFunction.applyMatrix(x.multiply(w0t)); hidden.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + biases[0].getEntry(0, column); } @Override public double end() { return 0; } }); RealMatrix scores = hidden.multiply(w1t); scores.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + biases[1].getEntry(0, column); } @Override public double end() { return 0; } }); RealMatrix probs = scores.copy(); int len = scores.getColumnDimension() - 1; for (int d = 0; d < configuration.window - 1; d++) { int startColumn = d * len / (configuration.window - 1); RealMatrix subMatrix = scores.getSubMatrix(0, scores.getRowDimension() - 1, startColumn, startColumn + x.getColumnDimension()); for (int sm = 0; sm < subMatrix.getRowDimension(); sm++) { probs.setSubMatrix(softmaxActivationFunction.applyMatrix(subMatrix.getRowMatrix(sm)).getData(), sm, startColumn); } } RealMatrix correctLogProbs = MatrixUtils.createRealMatrix(x.getRowDimension(), 1); correctLogProbs.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return -Math.log(probs.getEntry(row, getMaxIndex(y.getRow(row)))); } @Override public double end() { return 0; } }); double dataLoss = correctLogProbs.walkInOptimizedOrder(new RealMatrixPreservingVisitor() { private double d = 0; @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public void visit(int row, int column, double value) { d += value; } @Override public double end() { return d; } }) / samples.length; double reg = 0d; reg += weights[0].walkInOptimizedOrder(new RealMatrixPreservingVisitor() { private double d = 0d; @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public void visit(int row, int column, double value) { d += Math.pow(value, 2); } @Override public double end() { return d; } }); reg += weights[1].walkInOptimizedOrder(new RealMatrixPreservingVisitor() { private double d = 0d; @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public void visit(int row, int column, double value) { d += Math.pow(value, 2); } @Override public double end() { return d; } }); double regLoss = 0.5 * configuration.regularizationLambda * reg; double newCost = dataLoss + regLoss; if (iterations == 0) { // System.out.println("started with cost = " + dataLoss + " + " + regLoss + " = " + newCost); } if (Double.POSITIVE_INFINITY == newCost) { throw new Exception("failed to converge at iteration " + iterations + " with alpha " + configuration.alpha + " : cost going from " + cost + " to " + newCost); } else if (iterations > 1 && (newCost < configuration.threshold || iterations > configuration.maxIterations)) { cost = newCost; // System.out.println("successfully converged after " + (iterations - 1) + " iterations (alpha:" + configuration.alpha + ",threshold:" + configuration.threshold + ") with cost " + newCost); break; } else if (Double.isNaN(newCost)) { throw new Exception("failed to converge at iteration " + iterations + " with alpha " + configuration.alpha + " : cost calculation underflow"); } // update registered cost cost = newCost; // calculate the derivatives to update the parameters RealMatrix dscores = probs.copy(); dscores.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return (y.getEntry(row, column) == 1 ? (value - 1) : value) / samples.length; } @Override public double end() { return 0; } }); // get derivative on second layer RealMatrix dW2 = hidden.transpose().multiply(dscores); // regularize dw2 dW2.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + configuration.regularizationLambda * w1t.getEntry(row, column); } @Override public double end() { return 0; } }); RealMatrix db2 = MatrixUtils.createRealMatrix(biases[1].getRowDimension(), biases[1].getColumnDimension()); dscores.walkInOptimizedOrder(new RealMatrixPreservingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public void visit(int row, int column, double value) { db2.setEntry(0, column, db2.getEntry(0, column) + value); } @Override public double end() { return 0; } }); RealMatrix dhidden = dscores.multiply(weights[1]); dhidden.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value < 0 ? 0 : value; } @Override public double end() { return 0; } }); RealMatrix db = MatrixUtils.createRealMatrix(biases[0].getRowDimension(), biases[0].getColumnDimension()); dhidden.walkInOptimizedOrder(new RealMatrixPreservingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public void visit(int row, int column, double value) { db.setEntry(0, column, db.getEntry(0, column) + value); } @Override public double end() { return 0; } }); // get derivative on first layer RealMatrix dW = x.transpose().multiply(dhidden); // regularize dW.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + configuration.regularizationLambda * w0t.getEntry(row, column); } @Override public double end() { return 0; } }); RealMatrix dWt = dW.transpose(); RealMatrix dWt2 = dW2.transpose(); if (configuration.useNesterovMomentum) { // update nesterov momentum final RealMatrix vbPrev = vb.copy(); final RealMatrix vb2Prev = vb2.copy(); final RealMatrix vwPrev = vw.copy(); final RealMatrix vw2Prev = vw2.copy(); vb.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * db.getEntry(row, column); } @Override public double end() { return 0; } }); vb2.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * db2.getEntry(row, column); } @Override public double end() { return 0; } }); vw.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * dWt.getEntry(row, column); } @Override public double end() { return 0; } }); vw2.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * dWt2.getEntry(row, column); } @Override public double end() { return 0; } }); // update bias biases[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.mu * vbPrev.getEntry(row, column) + (1 + configuration.mu) * vb.getEntry(row, column); } @Override public double end() { return 0; } }); biases[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.mu * vb2Prev.getEntry(row, column) + (1 + configuration.mu) * vb2.getEntry(row, column); } @Override public double end() { return 0; } }); // update the weights weights[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.mu * vwPrev.getEntry(row, column) + (1 + configuration.mu) * vw.getEntry(row, column); } @Override public double end() { return 0; } }); weights[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.mu * vw2Prev.getEntry(row, column) + (1 + configuration.mu) * vw2.getEntry(row, column); } @Override public double end() { return 0; } }); } else if (configuration.useMomentum) { // update momentum vb.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * db.getEntry(row, column); } @Override public double end() { return 0; } }); vb2.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * db2.getEntry(row, column); } @Override public double end() { return 0; } }); vw.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * dWt.getEntry(row, column); } @Override public double end() { return 0; } }); vw2.walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return configuration.mu * value - configuration.alpha * dWt2.getEntry(row, column); } @Override public double end() { return 0; } }); // update bias biases[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + vb.getEntry(row, column); } @Override public double end() { return 0; } }); biases[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + vb2.getEntry(row, column); } @Override public double end() { return 0; } }); // update the weights weights[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + vw.getEntry(row, column); } @Override public double end() { return 0; } }); weights[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value + vw2.getEntry(row, column); } @Override public double end() { return 0; } }); } else { // standard parameter update // update bias biases[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.alpha * db.getEntry(row, column); } @Override public double end() { return 0; } }); biases[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.alpha * db2.getEntry(row, column); } @Override public double end() { return 0; } }); // update the weights weights[0].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.alpha * dWt.getEntry(row, column); } @Override public double end() { return 0; } }); weights[1].walkInOptimizedOrder(new RealMatrixChangingVisitor() { @Override public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) { } @Override public double visit(int row, int column, double value) { return value - configuration.alpha * dWt2.getEntry(row, column); } @Override public double end() { return 0; } }); } iterations++; } return cost; }
From source file:org.apache.solr.client.solrj.io.eval.DistanceEvaluator.java
private Matrix distance(DistanceMeasure distanceMeasure, Matrix matrix) { double[][] data = matrix.getData(); RealMatrix realMatrix = new Array2DRowRealMatrix(data); realMatrix = realMatrix.transpose(); data = realMatrix.getData();/*from w w w.j ava2 s.c o m*/ double[][] distanceMatrix = new double[data.length][data.length]; for (int i = 0; i < data.length; i++) { double[] row = data[i]; for (int j = 0; j < data.length; j++) { double[] row2 = data[j]; double dist = distanceMeasure.compute(row, row2); distanceMatrix[i][j] = dist; } } return new Matrix(distanceMatrix); }
From source file:org.apache.solr.client.solrj.io.eval.SumColumnsEvaluator.java
@Override public Object doWork(Object value) throws IOException { if (null == value) { return null; } else if (value instanceof Matrix) { //First transpose the matrix Matrix matrix = (Matrix) value;// www .ja va 2 s . c om double[][] data = matrix.getData(); RealMatrix realMatrix = new Array2DRowRealMatrix(data); realMatrix = realMatrix.transpose(); data = realMatrix.getData(); //Now sum the rows. List<Number> sums = new ArrayList(data.length); for (int i = 0; i < data.length; i++) { double sum = 0; double[] row = data[i]; for (int j = 0; j < row.length; j++) { sum += row[j]; } sums.add(sum); } return sums; } else { throw new IOException("Grand sum function only operates on a matrix"); } }
From source file:org.cirdles.calamari.algorithms.WeightedMeanCalculators.java
public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[] x, double[][] sigmaRhoY) { WeightedLinearCorrResults weightedLinearCorrResults = new WeightedLinearCorrResults(); RealMatrix omega = new BlockRealMatrix(convertCorrelationsToCovariances(sigmaRhoY)); RealMatrix invOmega = MatrixUtils.inverse(omega); int n = y.length; double mX = 0; double pX = 0; double pY = 0; double pXY = 0; double w = 0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { double invOm = invOmega.getEntry(i, j); w += invOm;/* ww w . j a va 2 s. co m*/ pX += (invOm * (x[i] + x[j])); pY += (invOm * (y[i] + y[j])); pXY += (invOm * (((x[i] * y[j]) + (x[j] * y[i])))); mX += (invOm * x[i] * x[j]); } } double slope = ((2 * pXY * w) - (pX * pY)) / ((4 * mX * w) - (pX * pX)); double intercept = (pY - (slope * pX)) / (2 * w); RealMatrix fischer = new BlockRealMatrix(new double[][] { { mX, pX / 2.0 }, { pX / 2.0, w } }); RealMatrix fischerInv = MatrixUtils.inverse(fischer); double slopeSig = StrictMath.sqrt(fischerInv.getEntry(0, 0)); double interceptSig = StrictMath.sqrt(fischerInv.getEntry(1, 1)); double slopeInterceptCov = fischerInv.getEntry(0, 1); RealMatrix resid = new BlockRealMatrix(n, 1); for (int i = 0; i < n; i++) { resid.setEntry(i, 0, y[i] - (slope * x[i]) - intercept); } RealMatrix residT = resid.transpose(); RealMatrix mM = residT.multiply(invOmega).multiply(resid); double sumSqWtdResids = mM.getEntry(0, 0); double mswd = sumSqWtdResids / (n - 2); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 2), 1E9); double prob = 1.0 - fdist.cumulativeProbability(mswd); weightedLinearCorrResults.setBad(false); weightedLinearCorrResults.setSlope(slope); weightedLinearCorrResults.setIntercept(intercept); weightedLinearCorrResults.setSlopeSig(slopeSig); weightedLinearCorrResults.setInterceptSig(interceptSig); weightedLinearCorrResults.setSlopeInterceptCov(slopeInterceptCov); weightedLinearCorrResults.setMswd(mswd); weightedLinearCorrResults.setProb(prob); return weightedLinearCorrResults; }
From source file:org.cirdles.calamari.algorithms.WeightedMeanCalculators.java
public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) { // assume varCov is variance-covariance matrix (i.e. SigRho = false) WtdAvCorrResults results = new WtdAvCorrResults(); int n = varCov.length; RealMatrix omegaInv = new BlockRealMatrix(varCov); RealMatrix omega = MatrixUtils.inverse(omegaInv); double numer = 0.0; double denom = 0.0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { numer += (values[i] + values[j]) * omega.getEntry(i, j); denom += omega.getEntry(i, j); }//from w ww .j a v a 2s . c o m } // test denom if (denom > 0.0) { double meanVal = numer / denom / 2.0; double meanValSigma = StrictMath.sqrt(1.0 / denom); double[][] unwtdResidsArray = new double[n][1]; for (int i = 0; i < n; i++) { unwtdResidsArray[i][0] = values[i] - meanVal; } RealMatrix unwtdResids = new BlockRealMatrix(unwtdResidsArray); RealMatrix transUnwtdResids = unwtdResids.transpose(); RealMatrix product = transUnwtdResids.multiply(omega); RealMatrix sumWtdResids = product.multiply(unwtdResids); double mswd = 0.0; double prob = 0.0; if (n > 1) { mswd = sumWtdResids.getEntry(0, 0) / (n - 1); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 1), 1E9); prob = 1.0 - fdist.cumulativeProbability(mswd); } results.setBad(false); results.setMeanVal(meanVal); results.setSigmaMeanVal(meanValSigma); results.setMswd(mswd); results.setProb(prob); } return results; }
From source file:org.cirdles.squid.algorithms.weightedMeans.WeightedMeanCalculators.java
/** * * @param y/* w w w . ja va 2s .co m*/ * @param x * @param sigmaRhoY * @return */ public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[] x, double[][] sigmaRhoY) { WeightedLinearCorrResults weightedLinearCorrResults = new WeightedLinearCorrResults(); RealMatrix omega = new BlockRealMatrix(convertCorrelationsToCovariances(sigmaRhoY)); RealMatrix invOmega = MatrixUtils.inverse(omega); int n = y.length; double mX = 0; double pX = 0; double pY = 0; double pXY = 0; double w = 0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { double invOm = invOmega.getEntry(i, j); w += invOm; pX += (invOm * (x[i] + x[j])); pY += (invOm * (y[i] + y[j])); pXY += (invOm * (((x[i] * y[j]) + (x[j] * y[i])))); mX += (invOm * x[i] * x[j]); } } double slope = ((2 * pXY * w) - (pX * pY)) / ((4 * mX * w) - (pX * pX)); double intercept = (pY - (slope * pX)) / (2 * w); RealMatrix fischer = new BlockRealMatrix(new double[][] { { mX, pX / 2.0 }, { pX / 2.0, w } }); RealMatrix fischerInv = MatrixUtils.inverse(fischer); double slopeSig = Math.sqrt(fischerInv.getEntry(0, 0)); double interceptSig = Math.sqrt(fischerInv.getEntry(1, 1)); double slopeInterceptCov = fischerInv.getEntry(0, 1); RealMatrix resid = new BlockRealMatrix(n, 1); for (int i = 0; i < n; i++) { resid.setEntry(i, 0, y[i] - (slope * x[i]) - intercept); } RealMatrix residT = resid.transpose(); RealMatrix mM = residT.multiply(invOmega).multiply(resid); double sumSqWtdResids = mM.getEntry(0, 0); double mswd = sumSqWtdResids / (n - 2); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 2), 1E9); double prob = 1.0 - fdist.cumulativeProbability(mswd); weightedLinearCorrResults.setBad(false); weightedLinearCorrResults.setSlope(slope); weightedLinearCorrResults.setIntercept(intercept); weightedLinearCorrResults.setSlopeSig(slopeSig); weightedLinearCorrResults.setInterceptSig(interceptSig); weightedLinearCorrResults.setSlopeInterceptCov(slopeInterceptCov); weightedLinearCorrResults.setMswd(mswd); weightedLinearCorrResults.setProb(prob); return weightedLinearCorrResults; }
From source file:org.cirdles.squid.algorithms.weightedMeans.WeightedMeanCalculators.java
/** * * @param values/*from ww w.j a v a2s. c o m*/ * @param varCov * @return */ public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) { // assume varCov is variance-covariance matrix (i.e. SigRho = false) WtdAvCorrResults results = new WtdAvCorrResults(); int n = varCov.length; RealMatrix omegaInv = new BlockRealMatrix(varCov); RealMatrix omega = MatrixUtils.inverse(omegaInv); double numer = 0.0; double denom = 0.0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { numer += (values[i] + values[j]) * omega.getEntry(i, j); denom += omega.getEntry(i, j); } } // test denom if (denom > 0.0) { double meanVal = numer / denom / 2.0; double meanValSigma = Math.sqrt(1.0 / denom); double[][] unwtdResidsArray = new double[n][1]; for (int i = 0; i < n; i++) { unwtdResidsArray[i][0] = values[i] - meanVal; } RealMatrix unwtdResids = new BlockRealMatrix(unwtdResidsArray); RealMatrix transUnwtdResids = unwtdResids.transpose(); RealMatrix product = transUnwtdResids.multiply(omega); RealMatrix sumWtdResids = product.multiply(unwtdResids); double mswd = 0.0; double prob = 0.0; if (n > 1) { mswd = sumWtdResids.getEntry(0, 0) / (n - 1); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 1), 1E9); prob = 1.0 - fdist.cumulativeProbability(mswd); } results.setBad(false); results.setMeanVal(meanVal); results.setSigmaMeanVal(meanValSigma); results.setMswd(mswd); results.setProb(prob); } return results; }
From source file:org.knime.al.util.noveltydetection.knfst.KNFST.java
public static RealMatrix projection(final RealMatrix kernelMatrix, final String[] labels) throws KNFSTException { final ClassWrapper[] classes = ClassWrapper.classes(labels); // check labels if (classes.length == 1) { throw new IllegalArgumentException( "not able to calculate a nullspace from data of a single class using KNFST (input variable \"labels\" only contains a single value)"); }/*from w w w .j a va2 s . c o m*/ // check kernel matrix if (!kernelMatrix.isSquare()) { throw new IllegalArgumentException("The KernelMatrix must be quadratic!"); } // calculate weights of orthonormal basis in kernel space final RealMatrix centeredK = centerKernelMatrix(kernelMatrix); final EigenDecomposition eig = new EigenDecomposition(centeredK); final double[] eigVals = eig.getRealEigenvalues(); final ArrayList<Integer> nonZeroEigValIndices = new ArrayList<Integer>(); for (int i = 0; i < eigVals.length; i++) { if (eigVals[i] > 1e-12) { nonZeroEigValIndices.add(i); } } int eigIterator = 0; final RealMatrix eigVecs = eig.getV(); RealMatrix basisvecs = null; try { basisvecs = MatrixUtils.createRealMatrix(eigVecs.getRowDimension(), nonZeroEigValIndices.size()); } catch (final Exception e) { throw new KNFSTException("Something went wrong. Try different parameters or a different kernel."); } for (final Integer index : nonZeroEigValIndices) { final double normalizer = 1 / Math.sqrt(eigVals[index]); final RealVector basisVec = eigVecs.getColumnVector(eigIterator).mapMultiply(normalizer); basisvecs.setColumnVector(eigIterator++, basisVec); } // calculate transformation T of within class scatter Sw: // T= B'*K*(I-L) and L a block matrix final RealMatrix L = kernelMatrix.createMatrix(kernelMatrix.getRowDimension(), kernelMatrix.getColumnDimension()); int start = 0; for (final ClassWrapper cl : classes) { final int count = cl.getCount(); L.setSubMatrix(MatrixFunctions.ones(count, count).scalarMultiply(1.0 / count).getData(), start, start); start += count; } // need Matrix M with all entries 1/m to modify basisvecs which allows // usage of // uncentered kernel values (eye(size(M)).M)*basisvecs final RealMatrix M = MatrixFunctions .ones(kernelMatrix.getColumnDimension(), kernelMatrix.getColumnDimension()) .scalarMultiply(1.0 / kernelMatrix.getColumnDimension()); final RealMatrix I = MatrixUtils.createRealIdentityMatrix(M.getColumnDimension()); // compute helper matrix H final RealMatrix H = ((I.subtract(M)).multiply(basisvecs)).transpose().multiply(kernelMatrix) .multiply(I.subtract(L)); // T = H*H' = B'*Sw*B with B=basisvecs final RealMatrix T = H.multiply(H.transpose()); // calculate weights for null space RealMatrix eigenvecs = MatrixFunctions.nullspace(T); if (eigenvecs == null) { final EigenDecomposition eigenComp = new EigenDecomposition(T); final double[] eigenvals = eigenComp.getRealEigenvalues(); eigenvecs = eigenComp.getV(); final int minId = MatrixFunctions.argmin(MatrixFunctions.abs(eigenvals)); final double[] eigenvecsData = eigenvecs.getColumn(minId); eigenvecs = MatrixUtils.createColumnRealMatrix(eigenvecsData); } // System.out.println("eigenvecs:"); // test.printMatrix(eigenvecs); // calculate null space projection final RealMatrix proj = ((I.subtract(M)).multiply(basisvecs)).multiply(eigenvecs); return proj; }
From source file:org.knime.al.util.noveltydetection.knfst.MultiClassKNFST.java
private NoveltyScores score(final RealMatrix kernelMatrix) { final RealMatrix projectionVectors = kernelMatrix.transpose().multiply(m_projection); // squared euclidean distances to target points: final RealMatrix squared_distances = squared_euclidean_distances(projectionVectors, m_targetPoints); // novelty scores as minimum distance to one of the target points final RealVector scoreVector = MatrixFunctions.sqrt(MatrixFunctions.rowMins(squared_distances)); return new NoveltyScores(scoreVector.toArray(), projectionVectors); }