Examples with Mean uk.ac.sussex.gdsc.core.math.Mean used on opensource projects

Example 1 with Mean

use of uk.ac.sussex.gdsc.core.math.Mean in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method getTracks.

/**
 * Gets the tracks. Each track has contiguous frames and the length is enough to fit
 * {@code minTrackLength} overlapping windows of the specified size:
 *
 * <pre>
 * length >= window + minTrackLength - 1
 * </pre>
 *
 * @param combinedResults the combined results
 * @param window the window size
 * @param minTrackLength the minimum track length (assumed to be {@code >= 1})
 * @return the tracks
 */
private static List<Trace> getTracks(List<MemoryPeakResults> combinedResults, int window, int minTrackLength) {
    final LocalList<Trace> tracks = new LocalList<>();
    final Statistics stats = new Statistics();
    final int minSize = window + Math.max(minTrackLength, 1) - 1;
    combinedResults.forEach(results -> {
        final int start = tracks.size();
        // Sort by id then frame
        results = results.copy();
        results.sort(IdFramePeakResultComparator.INSTANCE);
        final int size = results.size();
        // Skip IDs not associated with clustering
        int index = 0;
        while (index < size && results.get(index).getId() < 1) {
            index++;
        }
        // Initialise current id and frame
        int id = results.get(index).getId() - 1;
        int frame = results.get(index).getFrame();
        Trace track = new Trace();
        for (; index < size; index++) {
            final PeakResult result = results.get(index);
            // Same ID and contiguous frames
            if (result.getId() != id || result.getFrame() != frame + 1) {
                addTrack(minSize, tracks, track);
                track = new Trace();
            }
            id = result.getId();
            frame = result.getFrame();
            track.add(result);
        }
        addTrack(minSize, tracks, track);
        stats.reset();
        for (int i = start; i < tracks.size(); i++) {
            stats.add(tracks.unsafeGet(i).size());
        }
        final StringBuilder sb = new StringBuilder(256);
        TextUtils.formatTo(sb, "%s tracks=%d, length=%s +/- %s", results.getName(), stats.getN(), MathUtils.rounded(stats.getMean(), 3), MathUtils.rounded(stats.getStandardDeviation(), 3));
        // Limit of diffusion coefficient from the localisation precision.
        // Just use the entire dataset for simplicity (i.e. not the tracks of min length).
        final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
        try {
            pp.getPrecision();
            final Mean mean = new Mean();
            for (final double p : pp.precisions) {
                mean.add(p);
            }
            // 2nDt = MSD (n = number of dimensions)
            // D = MSD / 2nt
            final CalibrationReader reader = results.getCalibrationReader();
            final double t = reader.getExposureTime() / 1000.0;
            // Assume computed in nm. Convert to um.
            final double x = mean.getMean() / 1000;
            final double d = x * x / (2 * t);
            TextUtils.formatTo(sb, ", precision=%s nm, D limit=%s um^2/s", MathUtils.rounded(x * 1000, 4), MathUtils.rounded(d, 4));
        } catch (final DataException ex) {
        // No precision
        }
        IJ.log(sb.toString());
    });
    return tracks;
}

Example 2 with Mean

use of uk.ac.sussex.gdsc.core.math.Mean in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method fitGaussianMixture.

/**
 * Fit the Gaussian mixture to the data. The fitter with the highest likelihood from a number of
 * repeats is returned.
 *
 * @param data the data
 * @param sortDimension the sort dimension
 * @return the multivariate gaussian mixture
 */
private MultivariateGaussianMixtureExpectationMaximization fitGaussianMixture(final double[][] data, int sortDimension) {
    // Get the unmixed multivariate Guassian.
    MultivariateGaussianDistribution unmixed = MultivariateGaussianMixtureExpectationMaximization.createUnmixed(data);
    // Normalise the columns of the data
    // Get means and SD of each column
    final double[] means = unmixed.getMeans();
    final double[] sd = unmixed.getStandardDeviations();
    final int dimensions = means.length;
    for (final double[] value : data) {
        for (int i = 0; i < dimensions; i++) {
            value[i] = (value[i] - means[i]) / sd[i];
        }
    }
    // Repeat. The mean should be approximately 0 and std.dev. 1.
    unmixed = MultivariateGaussianMixtureExpectationMaximization.createUnmixed(data);
    // Record the likelihood of the unmixed model
    double logLikelihood = Arrays.stream(data).mapToDouble(unmixed::density).map(Math::log).sum();
    // x means, x*x covariances
    final int parametersPerGaussian = dimensions + dimensions * dimensions;
    double aic = MathUtils.getAkaikeInformationCriterion(logLikelihood, parametersPerGaussian);
    double bic = MathUtils.getBayesianInformationCriterion(logLikelihood, data.length, parametersPerGaussian);
    ImageJUtils.log("1 component log-likelihood=%s. AIC=%s. BIC=%s", logLikelihood, aic, bic);
    // Fit a mixed component model.
    // Increment the number of components up to a maximim or when the model does not improve.
    MultivariateGaussianMixtureExpectationMaximization mixed = null;
    for (int numComponents = 2; numComponents <= settings.maxComponents; numComponents++) {
        final MultivariateGaussianMixtureExpectationMaximization mixed2 = createMixed(data, dimensions, numComponents);
        if (mixed2 == null) {
            ImageJUtils.log("Failed to fit a %d component mixture model", numComponents);
            break;
        }
        final double logLikelihood2 = mixed2.getLogLikelihood();
        // n * (means, covariances, 1 weight) - 1
        // (Note: subtract 1 as the weights are constrained by summing to 1)
        final int param2 = numComponents * (parametersPerGaussian + 1) - 1;
        final double aic2 = MathUtils.getAkaikeInformationCriterion(logLikelihood2, param2);
        final double bic2 = MathUtils.getBayesianInformationCriterion(logLikelihood2, data.length, param2);
        // Log-likelihood ratio test statistic
        final double lambdaLr = -2 * (logLikelihood - logLikelihood2);
        // DF = difference in dimensionality from previous number of components
        // means, covariances, 1 weight
        final int degreesOfFreedom = parametersPerGaussian + 1;
        final double q = ChiSquaredDistributionTable.computeQValue(lambdaLr, degreesOfFreedom);
        ImageJUtils.log("%d component log-likelihood=%s. AIC=%s. BIC=%s. LLR significance=%s.", numComponents, logLikelihood2, aic2, bic2, MathUtils.rounded(q));
        final double[] weights = mixed2.getFittedModel().getWeights();
        // For consistency sort the mixture by the mean of the diffusion coefficient
        final double[] values = Arrays.stream(mixed2.getFittedModel().getDistributions()).mapToDouble(d -> d.getMeans()[sortDimension]).toArray();
        SortUtils.sortData(weights, values, false, false);
        ImageJUtils.log("Population weights: " + Arrays.toString(weights));
        if (MathUtils.min(weights) < settings.minWeight) {
            ImageJUtils.log("%d component model has population weight %s under minimum level %s", numComponents, MathUtils.min(weights), settings.minWeight);
            break;
        }
        if (aic <= aic2 || bic <= bic2 || q > 0.001) {
            ImageJUtils.log("%d component model is not significant", numComponents);
            break;
        }
        aic = aic2;
        bic = bic2;
        logLikelihood = logLikelihood2;
        mixed = mixed2;
    }
    return mixed;
}