Example usage for org.apache.commons.math3.stat.correlation StorelessCovariance getCovarianceMatrix

List of usage examples for org.apache.commons.math3.stat.correlation StorelessCovariance getCovarianceMatrix

Introduction

In this page you can find the example usage for org.apache.commons.math3.stat.correlation StorelessCovariance getCovarianceMatrix.

Prototype

@Override
public RealMatrix getCovarianceMatrix() throws NumberIsTooSmallException 

Source Link

Usage

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

/**
 * @param data The elements of 'data' will be modified.
 * @param label/*from   w w w.ja  v a2  s .c  o m*/
 * @param Nclasses
 * @param shrinkage_intensity
 */
public LinearDiscriminantAnalysis(final ArrayList<double[]> data, final int[] label, final int Nclasses,
        final double shrinkage) {
    assert (data.size() == label.length);

    final int Ndata = data.size();
    final int Ndim = data.get(0).length;

    // Partition data by class
    final ArrayList<ArrayList<double[]>> classes = new ArrayList<ArrayList<double[]>>(Nclasses);
    for (int i = 0; i < Nclasses; ++i) {
        classes.add(new ArrayList<double[]>());
    }
    for (int i = 0; i < data.size(); ++i) {
        classes.get(label[i]).add(data.get(i));
    }

    // Mean center the data

    final VectorMeanVarianceAccumulator mv = new VectorMeanVarianceAccumulator(Ndim);
    for (int i = 0; i < Ndata; ++i) {
        mv.add(data.get(i));
    }
    mean = mv.mean();
    // Subtract global mean
    for (final double[] x : data) {
        Fn.vminus_inplace(x, mean);
    }

    // Calculate class means and covariances
    final double[][] class_mean = new double[Nclasses][Ndim];
    final RealMatrix[] class_cov = new RealMatrix[Nclasses];

    for (int i = 0; i < Nclasses; ++i) {
        final ArrayList<double[]> Xc = classes.get(i);
        final VectorMeanVarianceAccumulator mv_i = new VectorMeanVarianceAccumulator(Ndim);
        final StorelessCovariance cov = new StorelessCovariance(Ndim);
        for (int j = 0; j < Xc.size(); ++j) {
            final double[] x = Xc.get(j);
            mv_i.add(x);
            cov.increment(x);
        }
        class_mean[i] = mv_i.mean();
        class_cov[i] = cov.getCovarianceMatrix();
    }

    // Between-class scatter.
    // Note that 'data' is mean-centered, so the global mean is 0.

    RealMatrix Sb_builder = new Array2DRowRealMatrix(Ndim, Ndim);
    for (int i = 0; i < Nclasses; ++i) {
        final RealVector mu_i = new ArrayRealVector(class_mean[i]);
        final RealMatrix xxt = mu_i.outerProduct(mu_i);
        Sb_builder = Sb_builder.add(xxt.scalarMultiply(classes.get(i).size() / ((double) Ndata - 1)));
    }
    this.Sb = Sb_builder;
    Sb_builder = null;

    // Within-class scatter with shrinkage estimate:
    // Sw = (1.0 - shrinkage) * \sum Sigma_i + shrinkage * I

    RealMatrix Sw_builder = new Array2DRowRealMatrix(Ndim, Ndim);
    for (int i = 0; i < Nclasses; ++i) {
        final RealMatrix Sigma_i = class_cov[i];
        final RealMatrix scaled = Sigma_i.scalarMultiply((1.0 - shrinkage) * (classes.get(i).size() - 1));
        Sw_builder = Sw_builder.add(scaled);
    }
    for (int i = 0; i < Ndim; ++i) {
        Sw_builder.setEntry(i, i, Sw_builder.getEntry(i, i) + shrinkage);
    }
    this.Sw = Sw_builder;
    Sw_builder = null;

    // Invert Sw
    System.out.println("[LDA] Sw inverse");
    final RealMatrix Sw_inv = new LUDecomposition(Sw).getSolver().getInverse();
    final RealMatrix F = Sw_inv.multiply(Sb);

    System.out.println("[LDA] Eigendecomposition");
    eigenvalues = new double[Nclasses - 1];
    eigenvectors = new ArrayList<RealVector>(Nclasses - 1);
    final EigenDecomposition evd = new EigenDecomposition(F);
    for (int j = 0; j < Nclasses - 1; ++j) {
        final double eigenvalue = evd.getRealEigenvalue(j);
        eigenvalues[j] = eigenvalue;
        //         final double scale = 1.0 / Math.sqrt( eigenvalue );
        //         eigenvectors.add( evd.getEigenvector( j ).mapMultiply( scale ) );
        eigenvectors.add(evd.getEigenvector(j));
    }
}

From source file:com.analog.lyric.dimple.test.solvers.sumproduct.TestSampledFactors.java

/**
 * Adapted from MATLAB test4 in tests/algoGaussian/testSampledFactors.m
 *//*from  w ww .  ja v  a 2 s. c o m*/
@Test
public void sampledComplexProduct() {
    // NOTE: test may fail if seed is changed! We keep the number of samples down so that the test doesn't
    // take too long. Increasing the samples produces better results.

    testRand.setSeed(42);

    try (CurrentModel cur = using(new FactorGraph())) {
        final Complex a = complex("a");
        final Complex b = complex("b");
        final Complex c = product(a, b);

        double[] aMean = new double[] { 10, 10 };
        RealMatrix aCovariance = randCovariance(2);
        a.setPrior(new MultivariateNormal(aMean, aCovariance.getData()));

        double[] bMean = new double[] { -20, 20 };
        RealMatrix bCovariance = randCovariance(2);
        b.setPrior(new MultivariateNormalParameters(bMean, bCovariance.getData()));

        GaussianRandomGenerator normalGenerator = new GaussianRandomGenerator(testRand);
        CorrelatedRandomVectorGenerator aGenerator = new CorrelatedRandomVectorGenerator(aMean, aCovariance,
                1e-12, normalGenerator);
        CorrelatedRandomVectorGenerator bGenerator = new CorrelatedRandomVectorGenerator(bMean, bCovariance,
                1e-12, normalGenerator);

        StorelessCovariance expectedCov = new StorelessCovariance(2);

        final int nSamples = 10000;

        RealVector expectedMean = MatrixUtils.createRealVector(new double[2]);
        double[] cSample = new double[2];

        for (int i = 0; i < nSamples; ++i) {
            double[] aSample = aGenerator.nextVector();
            double[] bSample = bGenerator.nextVector();

            // Compute complex product
            cSample[0] = aSample[0] * bSample[0] - aSample[1] * bSample[1];
            cSample[1] = aSample[0] * bSample[1] + aSample[1] * bSample[0];

            expectedMean.addToEntry(0, cSample[0]);
            expectedMean.addToEntry(1, cSample[1]);

            expectedCov.increment(cSample);
        }

        expectedMean.mapDivideToSelf(nSamples); // normalize

        SumProductSolverGraph sfg = requireNonNull(cur.graph.setSolverFactory(new SumProductSolver()));
        sfg.setOption(GibbsOptions.numSamples, nSamples);

        sfg.solve();

        MultivariateNormalParameters cBelief = requireNonNull(c.getBelief());

        RealVector observedMean = MatrixUtils.createRealVector(cBelief.getMean());
        double scaledMeanDistance = expectedMean.getDistance(observedMean) / expectedMean.getNorm();

        //         System.out.format("expectedMean = %s\n", expectedMean);
        //         System.out.format("observedMean = %s\n", observedMean);
        //         System.out.println(scaledMeanDistance);

        assertEquals(0.0, scaledMeanDistance, .02);

        RealMatrix expectedCovariance = expectedCov.getCovarianceMatrix();
        RealMatrix observedCovariance = MatrixUtils.createRealMatrix(cBelief.getCovariance());
        RealMatrix diffCovariance = expectedCovariance.subtract(observedCovariance);

        double scaledCovarianceDistance = diffCovariance.getNorm() / expectedCovariance.getNorm();

        //         System.out.println(expectedCovariance);
        //         System.out.println(expectedCovariance.getNorm());
        //         System.out.println(diffCovariance);
        //         System.out.println(diffCovariance.getNorm());
        //         System.out.println(diffCovariance.getNorm() / expectedCovariance.getNorm());

        assertEquals(0.0, scaledCovarianceDistance, .2);
    }
}