List of usage examples for org.apache.commons.math3.linear RealMatrix getEntry
double getEntry(int row, int column) throws OutOfRangeException;
From source file:tinfour.gwr.SurfaceGwr.java
/** * Gets a value equal to one half of the range of the prediction interval * on the observed response at the interpolation coordinates for the * most recent call to computeRegression(). * * @param alpha the significance level (typically 0..05, etc). * @return a positive value./*from w w w . jav a2 s.c om*/ */ public double getPredictionIntervalHalfRange(double alpha) { // TO DO: if the method is OLS, it would make sense to // use a OLS version of this calculation rather than // the more costly Leung version... Also, I am not 100 % // sure that they converge to the same answer, though they should computeVarianceAndHat(); //double effDegOfF = getEffectiveDegreesOfFreedom(); // should match delta1 double[][] input = computeDesignMatrix(model, xOffset, yOffset, nSamples, samples); RealMatrix mX = new BlockRealMatrix(input); RealMatrix mXT = mX.transpose(); // the weights array may not necessarily be of dimension nSamples, // so we need to copy it double[] rW = Arrays.copyOf(weights, nSamples); RealMatrix mW = new DiagonalMatrix(rW); RealMatrix design = mXT.multiply(mW).multiply(mX); RealMatrix vcm; try { QRDecomposition cd = new QRDecomposition(design); DecompositionSolver s = cd.getSolver(); vcm = s.getInverse(); } catch (SingularMatrixException npex) { return Double.NaN; } double nLeungDOF = (delta1 * delta1 / delta2); for (int i = 0; i < nSamples; i++) { rW[i] = weights[i] * weights[i]; } DiagonalMatrix mW2 = new DiagonalMatrix(rW); RealMatrix mS = vcm.multiply(mXT).multiply(mW2).multiply(mX).multiply(vcm); double pS = mS.getEntry(0, 0); double p = Math.sqrt(this.sigma2 * (1 + pS)); TDistribution td = new TDistribution(nLeungDOF); double ta = td.inverseCumulativeProbability(1.0 - alpha / 2.0); return ta * p; }
From source file:tinfour.gwr.SurfaceGwr.java
/** * Evaluates the AICc score for the specified coordinates. This method * does not change the state of any of the member elements of this class. * It is intended to be used in the automatic bandwidth selection * operations implemented by calling classes. * * @param xQuery the x coordinate of the query point for evaluation * @param yQuery the y coordinate of the query point for evaluation * @param nSamples the number of samples. * @param samples an array of nSamples samples including x, y, and z. * @param weights an array of nSamples weights for each sample * @param lambda the bandwidth parameter for evaluation * @return if successful, a valid AICc; if unsuccessful, a NaN *//*from w ww .j a v a2s. c om*/ double evaluateAICc(SurfaceModel sm, double xQuery, double yQuery, int nSamples, double[][] samples, double[] weights, double[][] sampleWeightsMatrix) { // RealMatrix xwx = computeXWX(xQuery, yQuery, nSamples, samples, weights); double[][] bigS = new double[nSamples][nSamples]; double[][] bigW = sampleWeightsMatrix; double[][] input = computeDesignMatrix(sm, xQuery, yQuery, nSamples, samples); RealMatrix mX = new BlockRealMatrix(input); RealMatrix mXT = mX.transpose(); // in the loop below, we compute // Tr(hat) and Tr(Hat' x Hat) // this second term is actually the square of the // Frobenius Norm. Internally, the Apache Commons classes // may provide a more numerically stable implementation of this operation. // This may be worth future investigation. double traceS = 0; for (int i = 0; i < nSamples; i++) { DiagonalMatrix mW = new DiagonalMatrix(bigW[i]); //NOPMD RealMatrix mXTW = mXT.multiply(mW); RealMatrix rx = mX.getRowMatrix(i); RealMatrix c = mXTW.multiply(mX); QRDecomposition cd = new QRDecomposition(c); // NOPMD DecompositionSolver cdSolver = cd.getSolver(); RealMatrix cInv; try { cInv = cdSolver.getInverse(); } catch (SingularMatrixException | NullPointerException badMatrix) { return Double.NaN; } RealMatrix r = rx.multiply(cInv).multiply(mXTW); double[] row = r.getRow(0); traceS += row[i]; System.arraycopy(row, 0, bigS[i], 0, nSamples); //NOPMD } RealMatrix mS = new BlockRealMatrix(bigS); // the Hat matrix double[][] zArray = new double[nSamples][1]; for (int i = 0; i < nSamples; i++) { zArray[i][0] = samples[i][2]; } RealMatrix mY = new BlockRealMatrix(zArray); RealMatrix mYH = mS.multiply(mY); double sse = 0; for (int i = 0; i < nSamples; i++) { double yHat = mYH.getEntry(i, 0); double e = zArray[i][0] - yHat; sse += e * e; } double mls2 = sse / nSamples; double lv = Math.log(mls2); // this is 2*log(sqrt(mls2)) double x = (nSamples + traceS) / (nSamples - 2 - traceS); return nSamples * (lv + log2PI + x); }
From source file:UI.MainStageController.java
/** * creates the data for the analysis pane and displays it * data created includes:/* w w w . j a va 2 s . c om*/ * highest positive/negative correlations * average counts */ private void displayAnalysisTextsAndGraphs() { //Display node with highest frequency double highestFrequency = AnalysisData.getHighestFrequency(); TaxonNode nodeWithHighestFrequency = AnalysisData.getNodeWithHighestFrequency(); dataStatText.setText("Highest Frequency:\n" + nodeWithHighestFrequency.getName() + " (" + String.format("%.3f", highestFrequency) + ")\n"); //Display nodes with highest positive/negative correlation RealMatrix correlationMatrix = AnalysisData.getCorrelationMatrix(); int[] highestPositiveCorrelationCoordinates = AnalysisData.getHighestPositiveCorrelationCoordinates(); int[] highestNegativeCorrelationCoordinates = AnalysisData.getHighestNegativeCorrelationCoordinates(); LinkedList<TaxonNode> taxonList = SampleComparison.getUnifiedTaxonList(LoadedData.getSamplesToAnalyze(), AnalysisData.getLevelOfAnalysis()); TaxonNode hPCNode1 = taxonList.get(highestPositiveCorrelationCoordinates[0]); TaxonNode hPCNode2 = taxonList.get(highestPositiveCorrelationCoordinates[1]); TaxonNode hNCNode1 = taxonList.get(highestNegativeCorrelationCoordinates[0]); TaxonNode hNCNode2 = taxonList.get(highestNegativeCorrelationCoordinates[1]); dataStatText.setText(dataStatText.getText() + "\nHighest Positive Correlation:\n" + hPCNode1.getName() + " - " + hPCNode2.getName() + " (" + String.format("%.3f", correlationMatrix.getEntry(highestPositiveCorrelationCoordinates[0], highestPositiveCorrelationCoordinates[1])) + ")\n"); dataStatText.setText(dataStatText.getText() + "\nHighest Negative Correlation:\n" + hNCNode1.getName() + " - " + hNCNode2.getName() + " (" + String.format("%.3f", correlationMatrix.getEntry(highestNegativeCorrelationCoordinates[0], highestNegativeCorrelationCoordinates[1])) + ")"); //Generate Data for the pie chart frequencyChart.getData().clear(); HashMap<TaxonNode, Double> averageCounts = SampleComparison .calcAverageCounts(LoadedData.getSamplesToAnalyze(), AnalysisData.getLevelOfAnalysis()); for (TaxonNode taxonNode : averageCounts.keySet()) { PieChart.Data data = new PieChart.Data(taxonNode.getName(), averageCounts.get(taxonNode)); frequencyChart.getData().add(data); } analysisPane.setVisible(true); }
From source file:us.levk.math.linear.EucledianDistanceClusterer.java
public Cluster eucledian(final RealMatrix original) throws IOException { try (HugeRealMatrix distances = new HugeRealMatrix(original.getRowDimension(), original.getRowDimension())) { final Map<Integer, Cluster> genehash = new HashMap<Integer, Cluster>() { private static final long serialVersionUID = 1L; {/*from w w w. ja v a2 s. co m*/ for (int index = original.getRowDimension(); --index >= 0; put(index, new Cluster(index))) ; } }; TreeMap<Double, int[]> sorted = new TreeMap<>(); log.debug("Populating distance matrix"); for (int i = 0; i < original.getRowDimension(); i++) { for (int j = i + 1; j < original.getRowDimension(); j++) { // Euclidean distance calculation. double total = 0; for (int k = 0; k < original.getColumnDimension(); k++) { double left = original.getEntry(i, k); double right = original.getEntry(j, k); if (!isNaN(left) && !isNaN(right) && !isInfinite(left) && !isInfinite(right)) total += Math.pow(left - right, 2); } double distance = Math.pow(total, 0.5); distances.setEntry(i, j, distance); distances.setEntry(j, i, distance); int[] genePair = { i, j }; // Enter the distance calculated and the genes measured into a // treemap. Will be automatically sorted. sorted.put(distance, genePair); } } log.debug("Initialized distances matrix " + distances); while (true) { // Get the first key of the TreeMap. Will be the shortest distance de // facto. final double minkey = (Double) sorted.firstKey(); int[] minValues = (int[]) sorted.firstEntry().getValue(); final int value1 = minValues[0], value2 = minValues[1]; // find Cluster cluster = new Cluster(genehash.get(value1), genehash.get(value2)) { { log.debug("Generating cluster from " + value1 + " and " + value2 + " in " + genehash); contains().addAll(genehash.get(value1).contains()); contains().addAll(genehash.get(value2).contains()); d(minkey); log.debug("Generated cluster " + this); } }; genehash.put(cluster.id(), cluster); genehash.remove(value1); genehash.remove(value2); if (genehash.size() <= 1) break; // Iterate over all the current clusters to remeasure distance with the // previously clustered group. for (Cluster c : genehash.values()) { // Skip measuring the new cluster with itself. if (c == cluster) continue; double distance = 0; int n = 0; // Get genes from each cluster. Distance is measured from each element // to every element. for (int current : c.contains()) for (int created : cluster.contains()) { distance += distances.getEntry(current, created); n++; } distance = distance / n; int[] valuePair = { c.id(), cluster.id() }; sorted.put(distance, valuePair); } // Get the shortest distance. // Check to make sure shortest distance does not include a gene pair // that // has already had its elements clustered. boolean minimized = false; while (!minimized) { double mk = sorted.firstKey(); minValues = sorted.firstEntry().getValue(); // If the gene pair is not present in the current gene set, remove // this // distance. if (!genehash.containsKey(minValues[0]) || !genehash.containsKey(minValues[1])) sorted.remove(mk); else minimized = true; } } return genehash.entrySet().iterator().next().getValue(); } }