org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution Java Examples

@Test(expected = ConvergenceException.class)
public void testConvergenceException() {
    // ConvergenceException thrown if fit terminates before threshold met
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    // 5 iterations not enough to meet convergence threshold
    fitter.fit(initialMix, 5, 1E-5);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testMaxIterationsPositive() {
    // Maximum iterations for fit must be positive integer
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 0, 1E-5);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testThresholdPositive() {
    // Maximum iterations for fit must be positive
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(
                data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 1000, 0);
}
 

@Test(expected = ConvergenceException.class)
public void testConvergenceException() {
    // ConvergenceException thrown if fit terminates before threshold met
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    // 5 iterations not enough to meet convergence threshold
    fitter.fit(initialMix, 5, 1E-5);
}
 

@Test(expected = DimensionMismatchException.class)
public void testIncompatibleIntialMixture() {
    // Data has 3 columns
    double[][] data = new double[][] {
            { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 }
    };
    double[] weights = new double[] { 0.5, 0.5 };

    // These distributions are compatible with 2-column data, not 3-column
    // data
    MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[2];

    mvns[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    mvns[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    // Create components and mixture
    List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[0], mvns[0]));
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[1], mvns[1]));

    MixtureMultivariateNormalDistribution badInitialMix
        = new MixtureMultivariateNormalDistribution(components);

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    fitter.fit(badInitialMix);
}
 

@Test(expected = DimensionMismatchException.class)
public void testIncompatibleIntialMixture() {
    // Data has 3 columns
    double[][] data = new double[][] {
            { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 }
    };
    double[] weights = new double[] { 0.5, 0.5 };

    // These distributions are compatible with 2-column data, not 3-column
    // data
    MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[2];

    mvns[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    mvns[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    // Create components and mixture
    List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[0], mvns[0]));
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[1], mvns[1]));

    MixtureMultivariateNormalDistribution badInitialMix
        = new MixtureMultivariateNormalDistribution(components);

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    fitter.fit(badInitialMix);
}
 

@Test(expected = ConvergenceException.class)
public void testConvergenceException() {
    // ConvergenceException thrown if fit terminates before threshold met
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    // 5 iterations not enough to meet convergence threshold
    fitter.fit(initialMix, 5, 1E-5);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testThresholdPositive() {
    // Maximum iterations for fit must be positive
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(
                data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 1000, 0);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testMaxIterationsPositive() {
    // Maximum iterations for fit must be positive integer
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 0, 1E-5);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testMaxIterationsPositive() {
    // Maximum iterations for fit must be positive integer
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 0, 1E-5);
}
 

@Ignore@Test
public void testInitialMixture() {
    // Testing initial mixture estimated from data
    double[] correctWeights = new double[] { 0.5, 0.5 };

    MultivariateNormalDistribution[] correctMVNs = new MultivariateNormalDistribution[2];

    correctMVNs[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    correctMVNs[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    final MixtureMultivariateNormalDistribution initialMix
        = MultivariateNormalMixtureExpectationMaximization.estimate(getTestSamples(), 2);

    int i = 0;
    for (Pair<Double, MultivariateNormalDistribution> component : initialMix
            .getComponents()) {
        Assert.assertEquals(correctWeights[i], component.getFirst(),
                Math.ulp(1d));
        Assert.assertEquals(correctMVNs[i], component.getSecond());
        i++;
    }
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testThresholdPositive() {
    // Maximum iterations for fit must be positive
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(
                data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 1000, 0);
}
 

@Test(expected = DimensionMismatchException.class)
public void testIncompatibleIntialMixture() {
    // Data has 3 columns
    double[][] data = new double[][] {
            { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 }
    };
    double[] weights = new double[] { 0.5, 0.5 };

    // These distributions are compatible with 2-column data, not 3-column
    // data
    MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[2];

    mvns[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    mvns[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    // Create components and mixture
    List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[0], mvns[0]));
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[1], mvns[1]));

    MixtureMultivariateNormalDistribution badInitialMix
        = new MixtureMultivariateNormalDistribution(components);

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    fitter.fit(badInitialMix);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testThresholdPositive() {
    // Maximum iterations for fit must be positive
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(
                data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 1000, 0);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testMaxIterationsPositive() {
    // Maximum iterations for fit must be positive integer
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 0, 1E-5);
}
 

@Test(expected = ConvergenceException.class)
public void testConvergenceException() {
    // ConvergenceException thrown if fit terminates before threshold met
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    // 5 iterations not enough to meet convergence threshold
    fitter.fit(initialMix, 5, 1E-5);
}
 

@Test(expected = DimensionMismatchException.class)
public void testIncompatibleIntialMixture() {
    // Data has 3 columns
    double[][] data = new double[][] {
            { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 }
    };
    double[] weights = new double[] { 0.5, 0.5 };

    // These distributions are compatible with 2-column data, not 3-column
    // data
    MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[2];

    mvns[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    mvns[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    // Create components and mixture
    List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[0], mvns[0]));
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[1], mvns[1]));

    MixtureMultivariateNormalDistribution badInitialMix
        = new MixtureMultivariateNormalDistribution(components);

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    fitter.fit(badInitialMix);
}
 

/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
 * <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
 * <p>Returns null if no pass filtering.  Please note that in these cases,
 * in the rest of this class, we use this to assume that a GMM is not a good model.</p>
 *
 * @param segments  -- segments with segment mean in log2 copy ratio space
 * @param targets -- targets with a log2 copy ratio estimate
 * @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
 *                      in the evaluation
 * @param numComponents -- number of components to use in the GMM.  Usually, this is 2.
 * @return  never {@code null}.  Fitting result with indications whether it converged or was even attempted.
 */
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments,
                                                                                          final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents){

    // For each target in a segment that contains enough targets, normalize the difference against the segment mean
    //  and collapse the filtered targets into the copy ratio estimates.
    final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);

    if (filteredTargetsSegDiff.size() < numComponents) {
        return new MixtureMultivariateNormalFitResult(null, false, false);
    }

    // Assume that Apache Commons wants data points in the first dimension.
    // Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
    final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];

    // Convert the filtered targets into 2d array (even if second dimension is length 1).  The second dimension is
    //  uncorrelated Gaussian.  This is only to get around funny API in Apache Commons, which will throw an
    //  exception if the length of the second dimension is < 2
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
    final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
    for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
        filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
        filteredTargetsSegDiff2d[i][1] = nd.sample();
    }

    final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
    final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);

    try {
        multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
    } catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
        // We are done, we cannot make a fitting.  We should return a result that we attempted a fitting, but it
        //  did not converge.  Include the model as it was when the exception was thrown.
        return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
    }
    return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testMaxIterationsPositive() {
    // Maximum iterations for fit must be positive integer
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 0, 1E-5);
}
 

@Test(expected = NotStrictlyPositiveException.class)
public void testThresholdPositive() {
    // Maximum iterations for fit must be positive
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter =
            new MultivariateNormalMixtureExpectationMaximization(
                data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    fitter.fit(initialMix, 1000, 0);
}
 

public double density(MixtureMultivariateNormalDistribution dist, double[] vals) {
    double dens = 0;
    for (int g=0; g<dist.getComponents().size(); g++) {
        double td = dist.getComponents().get(g).getSecond().density(vals);
        if (td>dens) {
            dens = td;
        }
    }
    return dens;
}
 

public double classify(MixtureMultivariateNormalDistribution dist, double[] vals) {
    double dens = 0;
    double c = 0;
    for (int g=0; g<dist.getComponents().size(); g++) {
        double td = dist.getComponents().get(g).getSecond().density(vals);
        if (td>dens) {
            c = g;
            dens = td;
        }
    }
    return c;
}
 

@Test(expected = ConvergenceException.class)
public void testConvergenceException() {
    // ConvergenceException thrown if fit terminates before threshold met
    double[][] data = getTestSamples();
    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution
        initialMix = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);

    // 5 iterations not enough to meet convergence threshold
    fitter.fit(initialMix, 5, 1E-5);
}
 

@Test(expected = DimensionMismatchException.class)
public void testIncompatibleIntialMixture() {
    // Data has 3 columns
    double[][] data = new double[][] {
            { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 }
    };
    double[] weights = new double[] { 0.5, 0.5 };

    // These distributions are compatible with 2-column data, not 3-column
    // data
    MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[2];

    mvns[0] = new MultivariateNormalDistribution(new double[] {
                    -0.0021722935000328823, 3.5432892936887908 },
                    new double[][] {
                            { 4.537422569229048, 3.5266152281729304 },
                            { 3.5266152281729304, 6.175448814169779 } });
    mvns[1] = new MultivariateNormalDistribution(new double[] {
                    5.090902706507635, 8.68540656355283 }, new double[][] {
                    { 2.886778573963039, 1.5257474543463154 },
                    { 1.5257474543463154, 3.3794567673616918 } });

    // Create components and mixture
    List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[0], mvns[0]));
    components.add(new Pair<Double, MultivariateNormalDistribution>(
            weights[1], mvns[1]));

    MixtureMultivariateNormalDistribution badInitialMix
        = new MixtureMultivariateNormalDistribution(components);

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    fitter.fit(badInitialMix);
}
 

@Test
public void testFit() {
    // Test that the loglikelihood, weights, and models are determined and
    // fitted correctly
    final double[][] data = getTestSamples();
    final double correctLogLikelihood = -4.292431006791994;
    final double[] correctWeights = new double[] { 0.2962324189652912, 0.7037675810347089 };
    
    final double[][] correctMeans = new double[][]{
        {-1.4213112715121132, 1.6924690505757753},
        {4.213612224374709, 7.975621325853645}
    };
    
    final RealMatrix[] correctCovMats = new Array2DRowRealMatrix[2];
    correctCovMats[0] = new Array2DRowRealMatrix(new double[][] {
        { 1.739356907285747, -0.5867644251487614 },
        { -0.5867644251487614, 1.0232932029324642 } }
            );
    correctCovMats[1] = new Array2DRowRealMatrix(new double[][] {
        { 4.245384898007161, 2.5797798966382155 },
        { 2.5797798966382155, 3.9200272522448367 } });
    
    final MultivariateNormalDistribution[] correctMVNs = new MultivariateNormalDistribution[2];
    correctMVNs[0] = new MultivariateNormalDistribution(correctMeans[0], correctCovMats[0].getData());
    correctMVNs[1] = new MultivariateNormalDistribution(correctMeans[1], correctCovMats[1].getData());

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution initialMix
        = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);
    fitter.fit(initialMix);
    MixtureMultivariateNormalDistribution fittedMix = fitter.getFittedModel();
    List<Pair<Double, MultivariateNormalDistribution>> components = fittedMix.getComponents();

    Assert.assertEquals(correctLogLikelihood,
                        fitter.getLogLikelihood(),
                        Math.ulp(1d));

    int i = 0;
    for (Pair<Double, MultivariateNormalDistribution> component : components) {
        final double weight = component.getFirst();
        final MultivariateNormalDistribution mvn = component.getSecond();
        final double[] mean = mvn.getMeans();
        final RealMatrix covMat = mvn.getCovariances();
        Assert.assertEquals(correctWeights[i], weight, Math.ulp(1d));
        Assert.assertTrue(Arrays.equals(correctMeans[i], mean));
        Assert.assertEquals(correctCovMats[i], covMat);
        i++;
    }
}
 

/**
 * Helper method to create a multivariate normal mixture model which can be
 * used to initialize {@link #fit(MixtureMultivariateNormalDistribution)}.
 *
 * This method uses the data supplied to the constructor to try to determine
 * a good mixture model at which to start the fit, but it is not guaranteed
 * to supply a model which will find the optimal solution or even converge.
 *
 * @param data Data to estimate distribution
 * @param numComponents Number of components for estimated mixture
 * @return Multivariate normal mixture model estimated from the data
 * @throws NumberIsTooLargeException if {@code numComponents} is greater
 * than the number of data rows.
 * @throws NumberIsTooSmallException if {@code numComponents < 2}.
 * @throws NotStrictlyPositiveException if data has less than 2 rows
 * @throws DimensionMismatchException if rows of data have different numbers
 *             of columns
 */
public static MixtureMultivariateNormalDistribution estimate(final double[][] data,
                                                             final int numComponents)
    throws NotStrictlyPositiveException,
           DimensionMismatchException {
    if (data.length < 2) {
        throw new NotStrictlyPositiveException(data.length);
    }
    if (numComponents < 2) {
        throw new NumberIsTooSmallException(numComponents, 2, true);
    }
    if (numComponents > data.length) {
        throw new NumberIsTooLargeException(numComponents, data.length, true);
    }

    final int numRows = data.length;
    final int numCols = data[0].length;

    // sort the data
    final DataRow[] sortedData = new DataRow[numRows];
    for (int i = 0; i < numRows; i++) {
        sortedData[i] = new DataRow(data[i]);
    }
    Arrays.sort(sortedData);

    // uniform weight for each bin
    final double weight = 1d / numComponents;

    // components of mixture model to be created
    final List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>(numComponents);

    // create a component based on data in each bin
    for (int binIndex = 0; binIndex < numComponents; binIndex++) {
        // minimum index (inclusive) from sorted data for this bin
        final int minIndex = (binIndex * numRows) / numComponents;

        // maximum index (exclusive) from sorted data for this bin
        final int maxIndex = ((binIndex + 1) * numRows) / numComponents;

        // number of data records that will be in this bin
        final int numBinRows = maxIndex - minIndex;

        // data for this bin
        final double[][] binData = new double[numBinRows][numCols];

        // mean of each column for the data in the this bin
        final double[] columnMeans = new double[numCols];

        // populate bin and create component
        for (int i = minIndex, iBin = 0; i < maxIndex; i++, iBin++) {
            for (int j = 0; j < numCols; j++) {
                final double val = sortedData[i].getRow()[j];
                columnMeans[j] += val;
                binData[iBin][j] = val;
            }
        }

        MathArrays.scaleInPlace(1d / numBinRows, columnMeans);

        // covariance matrix for this bin
        final double[][] covMat
            = new Covariance(binData).getCovarianceMatrix().getData();
        final MultivariateNormalDistribution mvn
            = new MultivariateNormalDistribution(columnMeans, covMat);

        components.add(new Pair<Double, MultivariateNormalDistribution>(weight, mvn));
    }

    return new MixtureMultivariateNormalDistribution(components);
}
 

@Ignore@Test
public void testFit() {
    // Test that the loglikelihood, weights, and models are determined and
    // fitted correctly
    double[][] data = getTestSamples();
    double correctLogLikelihood = -4.292431006791994;
    double[] correctWeights = new double[] { 0.2962324189652912, 0.7037675810347089 };
    MultivariateNormalDistribution[] correctMVNs = new MultivariateNormalDistribution[2];
    correctMVNs[0] = new MultivariateNormalDistribution(new double[] {
                    -1.4213112715121132, 1.6924690505757753 },
                    new double[][] {
                            { 1.739356907285747, -0.5867644251487614 },
                            { -0.5867644251487614, 1.0232932029324642 } });

    correctMVNs[1] = new MultivariateNormalDistribution(new double[] {
                    4.213612224374709, 7.975621325853645 },
                    new double[][] {
                            { 4.245384898007161, 2.5797798966382155 },
                            { 2.5797798966382155, 3.9200272522448367 } });

    MultivariateNormalMixtureExpectationMaximization fitter
        = new MultivariateNormalMixtureExpectationMaximization(data);

    MixtureMultivariateNormalDistribution initialMix
        = MultivariateNormalMixtureExpectationMaximization.estimate(data, 2);
    fitter.fit(initialMix);
    MixtureMultivariateNormalDistribution fittedMix = fitter.getFittedModel();
    List<Pair<Double, MultivariateNormalDistribution>> components = fittedMix.getComponents();

    Assert.assertEquals(correctLogLikelihood,
                        fitter.getLogLikelihood(),
                        Math.ulp(1d));

    int i = 0;
    for (Pair<Double, MultivariateNormalDistribution> component : components) {
        double weight = component.getFirst();
        MultivariateNormalDistribution mvn = component.getSecond();
        Assert.assertEquals(correctWeights[i], weight, Math.ulp(1d));
        Assert.assertEquals(correctMVNs[i], mvn);
        i++;
    }
}
 

/**
 * Helper method to create a multivariate normal mixture model which can be
 * used to initialize {@link #fit(MixtureMultivariateNormalDistribution)}.
 *
 * This method uses the data supplied to the constructor to try to determine
 * a good mixture model at which to start the fit, but it is not guaranteed
 * to supply a model which will find the optimal solution or even converge.
 *
 * @param data Data to estimate distribution
 * @param numComponents Number of components for estimated mixture
 * @return Multivariate normal mixture model estimated from the data
 * @throws NumberIsTooLargeException if {@code numComponents} is greater
 * than the number of data rows.
 * @throws NumberIsTooSmallException if {@code numComponents < 2}.
 * @throws NotStrictlyPositiveException if data has less than 2 rows
 * @throws DimensionMismatchException if rows of data have different numbers
 *             of columns
 */
public static MixtureMultivariateNormalDistribution estimate(final double[][] data,
                                                             final int numComponents)
    throws NotStrictlyPositiveException,
           DimensionMismatchException {
    if (data.length < 2) {
        throw new NotStrictlyPositiveException(data.length);
    }
    if (numComponents < 2) {
        throw new NumberIsTooSmallException(numComponents, 2, true);
    }
    if (numComponents > data.length) {
        throw new NumberIsTooLargeException(numComponents, data.length, true);
    }

    final int numRows = data.length;
    final int numCols = data[0].length;

    // sort the data
    final DataRow[] sortedData = new DataRow[numRows];
    for (int i = 0; i < numRows; i++) {
        sortedData[i] = new DataRow(data[i]);
    }
    Arrays.sort(sortedData);

    // uniform weight for each bin
    final double weight = 1d / numComponents;

    // components of mixture model to be created
    final List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();

    // create a component based on data in each bin
    for (int binIndex = 0; binIndex < numComponents; binIndex++) {
        // minimum index (inclusive) from sorted data for this bin
        final int minIndex = (binIndex * numRows) / numComponents;

        // maximum index (exclusive) from sorted data for this bin
        final int maxIndex = ((binIndex + 1) * numRows) / numComponents;

        // number of data records that will be in this bin
        final int numBinRows = maxIndex - minIndex;

        // data for this bin
        final double[][] binData = new double[numBinRows][numCols];

        // mean of each column for the data in the this bin
        final double[] columnMeans = new double[numCols];

        // populate bin and create component
        for (int i = minIndex, iBin = 0; i < maxIndex; i++, iBin++) {
            for (int j = 0; j < numCols; j++) {
                final double val = sortedData[i].getRow()[j];
                columnMeans[j] += val;
                binData[iBin][j] = val;
            }
        }

        MathArrays.scaleInPlace(1d / numBinRows, columnMeans);

        // covariance matrix for this bin
        final double[][] covMat
            = new Covariance(binData).getCovarianceMatrix().getData();
        final MultivariateNormalDistribution mvn
            = new MultivariateNormalDistribution(columnMeans, covMat);

        components.add(new Pair<Double, MultivariateNormalDistribution>(weight, mvn));
    }

    return new MixtureMultivariateNormalDistribution(components);
}
 

@Test
public void testInitialMixture() {
    // Testing initial mixture estimated from data
    final double[] correctWeights = new double[] { 0.5, 0.5 };

    final double[][] correctMeans = new double[][] {
        {-0.0021722935000328823, 3.5432892936887908},
        {5.090902706507635, 8.68540656355283},
    };

    final RealMatrix[] correctCovMats = new Array2DRowRealMatrix[2];

    correctCovMats[0] = new Array2DRowRealMatrix(new double[][] {
            { 4.537422569229048, 3.5266152281729304 },
            { 3.5266152281729304, 6.175448814169779 } });

    correctCovMats[1] = new Array2DRowRealMatrix( new double[][] {
            { 2.886778573963039, 1.5257474543463154 },
            { 1.5257474543463154, 3.3794567673616918 } });

    final MultivariateNormalDistribution[] correctMVNs = new
            MultivariateNormalDistribution[2];

    correctMVNs[0] = new MultivariateNormalDistribution(correctMeans[0],
            correctCovMats[0].getData());

    correctMVNs[1] = new MultivariateNormalDistribution(correctMeans[1],
            correctCovMats[1].getData());

    final MixtureMultivariateNormalDistribution initialMix
        = MultivariateNormalMixtureExpectationMaximization.estimate(getTestSamples(), 2);

    int i = 0;
    for (Pair<Double, MultivariateNormalDistribution> component : initialMix
            .getComponents()) {
        Assert.assertEquals(correctWeights[i], component.getFirst(),
                Math.ulp(1d));
        
        final double[] means = component.getValue().getMeans();
        Assert.assertTrue(Arrays.equals(correctMeans[i], means));
        
        final RealMatrix covMat = component.getValue().getCovariances();
        Assert.assertEquals(correctCovMats[i], covMat);
        i++;
    }
}
 

/**
 * Helper method to create a multivariate normal mixture model which can be
 * used to initialize {@link #fit(MixtureMultivariateRealDistribution)}.
 *
 * This method uses the data supplied to the constructor to try to determine
 * a good mixture model at which to start the fit, but it is not guaranteed
 * to supply a model which will find the optimal solution or even converge.
 *
 * @param data Data to estimate distribution
 * @param numComponents Number of components for estimated mixture
 * @return Multivariate normal mixture model estimated from the data
 * @throws NumberIsTooLargeException if {@code numComponents\ is greater
 * than the number of data rows.
 * @throws NumberIsTooSmallException if {@code numComponents < 2}.
 * @throws NotStrictlyPositiveException if data has less than 2 rows
 * @throws DimensionMismatchException if rows of data have different numbers
 *             of columns
 * @see #fit
 */
public static MixtureMultivariateNormalDistribution estimate(final double[][] data,
                                                             final int numComponents)
    throws NotStrictlyPositiveException,
           DimensionMismatchException {
    if (data.length < 2) {
        throw new NotStrictlyPositiveException(data.length);
    }
    if (numComponents < 2) {
        throw new NumberIsTooSmallException(numComponents, 2, true);
    }
    if (numComponents > data.length) {
        throw new NumberIsTooLargeException(numComponents, data.length, true);
    }

    final int numRows = data.length;
    final int numCols = data[0].length;

    // sort the data
    final DataRow[] sortedData = new DataRow[numRows];
    for (int i = 0; i < numRows; i++) {
        sortedData[i] = new DataRow(data[i]);
    }
    Arrays.sort(sortedData);

    final int totalBins = numComponents;

    // uniform weight for each bin
    final double weight = 1d / totalBins;

    // components of mixture model to be created
    final List<Pair<Double, MultivariateNormalDistribution>> components =
            new ArrayList<Pair<Double, MultivariateNormalDistribution>>();

    // create a component based on data in each bin
    for (int binNumber = 1; binNumber <= totalBins; binNumber++) {
        // minimum index from sorted data for this bin
        final int minIndex
            = (int) FastMath.max(0,
                                 FastMath.floor((binNumber - 1) * numRows / totalBins));

        // maximum index from sorted data for this bin
        final int maxIndex
            = (int) FastMath.ceil(binNumber * numRows / numComponents) - 1;

        // number of data records that will be in this bin
        final int numBinRows = maxIndex - minIndex + 1;

        // data for this bin
        final double[][] binData = new double[numBinRows][numCols];

        // mean of each column for the data in the this bin
        final double[] columnMeans = new double[numCols];

        // populate bin and create component
        for (int i = minIndex, iBin = 0; i <= maxIndex; i++, iBin++) {
            for (int j = 0; j < numCols; j++) {
                final double val = sortedData[i].getRow()[j];
                columnMeans[j] += val;
                binData[iBin][j] = val;
            }
        }

        MathArrays.scaleInPlace(1d / numBinRows, columnMeans);

        // covariance matrix for this bin
        final double[][] covMat
            = new Covariance(binData).getCovarianceMatrix().getData();
        final MultivariateNormalDistribution mvn
            = new MultivariateNormalDistribution(columnMeans, covMat);

        components.add(new Pair<Double, MultivariateNormalDistribution>(weight, mvn));
    }

    return new MixtureMultivariateNormalDistribution(components);
}