List of usage examples for org.apache.commons.math3.linear BlockRealMatrix BlockRealMatrix
public BlockRealMatrix(final double[][] rawData) throws DimensionMismatchException, NotStrictlyPositiveException
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 . j a va2s .co m */ 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 *//* ww w. j a va 2s.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:USVProsjekt.ThrustAllocator.java
private void setUp(double Lx1, double Lx2, double Ly1, double Ly2) { Phi = new double[][] { { 1., 0., 0., 0., 0., 0., 0. }, { 0., 1., 0., 0., 0., 0., 0., }, { 0., 0., 1., 0., 0., 0., 0. }, { 0., 0., 0., 1., 0., 0., 0. }, { 0., 0., 0., 0., 10., 0., 0. }, { 0., 0., 0., 0., 0., 10., 0. }, { 0., 0., 0., 0., 0., 0., 10. } }; A1 = new double[][] { { 1., 1., 0., 0., -1., 0., 0. }, { 0., 0., -1., -1., 0., -1., 0. }, { -Ly1, -Ly2, -Lx1, -Lx2, 0., 0., -1. } }; A2 = new BlockRealMatrix(new double[][] { { -1., 0., 0., 0., 0., 0., 0. }, { 0., -1., 0., 0., 0., 0., 0. }, { 0., 0., -1., 0., 0., 0., 0. }, { 0., 0., 0., -1., 0., 0., 0. }, { 1., 0., 0., 0., 0., 0., 0. }, { 0., 1., 0., 0., 0., 0., 0. }, { 0., 0., 1., 0., 0., 0., 0. }, { 0., 0., 0., 1., 0., 0., 0. } }); C2 = new BlockRealMatrix(new double[][] { { 0., 0., 0., -1., 0., 0., 0., 0., 0., 0., 0. }, { 0., 0., 0., 0., -1., 0., 0., 0., 0., 0., 0. }, { 0., 0., 0., 0., 0., -1., 0., 0., 0., 0., 0. }, { 0., 0., 0., 0., 0., 0., -1., 0., 0., 0., 0. }, { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0. }, { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0. }, { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0. }, { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1. } }); p = new ArrayRealVector(new double[] { 0., 0., 0., -29., -29., -29., -29., 34., 34., 34., 34. }); objectiveFunction = new PDQuadraticMultivariateRealFunction(Phi, null, 0); inequalities = new ConvexMultivariateRealFunction[8]; opt = new JOptimizer(); }