MultivariateNormalMixtureExpectationMaximization { 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; final DataRow[] sortedData = new DataRow[numRows]; for (int i = 0; i < numRows; i++) { sortedData[i] = new DataRow(data[i]); } Arrays.sort(sortedData); final double weight = 1d / numComponents; final List<Pair<Double, MultivariateNormalDistribution>> components = new ArrayList<Pair<Double, MultivariateNormalDistribution>>(numComponents); for (int binIndex = 0; binIndex < numComponents; binIndex++) { final int minIndex = (binIndex * numRows) / numComponents; final int maxIndex = ((binIndex + 1) * numRows) / numComponents; final int numBinRows = maxIndex - minIndex; final double[][] binData = new double[numBinRows][numCols]; final double[] columnMeans = new double[numCols]; 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); 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); } MultivariateNormalMixtureExpectationMaximization(double[][] data); void fit(final MixtureMultivariateNormalDistribution initialMixture,
                    final int maxIterations,
                    final double threshold); void fit(MixtureMultivariateNormalDistribution initialMixture); static MixtureMultivariateNormalDistribution estimate(final double[][] data,
                                                                 final int numComponents); double getLogLikelihood(); MixtureMultivariateNormalDistribution getFittedModel();  }
@Test public void testInitialMixture() { 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++; } }