List of usage examples for org.apache.commons.math3.linear RealVector getEntry
public abstract double getEntry(int index) throws OutOfRangeException;
From source file:com.datumbox.framework.machinelearning.regression.MatrixLinearRegression.java
@Override protected void estimateModelParameters(Dataset trainingData) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = trainingData.size(); int d = trainingData.getColumnSize() + 1;//plus one for the constant //initialization modelParameters.setN(n);//from w w w.ja v a2s. co m modelParameters.setD(d); Map<Object, Double> thitas = modelParameters.getThitas(); Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); MatrixDataset matrixDataset = MatrixDataset.newInstance(trainingData, true, featureIds); RealVector Y = matrixDataset.getY(); RealMatrix X = matrixDataset.getX(); //(X'X)^-1 RealMatrix Xt = X.transpose(); LUDecomposition lud = new LUDecomposition(Xt.multiply(X)); RealMatrix XtXinv = lud.getSolver().getInverse(); lud = null; //(X'X)^-1 * X'Y RealVector coefficients = XtXinv.multiply(Xt).operate(Y); Xt = null; //put the features coefficients in the thita map thitas.put(Dataset.constantColumnName, coefficients.getEntry(0)); for (Map.Entry<Object, Integer> entry : featureIds.entrySet()) { Object feature = entry.getKey(); Integer featureId = entry.getValue(); thitas.put(feature, coefficients.getEntry(featureId)); } if (knowledgeBase.getTrainingParameters().getCalculatePvalue()) { //calculate them only if //get the predictions and subtact the Y vector. Sum the squared differences to get the error double SSE = 0.0; for (double v : X.operate(coefficients).subtract(Y).toArray()) { SSE += v * v; } Y = null; //standard error matrix double MSE = SSE / (n - d); //mean square error = SSE / dfResidual RealMatrix SE = XtXinv.scalarMultiply(MSE); XtXinv = null; //creating a flipped map of ids to features Map<Integer, Object> idsFeatures = PHPfunctions.array_flip(featureIds); Map<Object, Double> pvalues = new HashMap<>(); //This is not small, but it does not make sense to store it in the db for (int i = 0; i < d; ++i) { double error = SE.getEntry(i, i); Object feature = idsFeatures.get(i); if (error <= 0.0) { //double tstat = Double.MAX_VALUE; pvalues.put(feature, 0.0); } else { double tstat = coefficients.getEntry(i) / Math.sqrt(error); pvalues.put(feature, 1.0 - ContinuousDistributions.StudentsCdf(tstat, n - d)); //n-d degrees of freedom } } SE = null; coefficients = null; idsFeatures = null; matrixDataset = null; modelParameters.setFeaturePvalues(pvalues); } }
From source file:edu.oregonstate.eecs.mcplan.ml.KulisLowRankKernelLearner.java
private void choleskyUpdate(final int r, final double alpha, final RealVector x, final RealMatrix B) { final double[] a = new double[r + 1]; final double[] h = new double[r]; a[0] = alpha;//from w w w . j av a 2 s . co m for (int i = 0; i < r; ++i) { final double xi = x.getEntry(i); double t = 1 + a[i] * xi * xi; h[i] = Math.sqrt(t); a[i + 1] = a[i] / t; t = B.getEntry(i, i); double s = 0.0; B.setEntry(i, i, t * h[i]); // t == B_ii for (int j = i - 1; j >= 0; --j) { s += t * x.getEntry(j + 1); t = B.getEntry(i, j); B.setEntry(i, j, (t + a[j + 1] * x.getEntry(j) * s) * h[j]); // t == B_ij } } }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the standardized adjusted residuals (according to the same definition used by MATLAB) of the data points for fitting. * * @param indVarValues The values of the independent variable used for the fitting. * @param depVarValues The values of the dependent variable used for the fitting. * @param leverages the leverages of the independent variables, as compted by {@link #calculateLeverages(RealVector)} * @param fitParams the results of a (possibly weighted) least squares fit to the data, containing one or two components: a slope and an optional y-intercept. * @return a RealVector containing an adjusted residual value for each data point */// w w w . jav a 2s. c o m protected RealVector calculateStandardizedAdjustedResiduals(RealVector indVarValues, RealVector depVarValues, RealVector leverages, RealVector fitParams) { RealVector predictedValues = indVarValues.mapMultiply(fitParams.getEntry(0)); RealVector denom = leverages.mapMultiply(-1.0).mapAddToSelf(1 + this.CLOSE_TO_ZERO) .mapToSelf(new org.apache.commons.math3.analysis.function.Sqrt()); if (!this.noIntercept) { predictedValues = predictedValues.mapAdd(fitParams.getEntry(1)); } double stddev = 0; double mean = 0; for (int i = 0; i < depVarValues.getDimension(); i++) { mean += depVarValues.getEntry(i); } mean /= depVarValues.getDimension(); stddev = depVarValues.mapSubtract(mean).getNorm() * (depVarValues.getDimension() * 1.0 / (depVarValues.getDimension() - 1)); RealVector residuals = depVarValues.subtract(predictedValues).ebeDivide(denom); RealVector absDev = residuals.map(new org.apache.commons.math3.analysis.function.Abs()); int smallerDim = 2; if (this.noIntercept) { smallerDim = 1; } double[] resArray = residuals.map(new org.apache.commons.math3.analysis.function.Abs()).toArray(); java.util.Arrays.sort(resArray); RealVector partialRes = new ArrayRealVector(absDev.getDimension() - smallerDim + 1, 0.0); for (int i = smallerDim - 1; i < resArray.length; i++) { partialRes.setEntry(i - smallerDim + 1, resArray[i]); } double resMAD = 0; if (partialRes.getDimension() % 2 == 0) { resMAD = LocalBackgroundEstimationFilter.quickFindKth(partialRes.getDimension() / 2, partialRes) + LocalBackgroundEstimationFilter.quickFindKth(partialRes.getDimension() / 2 - 1, partialRes); resMAD /= 2.0; } else { resMAD = LocalBackgroundEstimationFilter.quickFindKth((partialRes.getDimension() - 1) / 2, partialRes); } resMAD /= 0.6745; if (resMAD < stddev * CLOSE_TO_ZERO) { resMAD = stddev * CLOSE_TO_ZERO; } return residuals.mapDivide(DEFAULT_TUNING_CONST * resMAD); }
From source file:ellipsoidFit.FitPoints.java
/** * Create a matrix in the algebraic form of the polynomial Ax^2 + By^2 + * Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1. * /*from w ww . j ava2 s . c o m*/ * @param v * the vector polynomial. * @return the matrix of the algebraic form of the polynomial. */ private RealMatrix formAlgebraicMatrix(RealVector v) { // a = // [ Ax^2 2Dxy 2Exz 2Gx ] // [ 2Dxy By^2 2Fyz 2Hy ] // [ 2Exz 2Fyz Cz^2 2Iz ] // [ 2Gx 2Hy 2Iz -1 ] ] RealMatrix a = new Array2DRowRealMatrix(4, 4); a.setEntry(0, 0, v.getEntry(0)); a.setEntry(0, 1, v.getEntry(3)); a.setEntry(0, 2, v.getEntry(4)); a.setEntry(0, 3, v.getEntry(6)); a.setEntry(1, 0, v.getEntry(3)); a.setEntry(1, 1, v.getEntry(1)); a.setEntry(1, 2, v.getEntry(5)); a.setEntry(1, 3, v.getEntry(7)); a.setEntry(2, 0, v.getEntry(4)); a.setEntry(2, 1, v.getEntry(5)); a.setEntry(2, 2, v.getEntry(2)); a.setEntry(2, 3, v.getEntry(8)); a.setEntry(3, 0, v.getEntry(6)); a.setEntry(3, 1, v.getEntry(7)); a.setEntry(3, 2, v.getEntry(8)); a.setEntry(3, 3, -1); return a; }
From source file:com.joptimizer.util.MPSParserNetlibTest.java
/** * Tests the parsing of a netlib problem. *//*from w ww. ja v a 2 s . c o m*/ public void xxxtestSingleNetlib() throws Exception { log.debug("testSingleNetlib"); //String problemName = "afiro"; //String problemName = "afiroPresolved"; //String problemName = "adlittle"; //String problemName = "kb2"; //String problemName = "sc50a"; //String problemName = "sc50b"; //String problemName = "blend"; //String problemName = "scorpion"; //String problemName = "recipe"; //String problemName = "recipePresolved"; //String problemName = "sctap1"; //String problemName = "fit1d"; //String problemName = "israel"; //String problemName = "grow15"; //String problemName = "etamacro"; //String problemName = "pilot"; //String problemName = "pilot4"; //String problemName = "osa-14"; //String problemName = "brandyPresolved"; String problemName = "maros"; File f = Utils.getClasspathResourceAsFile("lp" + File.separator + "netlib" + File.separator + problemName + File.separator + problemName + ".mps"); MPSParser mpsParser = new MPSParser(); mpsParser.parse(f); Properties expectedSolProps = null; try { //this is the solution of the mps problem given by Mathematica expectedSolProps = load(Utils.getClasspathResourceAsFile( "lp" + File.separator + "netlib" + File.separator + problemName + File.separator + "sol.txt")); } catch (Exception e) { } log.debug("name: " + mpsParser.getName()); log.debug("n : " + mpsParser.getN()); log.debug("meq : " + mpsParser.getMeq()); log.debug("mieq: " + mpsParser.getMieq()); log.debug("meq+mieq: " + (mpsParser.getMeq() + mpsParser.getMieq())); List<String> variablesNames = mpsParser.getVariablesNames(); log.debug("x: " + ArrayUtils.toString(variablesNames)); // log.debug("c: " + ArrayUtils.toString(p.getC())); // log.debug("G: " + ArrayUtils.toString(p.getG())); // log.debug("h: " + ArrayUtils.toString(p.getH())); // log.debug("A: " + ArrayUtils.toString(p.getA())); // log.debug("b: " + ArrayUtils.toString(p.getB())); // log.debug("lb:" + ArrayUtils.toString(p.getLb())); // log.debug("ub:" + ArrayUtils.toString(p.getUb())); //check consistency: if the problem was correctly parsed, the expectedSol must be its solution double delta = 1.e-7; if (expectedSolProps != null) { //key = variable name //value = sol value assertEquals(expectedSolProps.size(), variablesNames.size()); RealVector expectedSol = new ArrayRealVector(variablesNames.size()); for (int i = 0; i < variablesNames.size(); i++) { expectedSol.setEntry(i, Double.parseDouble(expectedSolProps.getProperty(variablesNames.get(i)))); } log.debug("expectedSol: " + ArrayUtils.toString(expectedSol.toArray())); //check objective function value Map<String, LPNetlibProblem> problemsMap = LPNetlibProblem.loadAllProblems(); LPNetlibProblem problem = problemsMap.get(problemName); RealVector c = new ArrayRealVector(mpsParser.getC().toArray()); double value = c.dotProduct(expectedSol); log.debug("optimalValue: " + problem.optimalValue); log.debug("value : " + value); assertEquals(problem.optimalValue, value, delta); //check G.x < h if (mpsParser.getG() != null) { RealMatrix G = new Array2DRowRealMatrix(mpsParser.getG().toArray()); RealVector h = new ArrayRealVector(mpsParser.getH().toArray()); RealVector Gxh = G.operate(expectedSol).subtract(h); double maxGxh = -Double.MAX_VALUE; for (int i = 0; i < Gxh.getDimension(); i++) { //log.debug(i); maxGxh = Math.max(maxGxh, Gxh.getEntry(i)); assertTrue(Gxh.getEntry(i) <= 0); } log.debug("max(G.x - h): " + maxGxh); } //check A.x = b if (mpsParser.getA() != null) { RealMatrix A = new Array2DRowRealMatrix(mpsParser.getA().toArray()); RealVector b = new ArrayRealVector(mpsParser.getB().toArray()); RealVector Axb = A.operate(expectedSol).subtract(b); double norm = Axb.getNorm(); log.debug("||A.x -b||: " + norm); assertEquals(0., norm, delta * mpsParser.getN());//some more tolerance } //check upper and lower bounds for (int i = 0; i < mpsParser.getLb().size(); i++) { double di = Double.isNaN(mpsParser.getLb().getQuick(i)) ? -Double.MAX_VALUE : mpsParser.getLb().getQuick(i); assertTrue(di <= expectedSol.getEntry(i)); } for (int i = 0; i < mpsParser.getUb().size(); i++) { double di = Double.isNaN(mpsParser.getUb().getQuick(i)) ? Double.MAX_VALUE : mpsParser.getUb().getQuick(i); assertTrue(di >= expectedSol.getEntry(i)); } } Utils.writeDoubleArrayToFile(mpsParser.getC().toArray(), "target" + File.separator + "c.txt"); Utils.writeDoubleMatrixToFile(mpsParser.getG().toArray(), "target" + File.separator + "G.csv"); Utils.writeDoubleArrayToFile(mpsParser.getH().toArray(), "target" + File.separator + "h.txt"); Utils.writeDoubleMatrixToFile(mpsParser.getA().toArray(), "target" + File.separator + "A.csv"); Utils.writeDoubleArrayToFile(mpsParser.getB().toArray(), "target" + File.separator + "b.txt"); Utils.writeDoubleArrayToFile(mpsParser.getLb().toArray(), "target" + File.separator + "lb.txt"); Utils.writeDoubleArrayToFile(mpsParser.getUb().toArray(), "target" + File.separator + "ub.txt"); }
From source file:edu.washington.gs.skyline.model.quantification.LinearModel.java
public List<LinearFitResult> fit(double[] observations) { if (observations.length != designMatrix.getRowDimension()) { throw new IllegalArgumentException("Wrong number of rows"); }/*from ww w . ja v a 2 s .c om*/ RealVector coefficients = qrDecomposition.getSolver().solve(new ArrayRealVector(observations)); RealVector fittedValues = new ArrayRealVector(observations.length); RealVector residuals = new ArrayRealVector(observations.length); for (int iRow = 0; iRow < observations.length; iRow++) { RealVector designRow = designMatrix.getRowVector(iRow); fittedValues.setEntry(iRow, designRow.dotProduct(coefficients)); residuals.setEntry(iRow, observations[iRow] - fittedValues.getEntry(iRow)); } double rss = residuals.dotProduct(residuals); int degreesOfFreedom = observations.length - qrDecomposition.getR().getColumnDimension(); double resVar = rss / degreesOfFreedom; double sigma = Math.sqrt(resVar); RealMatrix covarianceUnscaled = matrixCrossproductInverse; RealMatrix scaledCovariance = covarianceUnscaled.scalarMultiply(sigma * sigma); List<LinearFitResult> results = new ArrayList<>(); for (int iContrastRow = 0; iContrastRow < contrastValues.getRowDimension(); iContrastRow++) { RealVector contrastRow = contrastValues.getRowVector(iContrastRow); double standardError = 0; for (int iRow = 0; iRow < independentColumnIndices.length; iRow++) { for (int iCol = 0; iCol < independentColumnIndices.length; iCol++) { standardError = contrastRow.getEntry(independentColumnIndices[iRow]) * scaledCovariance.getEntry(iRow, iCol) * contrastRow.getEntry(independentColumnIndices[iCol]); } } standardError = Math.sqrt(standardError); double foldChange = coefficients.dotProduct(contrastRow); LinearFitResult linearFitResult = new LinearFitResult(foldChange); double tValue = foldChange / standardError; linearFitResult.setTValue(tValue); linearFitResult.setStandardError(standardError); linearFitResult.setDegreesOfFreedom(degreesOfFreedom); if (0 == degreesOfFreedom) { linearFitResult.setPValue(1.0); } else { TDistribution tDistribution = new TDistribution(degreesOfFreedom); double pValue = (1 - tDistribution.cumulativeProbability(Math.abs(tValue))) * 2; linearFitResult.setPValue(pValue); } results.add(linearFitResult); } return results; }
From source file:edu.oregonstate.eecs.mcplan.ml.InformationTheoreticMetricLearner.java
/** * Execute the algorithm./* w w w. j a v a 2 s . co m*/ * @see java.lang.Runnable#run() */ @Override public void run() { System.out.println("\tITML run()"); final RealMatrix A = A0_.copy(); final int[] idx = Fn.range(0, Nc_); int iter_count = 0; final int logging_interval = 1; double cumulative_update = 0.0; while (iter_count++ < iteration_limit_) { if (iter_count % logging_interval == 0) { System.out.println("\tIteration " + iter_count); } Fn.shuffle(rng_, idx); double update_magnitude = 0.0; for (int c = 0; c < Nc_; ++c) { final int[] constraint; final RealVector xdiff; final int delta; if (c < S_.size()) { // constraint = S_.get( c ); xdiff = new ArrayRealVector(S_.get(c)); delta = 1; } else { // constraint = D_.get( c - S_.size() ); xdiff = new ArrayRealVector(D_.get(c - S_.size())); delta = -1; } // final int i = constraint[0]; // final int j = constraint[1]; // final RealVector xdiff = X_.get( i ).subtract( X_.get( j ) ); // final double p = xdiff.dotProduct( A.operate( xdiff ) ); // if( p == 0.0 ) { // System.out.println( "! p == 0.0" ); // } // if( xi_[c] == 0.0 ) { // System.out.println( "! xi_[c] == 0.0" ); // } final double p = HilbertSpace.inner_prod(xdiff, A, xdiff); final double alpha; if (p == 0.0) { alpha = lambda_[c]; } else { alpha = Math.min(lambda_[c], (delta / 2.0) * ((1.0 / p) - (gamma_ / xi_[c]))); } final double beta = (delta * alpha) / (1 - delta * alpha * p); xi_[c] = (gamma_ * xi_[c]) / (gamma_ + delta * alpha * xi_[c]); lambda_[c] -= alpha; // This implements: A += beta * A xdiff xdiff^T A final RealVector Ax = A.operate(xdiff); // final RealMatrix outer_prod = Ax.outerProduct( Ax ); for (int row = 0; row < A.getRowDimension(); ++row) { final double axi = Ax.getEntry(row); for (int col = 0; col < A.getColumnDimension(); ++col) { final double a = axi * Ax.getEntry(col); A.addToEntry(row, col, a * beta); } } // final RealMatrix outer_prod = xdiff.outerProduct( xdiff ); // final RealMatrix update = A.multiply( outer_prod ).multiply( A ).scalarMultiply( beta ); // A = A.add( update ); update_magnitude += alpha * alpha; //update.getFrobeniusNorm(); } cumulative_update += update_magnitude; if (iter_count % logging_interval == 0) { System.out.println("\tupdate = " + (cumulative_update / logging_interval)); cumulative_update = 0.0; } // Declare convergence if all updates were small. if (Math.sqrt(update_magnitude) < convergence_tolerance_) { break; } } A_ = A; // A_ = MatrixAlgorithms.makePositiveDefinite( A, 1e-4 ); // System.out.println( "Metric:" ); // for( int i = 0; i < A_.getRowDimension(); ++i ) { // System.out.println( A_.getRowVector( i ) ); // } // Check for positive-definiteness // final EigenDecomposition eig = new EigenDecomposition( A_ ); // final double det = eig.getDeterminant(); // assert( det > 0.0 ); }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the weight for the next weighted least squares iteration using the bisquare weighting function. * @param stdAdjR the standardized adjusted residuals, as computed by {@link #calculateStandardizedAdjustedResiduals(RealVector, RealVector, RealVector, RealVector)} * @return a RealVector containing weights for each data point suitable for weighted least squares fitting. *//*from w w w .j a v a2s . co m*/ protected RealVector calculateBisquareWeights(RealVector stdAdjR) { RealVector bisquareWeights = new ArrayRealVector(stdAdjR.getDimension(), 0.0); for (int i = 0; i < bisquareWeights.getDimension(); i++) { if (Math.abs(stdAdjR.getEntry(i)) < 1) { bisquareWeights.setEntry(i, Math.pow(1 - Math.pow(stdAdjR.getEntry(i), 2), 2)); } } return bisquareWeights; }
From source file:GeMSE.GS.Analysis.Stats.OneSamplePCAPanel.java
private void computePrincipalComponents() { RealMatrix realMatrix = MatrixUtils.createRealMatrix(_data); Covariance covariance = new Covariance(realMatrix); _covariance = covariance.getCovarianceMatrix(); EigenDecomposition ed = new EigenDecomposition(_covariance); double[] realEigenvalues = ed.getRealEigenvalues(); int pcaCols = numPCAIndices(realEigenvalues, _level); int eigenCount = realEigenvalues.length; _principalComponents = new Array2DRowRealMatrix(eigenCount, pcaCols); _variance = new ArrayRealVector(pcaCols); for (int i = 0; i < pcaCols; i++) { RealVector eigenVec = ed.getEigenvector(i); for (int j = 0; j < eigenCount; j++) _principalComponents.setEntry(j, i, eigenVec.getEntry(j)); _variance.setEntry(i, realEigenvalues[i]); }//from ww w. j a v a2s . c o m }
From source file:edu.oregonstate.eecs.mcplan.ml.ConstrainedKMeans.java
private boolean mstep() { final double[][] mu_prime = new double[k][d]; final int[] nh = new int[k]; for (int i = 0; i < N; ++i) { final int h = assignments_[i]; nh[h] += 1;/*from w ww. j a v a 2 s . c o m*/ final RealVector xi = X_.get(i); // mu_prime[h] = mu_prime[h].add( xi.subtract( mu_prime[h] ).mapDivideToSelf( i + 1 ) ); for (int j = 0; j < d; ++j) { // Incremental mean calculation mu_prime[h][j] += (xi.getEntry(j) - mu_prime[h][j]) / nh[h]; } } double Jprime = 0.0; for (int i = 0; i < N; ++i) { final int h = assignments_[i]; Jprime += J(i, h); } double delta = 0.0; for (int h = 0; h < k; ++h) { final RealVector v = new ArrayRealVector(mu_prime[h], false); // No copy // System.out.println( mu_[h] ); // System.out.println( v ); final RealVector diff = mu_[h].subtract(v); delta += diff.getNorm(); // delta += mu_[h].getDistance( v ); // delta += distance( mu_[h], v ); mu_[h] = v; } // System.out.println( "\tdelta = " + delta ); final double deltaJ = J_ - Jprime; // System.out.println( "\tdeltaJ = " + deltaJ ); J_ = Jprime; return deltaJ > convergence_tolerance_; }