MultivariateNormalMixtureExpectationMaximization { public void fit(final MixtureMultivariateNormalDistribution initialMixture, final int maxIterations, final double threshold) throws SingularMatrixException, NotStrictlyPositiveException, DimensionMismatchException { if (maxIterations < 1) { throw new NotStrictlyPositiveException(maxIterations); } if (threshold < Double.MIN_VALUE) { throw new NotStrictlyPositiveException(threshold); } final int n = data.length; final int numCols = data[0].length; final int k = initialMixture.getComponents().size(); final int numMeanColumns = initialMixture.getComponents().get(0).getSecond().getMeans().length; if (numMeanColumns != numCols) { throw new DimensionMismatchException(numMeanColumns, numCols); } int numIterations = 0; double previousLogLikelihood = 0d; logLikelihood = Double.NEGATIVE_INFINITY; fittedModel = new MixtureMultivariateNormalDistribution(initialMixture.getComponents()); while (numIterations++ <= maxIterations && FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { previousLogLikelihood = logLikelihood; double sumLogLikelihood = 0d; final List<Pair<Double, MultivariateNormalDistribution>> components = fittedModel.getComponents(); final double[] weights = new double[k]; final MultivariateNormalDistribution[] mvns = new MultivariateNormalDistribution[k]; for (int j = 0; j < k; j++) { weights[j] = components.get(j).getFirst(); mvns[j] = components.get(j).getSecond(); } final double[][] gamma = new double[n][k]; final double[] gammaSums = new double[k]; final double[][] gammaDataProdSums = new double[k][numCols]; for (int i = 0; i < n; i++) { final double rowDensity = fittedModel.density(data[i]); sumLogLikelihood += FastMath.log(rowDensity); for (int j = 0; j < k; j++) { gamma[i][j] = weights[j] * mvns[j].density(data[i]) / rowDensity; gammaSums[j] += gamma[i][j]; for (int col = 0; col < numCols; col++) { gammaDataProdSums[j][col] += gamma[i][j] * data[i][col]; } } } logLikelihood = sumLogLikelihood / n; final double[] newWeights = new double[k]; final double[][] newMeans = new double[k][numCols]; for (int j = 0; j < k; j++) { newWeights[j] = gammaSums[j] / n; for (int col = 0; col < numCols; col++) { newMeans[j][col] = gammaDataProdSums[j][col] / gammaSums[j]; } } final RealMatrix[] newCovMats = new RealMatrix[k]; for (int j = 0; j < k; j++) { newCovMats[j] = new Array2DRowRealMatrix(numCols, numCols); } for (int i = 0; i < n; i++) { for (int j = 0; j < k; j++) { final RealMatrix vec = new Array2DRowRealMatrix(MathArrays.ebeSubtract(data[i], newMeans[j])); final RealMatrix dataCov = vec.multiply(vec.transpose()).scalarMultiply(gamma[i][j]); newCovMats[j] = newCovMats[j].add(dataCov); } } final double[][][] newCovMatArrays = new double[k][numCols][numCols]; for (int j = 0; j < k; j++) { newCovMats[j] = newCovMats[j].scalarMultiply(1d / gammaSums[j]); newCovMatArrays[j] = newCovMats[j].getData(); } fittedModel = new MixtureMultivariateNormalDistribution(newWeights, newMeans, newCovMatArrays); } if (FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { throw new ConvergenceException(); } } 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(expected = DimensionMismatchException.class) public void testIncompatibleIntialMixture() { double[][] data = new double[][] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } }; double[] weights = new double[] { 0.5, 0.5 }; 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 } }); 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() { 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++; } }