List of usage examples for org.apache.commons.math3.linear RealMatrix getRow
double[] getRow(int row) throws OutOfRangeException;
From source file:se.toxbee.sleepfighter.challenge.math.LinearEquationProblem.java
private static RealVector multiply(RealMatrix matrix, RealVector vector) { double[] m1 = matrix.getRow(0); double[] m2 = matrix.getRow(1); double[] m3 = matrix.getRow(2); double[] v = vector.toArray(); return new ArrayRealVector(new double[] { m1[0] * v[0] + m1[1] * v[1] + m1[2] * v[2], m2[0] * v[0] + m2[1] * v[1] + m2[2] * v[2], m3[0] * v[0] + m3[1] * v[1] + m3[2] * v[2], });/*from w w w .j ava 2 s . com*/ }
From source file:tinfour.gwr.SurfaceGwr.java
public void computeVarianceAndHat() { if (areVarianceAndHatPrepped) { return;//from w ww .j ava2 s . co m } areVarianceAndHatPrepped = true; if (sampleWeightsMatrix == null) { throw new NullPointerException("Null specification for sampleWeightsMatrix"); } else if (sampleWeightsMatrix.length != nSamples) { throw new IllegalArgumentException("Incorrectly specified sampleWeightsMatrix"); } double[][] bigS = new double[nSamples][nSamples]; double[][] bigW = sampleWeightsMatrix; double[][] input = computeDesignMatrix(model, xOffset, yOffset, 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 sTrace = 0; double sTrace2 = 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 = cdSolver.getInverse(); RealMatrix r = rx.multiply(cInv).multiply(mXTW); double[] row = r.getRow(0); sTrace += row[i]; System.arraycopy(row, 0, bigS[i], 0, nSamples); for (int j = 0; j < nSamples; j++) { sTrace2 += row[j] * row[j]; } } hat = new BlockRealMatrix(bigS); traceHat = sTrace; traceHat2 = sTrace2; 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 = hat.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; } rss = sse; double d1 = nSamples - (2 * traceHat - sTrace2); sigma2 = rss / d1; mlSigma2 = rss / nSamples; RealMatrix mIL = hat.copy(); for (int i = 0; i < nSamples; i++) { double c = 1.0 - mIL.getEntry(i, i); mIL.setEntry(i, i, c); } RealMatrix mILT = mIL.transpose().multiply(mIL); delta1 = mILT.getTrace(); delta2 = (mILT.multiply(mILT)).getTrace(); }
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 */// w w w .j av a 2 s.co m 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); }