Examples with DoubleDoubleBiPredicate uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate used on opensource projects

Example 1 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method createMixed.

/**
 * Creates the multivariate gaussian mixture as the best of many repeats of the expectation
 * maximisation algorithm.
 *
 * @param data the data
 * @param dimensions the dimensions
 * @param numComponents the number of components
 * @return the multivariate gaussian mixture expectation maximization
 */
private MultivariateGaussianMixtureExpectationMaximization createMixed(final double[][] data, int dimensions, int numComponents) {
    // Fit a mixed multivariate Gaussian with different repeats.
    final UnitSphereSampler sampler = UnitSphereSampler.of(UniformRandomProviders.create(Mixers.stafford13(settings.seed++)), dimensions);
    final LocalList<CompletableFuture<MultivariateGaussianMixtureExpectationMaximization>> results = new LocalList<>(settings.repeats);
    final DoubleDoubleBiPredicate test = createConvergenceTest(settings.relativeError);
    if (settings.debug) {
        ImageJUtils.log("  Fitting %d components", numComponents);
    }
    final Ticker ticker = ImageJUtils.createTicker(settings.repeats, 2, "Fitting...");
    final AtomicInteger failures = new AtomicInteger();
    for (int i = 0; i < settings.repeats; i++) {
        final double[] vector = sampler.sample();
        results.add(CompletableFuture.supplyAsync(() -> {
            final MultivariateGaussianMixtureExpectationMaximization fitter = new MultivariateGaussianMixtureExpectationMaximization(data);
            try {
                // This may also throw the same exceptions due to inversion of the covariance matrix
                final MixtureMultivariateGaussianDistribution initialMixture = MultivariateGaussianMixtureExpectationMaximization.estimate(data, numComponents, point -> {
                    double dot = 0;
                    for (int j = 0; j < dimensions; j++) {
                        dot += vector[j] * point[j];
                    }
                    return dot;
                });
                final boolean result = fitter.fit(initialMixture, settings.maxIterations, test);
                // Log the result. Note: The ImageJ log is synchronized.
                if (settings.debug) {
                    ImageJUtils.log("  Fit: log-likelihood=%s, iter=%d, converged=%b", fitter.getLogLikelihood(), fitter.getIterations(), result);
                }
                return result ? fitter : null;
            } catch (NonPositiveDefiniteMatrixException | SingularMatrixException ex) {
                failures.getAndIncrement();
                if (settings.debug) {
                    ImageJUtils.log("  Fit failed during iteration %d. No variance in a sub-population " + "component (check alpha is not always 1.0).", fitter.getIterations());
                }
            } finally {
                ticker.tick();
            }
            return null;
        }));
    }
    ImageJUtils.finished();
    if (failures.get() != 0 && settings.debug) {
        ImageJUtils.log("  %d component fit failed %d/%d", numComponents, failures.get(), settings.repeats);
    }
    // Collect results and return the best model.
    return results.stream().map(f -> f.join()).filter(f -> f != null).sorted((f1, f2) -> Double.compare(f2.getLogLikelihood(), f1.getLogLikelihood())).findFirst().orElse(null);
}

Example 2 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class PoissonFunctionTest method probabilityMatchesPoissonWithNoGain.

private static void probabilityMatchesPoissonWithNoGain(final double mu) {
    final double o = mu;
    final PoissonFunction f = new PoissonFunction(1.0);
    final PoissonDistribution pd = new PoissonDistribution(mu);
    final double p = 0;
    final int[] range = getRange(1, mu);
    final int min = range[0];
    final int max = range[1];
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-8, 0);
    for (int x = min; x <= max; x++) {
        final double v1 = f.likelihood(x, o);
        final double v2 = pd.probability(x);
        TestAssertions.assertTest(v1, v2, predicate, FunctionUtils.getSupplier("g=%f, mu=%f, x=%d", gain, mu, x));
    }
}

Example 3 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class MultivariateGaussianMixtureExpectationMaximizationTest method canCreateMixedMultivariateGaussianDistribution.

@SeededTest
void canCreateMixedMultivariateGaussianDistribution(RandomSeed seed) {
    // Will be normalised
    final double[] weights = { 1, 1 };
    final double[][] means = new double[2][];
    final double[][][] covariances = new double[2][][];
    final double[][] data1 = { { 1, 2 }, { 2.5, 1.5 }, { 3.5, 1.0 } };
    final double[][] data2 = { { 4, 2 }, { 3.5, -1.5 }, { -3.5, 1.0 } };
    means[0] = getColumnMeans(data1);
    covariances[0] = getCovariance(data1);
    means[1] = getColumnMeans(data2);
    covariances[1] = getCovariance(data2);
    // Create components. This does not have to be zero based.
    final LocalList<double[]> list = new LocalList<>();
    list.addAll(Arrays.asList(data1));
    list.addAll(Arrays.asList(data2));
    final double[][] data = list.toArray(new double[0][]);
    final int[] components = { -1, -1, -1, 3, 3, 3 };
    // Randomise the data
    for (int n = 0; n < 3; n++) {
        final long start = n + seed.getSeedAsLong();
        // This relies on the shuffle being the same
        RandomUtils.shuffle(data, RngUtils.create(start));
        RandomUtils.shuffle(components, RngUtils.create(start));
        final MixtureMultivariateGaussianDistribution dist = MultivariateGaussianMixtureExpectationMaximization.createMixed(data, components);
        Assertions.assertArrayEquals(new double[] { 0.5, 0.5 }, dist.getWeights());
        final MultivariateGaussianDistribution[] distributions = dist.getDistributions();
        Assertions.assertEquals(weights.length, distributions.length);
        final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-8);
        for (int i = 0; i < means.length; i++) {
            TestAssertions.assertArrayTest(means[i], distributions[i].getMeans(), test);
            TestAssertions.assertArrayTest(covariances[i], distributions[i].getCovariances(), test);
        }
    }
}

Example 4 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class MultivariateGaussianMixtureExpectationMaximization method fit.

/**
 * Fit a mixture model to the data supplied to the constructor.
 *
 * <p>The quality of the fit depends on the concavity of the data provided to the constructor and
 * the initial mixture provided to this function. If the data has many local optima, multiple runs
 * of the fitting function with different initial mixtures may be required to find the optimal
 * solution. If a SingularMatrixException is encountered, it is possible that another
 * initialisation would work.
 *
 * @param initialMixture model containing initial values of weights and multivariate normals
 * @param maxIterations maximum iterations allowed for fit
 * @param convergencePredicate convergence predicated used to test the logLikelihoods between
 *        successive iterations
 * @return true if converged within the threshold
 * @throws IllegalArgumentException if maxIterations is less than one or if initialMixture mean
 *         vector and data number of columns are not equal
 * @throws SingularMatrixException if any component's covariance matrix is singular during
 *         fitting.
 * @throws NonPositiveDefiniteMatrixException if any component's covariance matrix is not positive
 *         definite during fitting.
 */
public boolean fit(MixtureMultivariateGaussianDistribution initialMixture, int maxIterations, DoubleDoubleBiPredicate convergencePredicate) {
    ValidationUtils.checkStrictlyPositive(maxIterations, "maxIterations");
    ValidationUtils.checkNotNull(convergencePredicate, "convergencePredicate");
    ValidationUtils.checkNotNull(initialMixture, "initialMixture");
    final int n = data.length;
    final int k = initialMixture.weights.length;
    // Number of data columns.
    final int numCols = data[0].length;
    final int numMeanColumns = initialMixture.distributions[0].means.length;
    ValidationUtils.checkArgument(numCols == numMeanColumns, "Mixture model dimension mismatch with data columns: %d != %d", numCols, numMeanColumns);
    logLikelihood = -Double.MAX_VALUE;
    iterations = 0;
    // Initialize model to fit to initial mixture.
    fittedModel = initialMixture;
    while (iterations++ <= maxIterations) {
        final double previousLogLikelihood = logLikelihood;
        double sumLogLikelihood = 0;
        // Weight and distribution of each component
        final double[] weights = fittedModel.weights;
        final MultivariateGaussianDistribution[] mvns = fittedModel.distributions;
        // E-step: compute the data dependent parameters of the expectation function.
        // The percentage of row's total density between a row and a component
        final double[][] gamma = new double[n][k];
        // Sum of gamma for each component
        final double[] gammaSums = new double[k];
        // Sum of gamma times its row for each each component
        final double[][] gammaDataProdSums = new double[k][numCols];
        // Cache for the weight multiplied by the distribution density
        final double[] mvnDensity = new double[k];
        for (int i = 0; i < n; i++) {
            final double[] point = data[i];
            // Compute densities for each component and the row density
            double rowDensity = 0;
            for (int j = 0; j < k; j++) {
                final double d = weights[j] * mvns[j].density(point);
                mvnDensity[j] = d;
                rowDensity += d;
            }
            sumLogLikelihood += Math.log(rowDensity);
            for (int j = 0; j < k; j++) {
                gamma[i][j] = mvnDensity[j] / rowDensity;
                gammaSums[j] += gamma[i][j];
                for (int col = 0; col < numCols; col++) {
                    gammaDataProdSums[j][col] += gamma[i][j] * point[col];
                }
            }
        }
        logLikelihood = sumLogLikelihood;
        // M-step: compute the new parameters based on the expectation function.
        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];
            }
        }
        // Compute new covariance matrices.
        // These are symmetric so we compute the triangular half.
        final double[][][] newCovMats = new double[k][numCols][numCols];
        final double[] vec = new double[numCols];
        for (int i = 0; i < n; i++) {
            final double[] point = data[i];
            for (int j = 0; j < k; j++) {
                subtract(point, newMeans[j], vec);
                final double g = gamma[i][j];
                // covariance = vec * vecT
                // covariance[ii][jj] = vec[ii] * vec[jj] * gamma[i][j]
                final double[][] covar = newCovMats[j];
                for (int ii = 0; ii < numCols; ii++) {
                    // pre-compute
                    final double vig = vec[ii] * g;
                    final double[] covari = covar[ii];
                    for (int jj = 0; jj <= ii; jj++) {
                        covari[jj] += vig * vec[jj];
                    }
                }
            }
        }
        // Converting to arrays for use by fitted model
        final MultivariateGaussianDistribution[] distributions = new MultivariateGaussianDistribution[k];
        for (int j = 0; j < k; j++) {
            // Make symmetric and normalise by gamma sum
            final double norm = 1.0 / gammaSums[j];
            final double[][] covar = newCovMats[j];
            for (int ii = 0; ii < numCols; ii++) {
                // diagonal
                covar[ii][ii] *= norm;
                // elements
                for (int jj = 0; jj < ii; jj++) {
                    final double tmp = covar[ii][jj] * norm;
                    covar[ii][jj] = tmp;
                    covar[jj][ii] = tmp;
                }
            }
            distributions[j] = new MultivariateGaussianDistribution(newMeans[j], covar);
        }
        // Update current model
        fittedModel = MixtureMultivariateGaussianDistribution.create(newWeights, distributions);
        // Check convergence
        if (convergencePredicate.test(previousLogLikelihood, logLikelihood)) {
            return true;
        }
    }
    // No convergence
    return false;
}

Example 5 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class MultivariateGaussianMixtureExpectationMaximizationTest method canCreateUnmixedMultivariateGaussianDistribution.

@Test
void canCreateUnmixedMultivariateGaussianDistribution() {
    final double[][] data = { { 1, 2 }, { 2.5, 1.5 }, { 3.5, 1.0 } };
    final double[] means = getColumnMeans(data);
    final double[][] covariances = getCovariance(data);
    final MultivariateGaussianDistribution exp = MultivariateGaussianDistribution.create(means, covariances);
    final MultivariateGaussianDistribution obs = MultivariateGaussianMixtureExpectationMaximization.createUnmixed(data);
    final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-6);
    TestAssertions.assertArrayTest(exp.getMeans(), obs.getMeans(), test);
    TestAssertions.assertArrayTest(exp.getCovariances(), obs.getCovariances(), test);
}