Example usage for org.apache.commons.math3.linear RealVector subtract

List of usage examples for org.apache.commons.math3.linear RealVector subtract

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.linear.*
  5. RealVector
  6. subtract

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear RealVector subtract.

Prototype

public RealVector subtract(RealVector v) throws DimensionMismatchException

Source Link

Document

Subtract v from this vector.

Usage

From source file:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java

@Override
public void run() {
    init();/*from   w w w.j a v a 2 s.c  om*/
    System.out.println("Init");
    for (int i = 0; i < mu().length; ++i) {
        System.out.println("Mu " + i + ": " + mu()[i]);
        System.out.println("Sigma " + i + ": " + Sigma()[i]);
    }

    int iterations = 0;
    while (!converged_ && iterations++ < max_iterations_) {
        // Expectation
        makeMixture();
        for (int i = 0; i < n_; ++i) {
            for (int j = 0; j < k_; ++j) {
                c_[i][j] = posterior(data_[i], j);
            }
            Fn.normalize_inplace(c_[i]);
        }

        // Maximization
        for (int j = 0; j < k_; ++j) {
            double Z = 0.0;
            final RealVector mu_j = new ArrayRealVector(d_);
            RealMatrix Sigma_j = new Array2DRowRealMatrix(d_, d_);
            for (int i = 0; i < n_; ++i) {
                final double c_ij = c_[i][j];
                Z += c_ij;
                final RealVector x_i = new ArrayRealVector(data_[i]);
                // mu_j += c_ij * x_i
                mu_j.combineToSelf(1.0, 1.0, x_i.mapMultiply(c_ij));
                final RealVector v = x_i.subtract(mu_[j]);
                // Sigma_j += c_ij * |v><v|
                Sigma_j = Sigma_j.add(v.outerProduct(v).scalarMultiply(c_ij));
            }
            final double Zinv = 1.0 / Z;
            final double pi_j = Z / n_;
            mu_j.mapMultiplyToSelf(Zinv);
            Sigma_j = Sigma_j.scalarMultiply(Zinv);
            //            converged &= hasConverged( j, pi_j, mu_j, Sigma_j );
            pi_[j] = pi_j;
            mu_[j] = mu_j;
            Sigma_[j] = Sigma_j;
        }
        //         debug();

        final double log_likelihood = logLikelihood();
        if (Math.abs(log_likelihood - old_log_likelihood_) < epsilon_) {
            converged_ = true;
        }
        old_log_likelihood_ = log_likelihood;
    }
}

From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java

@Override
public GaussianFactor division(GaussianFactor other) {
    if (!scope.contains(other.getVariables().getVariableIds())) {
        throw new FactorOperationException("Divisor scope " + other.getVariables() + " is not a subset of"
                + " this factor's scope " + scope);
    }/* w ww  . j ava  2s  .  c  o  m*/

    int[] otherMapping = scope.createContinuousVariableMapping(other.getVariables());

    RealMatrix newPrecisionMatrix = precisionMatrix.copy();

    RealMatrix otherPrecisionMatrix = other.getPrecisionMatrix();

    for (int i = 0; i < scope.size(); i++) {
        RealVector column = newPrecisionMatrix.getColumnVector(i);

        if (otherMapping[i] >= 0) {
            column = column.subtract(padVector(otherPrecisionMatrix.getColumnVector(otherMapping[i]),
                    scope.size(), otherMapping));
            newPrecisionMatrix.setColumnVector(i, column);
        }
    }

    RealVector newScaledMeanVector = scaledMeanVector.copy();
    RealVector otherScaledMeanVector = other.getScaledMeanVector();
    newScaledMeanVector = newScaledMeanVector
            .subtract(padVector(otherScaledMeanVector, scope.size(), otherMapping));

    double newNormalizationConstant = normalizationConstant - other.getNormalizationConstant();

    return new CanonicalGaussianFactor(scope, newPrecisionMatrix, newScaledMeanVector,
            newNormalizationConstant);
}

From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java

@Override
public GaussianFactor marginal(Scope marginalizationScope) {
    // the following assumes that the precision matrix (and mean vector) can be restructured as follows:
    // ( SUBMATRIX_XX SUBMATRIX_XY )
    // ( SUBMATRIX_YX SUBMATRIX_YY )
    // where X indicates the entries for the variables that are kept (i.e. in the scope argument) and
    // Y the variables that are marginalized out

    if (marginalizationScope.contains(scope)) {
        return this;
    }//from  ww w  . j av  a  2  s .  c  o m

    Scope newScope = scope.intersect(marginalizationScope);
    Scope scopeToMarginalize = scope.reduceBy(newScope);

    int[] xMapping = newScope.createContinuousVariableMapping(scope);
    RealMatrix xxMatrix = precisionMatrix.getSubMatrix(xMapping, xMapping);

    int[] yMapping = scopeToMarginalize.createContinuousVariableMapping(scope);
    RealMatrix yyMatrix = precisionMatrix.getSubMatrix(yMapping, yMapping);

    RealMatrix xyMatrix = precisionMatrix.getSubMatrix(xMapping, yMapping);
    RealMatrix yxMatrix = precisionMatrix.getSubMatrix(yMapping, xMapping);

    MathUtil yyUtil = new MathUtil(yyMatrix);
    RealMatrix yyInverse = yyUtil.invert();
    RealMatrix newPrecisionMatrix = xxMatrix.subtract(xyMatrix.multiply(yyInverse.multiply(yxMatrix)));

    RealVector xVector = getSubVector(scaledMeanVector, xMapping);
    RealVector yVector = getSubVector(scaledMeanVector, yMapping);

    RealVector newScaledMeanVector = xVector.subtract(xyMatrix.operate(yyInverse.operate(yVector)));

    MathUtil inverseUtil = new MathUtil(yyInverse.scalarMultiply(2.0d * Math.PI));
    double newNormalizationConstant = normalizationConstant
            + 0.5d * (Math.log(inverseUtil.determinant()) + yVector.dotProduct(yyInverse.operate(yVector)));

    return new CanonicalGaussianFactor(newScope, newPrecisionMatrix, newScaledMeanVector,
            newNormalizationConstant);
}

From source file:hivemall.utils.math.StatsUtils.java

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv() * (x-x_hat)T) / ( 2^0.5d * det()^0.5)
 * /*from   ww  w. j  a v  a  2 s . c o  m*/
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim,
            "|x| != |x_hat|, |x|=" + dim + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim,
            "|x| != |sigma|, |x|=" + dim + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}

From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java

/**
* Performs a robust least squares fit with bisquare weights to the supplied data.
* 
* @param indVarValues A RealVector containing the values of the independent variable.
* @param depVarValues A RealVector containing the values of the dependent variable.
* @return a RealVector containing two elements: the slope of the fit and the y-intercept of the fit.
*///www .  ja v  a 2  s  . com
public RealVector fit(RealVector indVarValues, RealVector depVarValues) {

    RealVector uniformWeights = new ArrayRealVector(indVarValues.getDimension(), 1.0);

    RealVector lastParams = new ArrayRealVector(2, Double.MAX_VALUE);

    RealVector currParams = wlsFit(indVarValues, depVarValues, uniformWeights);

    RealVector weights = uniformWeights;

    RealVector leverages = this.calculateLeverages(indVarValues);

    int c = 0;

    double norm_mult = 1.0;

    if (!this.noIntercept) {
        norm_mult = 2.0;
    }

    int maxiter = 10000;

    while (lastParams.subtract(currParams).getNorm() > CONV_NORM * norm_mult && c++ < maxiter) {

        lastParams = currParams;

        RealVector stdAdjR = this.calculateStandardizedAdjustedResiduals(indVarValues, depVarValues, leverages,
                currParams);

        weights = calculateBisquareWeights(stdAdjR);

        currParams = wlsFit(indVarValues, depVarValues, weights);

    }

    return currParams;

}

From source file:edu.stanford.cfuller.imageanalysistools.fitting.ImageObject.java

/**
 * Gets the vector difference between the position of the object in two channels.
 * //from  w  w  w .  ja v  a 2  s .c o m
 * Note that there is no unit conversion here, and the distance is returned in image units of pixels or sections.
 * 
 * @param channel0   The index of one channel to use for the difference.
 * @param channel1   The index of the other channel to use for the difference.
 * @return         The vector difference between the two channels, as channel1 - channel0, or null if either channel is out of range or has not yet been fit.
 */
public RealVector getVectorDifferenceBetweenChannels(int channel0, int channel1) {

    int key = this.numberOfChannels * channel0 + channel1;

    if (!this.vectorDifferencesBetweenChannels.containsKey(key)) {

        RealVector ch0 = this.getPositionForChannel(channel0);
        RealVector ch1 = this.getPositionForChannel(channel1);

        if (ch0 == null || ch1 == null) {
            return null;
        }

        this.vectorDifferencesBetweenChannels.put(key, ch1.subtract(ch0));
    }

    return this.vectorDifferencesBetweenChannels.get(key);

}

From source file:edu.stanford.cfuller.imageanalysistools.fitting.ImageObject.java

/**
 * Gets the vector difference between the corrected positions of the object in two channels.
 * /*from   ww w .  ja va 2s. c  o  m*/
 * If there is no corrected position for a channel, its uncorrected position is used. 
 * 
 * Note that there is no unit conversion here, and the distance is returned in image units of pixels or sections.
 * 
 * @param channel0   The index of one channel to use for the difference.
 * @param channel1   The index of the other channel to use for the difference.
 * @return         The vector difference between the two channels, as channel1 - channel0, or null if either channel is out of range or has not yet been fit.
 */
public RealVector getCorrectedVectorDifferenceBetweenChannels(int channel0, int channel1) {

    RealVector c0 = this.correctedPositionsByChannel.get(channel0);
    RealVector c1 = this.correctedPositionsByChannel.get(channel1);

    if (c0 == null)
        c0 = this.positionsByChannel.get(channel0);
    if (c1 == null)
        c1 = this.positionsByChannel.get(channel1);

    return c1.subtract(c0);

}

From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java

/**
* Applies an in situ cellular aberration correction to a dataset.
* 
* @param toCorrect a list of ImageObjects to be corrected.
* @param slopes the slopes of a linear fit between the channel pair used to correct and the channel pair being corrected (as calculated by {@link #determineInSituAberrationCorrection()}, for instance).
* @return a List of RealVectors that are the corrected vector distances between each object's images in the channel pair being corrected.
*//*from   ww w . j av  a2  s  . c  o  m*/
public java.util.List<RealVector> applyInSituAberrationCorrection(java.util.List<ImageObject> toCorrect,
        RealVector slopes) {

    int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM);
    int inSituAberrCorrChannel = this.parameters.getIntValueForKey(IN_SITU_ABERR_SECOND_CH_PARAM);
    int measurementChannel = this.parameters.getIntValueForKey(CORR_CH_PARAM);

    java.util.List<RealVector> correctedDifferences = new java.util.ArrayList<RealVector>();

    for (ImageObject iobj : toCorrect) {
        RealVector corrDiff = iobj
                .getCorrectedVectorDifferenceBetweenChannels(referenceChannel, inSituAberrCorrChannel)
                .ebeMultiply(this.pixelToDistanceConversions);
        RealVector exptDiff = iobj
                .getCorrectedVectorDifferenceBetweenChannels(referenceChannel, measurementChannel)
                .ebeMultiply(this.pixelToDistanceConversions);

        RealVector correction = corrDiff.ebeMultiply(slopes);

        correctedDifferences.add(exptDiff.subtract(correction));
    }

    return correctedDifferences;

}

From source file:com.joptimizer.algebra.CholeskyFactorizationTest.java

/**
 * This test shows that the correct check of the inversion accuracy must be done with
 * the scaled residual, not with the simple norm ||A.x-b||
 *//*from  w ww . j  av  a 2s.  c  o m*/
public void testScaledResidual() throws Exception {
    log.debug("testScaledResidual");

    String matrixId = "1";
    double[][] A = Utils
            .loadDoubleMatrixFromFile("factorization" + File.separator + "matrix" + matrixId + ".csv");
    RealMatrix Q = MatrixUtils.createRealMatrix(A);
    int dim = Q.getRowDimension();

    RealVector b = new ArrayRealVector(new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 });

    CholeskyFactorization cs = new CholeskyFactorization(DoubleFactory2D.dense.make(Q.getData()));
    cs.factorize();
    RealVector x = new ArrayRealVector(cs.solve(DoubleFactory1D.dense.make(b.toArray())).toArray());

    //scaledResidual = ||Ax-b||_oo/( ||A||_oo . ||x||_oo + ||b||_oo )
    // with ||x||_oo = max(x[i])
    double scaledResidual = Utils.calculateScaledResidual(Q.getData(), x.toArray(), b.toArray());
    log.debug("scaledResidual: " + scaledResidual);
    assertTrue(scaledResidual < Utils.getDoubleMachineEpsilon());

    //b - A.x
    //checking the simple norm, this will fail
    double n1 = b.subtract(Q.operate(x)).getNorm();
    log.debug("||b - A.x||: " + n1);
    //assertTrue(n1 < 1.E-8);
}

From source file:com.joptimizer.algebra.CholeskySparseFactorizationTest.java

/**
 * This test shows that the correct check of the inversion accuracy must be done with
 * the scaled residual, not with the simple norm ||A.x-b||
 *///  www  . j  a  v  a  2  s .  c o m
public void testScaledResidual() throws Exception {
    log.debug("testScaledResidual");

    String matrixId = "1";
    double[][] A = Utils
            .loadDoubleMatrixFromFile("factorization" + File.separator + "matrix" + matrixId + ".csv");
    RealMatrix Q = MatrixUtils.createRealMatrix(A);
    int dim = Q.getRowDimension();

    RealVector b = new ArrayRealVector(new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 });

    CholeskySparseFactorization cs = new CholeskySparseFactorization(new SparseDoubleMatrix2D(Q.getData()));
    cs.factorize();
    RealVector x = new ArrayRealVector(cs.solve(F1.make(b.toArray())).toArray());

    //scaledResidual = ||Ax-b||_oo/( ||A||_oo . ||x||_oo + ||b||_oo )
    // with ||x||_oo = max(x[i])
    double residual = Utils.calculateScaledResidual(A, x.toArray(), b.toArray());
    log.debug("residual: " + residual);
    assertTrue(residual < Utils.getDoubleMachineEpsilon());

    //b - A.x
    //checking the simple norm, this will fail
    double n1 = b.subtract(Q.operate(x)).getNorm();
    log.debug("||b - A.x||: " + n1);
    //assertTrue(n1 < 1.E-8);
}