Examples with Gaussian org.apache.commons.math3.analysis.function.Gaussian used on opensource projects

Example 46 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project metron by apache.

the class StatisticalBinningPerformanceDriver method main.

public static void main(String... argv) {
    DescriptiveStatistics perfStats = new DescriptiveStatistics();
    OnlineStatisticsProvider statsProvider = new OnlineStatisticsProvider();
    List<Double> values = new ArrayList<>();
    GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(0L));
    for (int i = 0; i < NUM_DATA_POINTS; ++i) {
        // get the data point out of the [0,1] range
        double d = 1000 * gaussian.nextNormalizedDouble();
        values.add(d);
        statsProvider.addValue(d);
    }
    for (int perfRun = 0; perfRun < NUM_RUNS; ++perfRun) {
        StellarStatisticsFunctions.StatsBin bin = new StellarStatisticsFunctions.StatsBin();
        long start = System.currentTimeMillis();
        Random r = new Random(0);
        for (int i = 0; i < TRIALS_PER_RUN; ++i) {
            // grab a random value and fuzz it a bit so we make sure there's no cheating via caching in t-digest.
            bin.apply(ImmutableList.of(statsProvider, values.get(r.nextInt(values.size())) - 3.5, PERCENTILES));
        }
        perfStats.addValue(System.currentTimeMillis() - start);
    }
    System.out.println("Min/25th/50th/75th/Max Milliseconds: " + perfStats.getMin() + " / " + perfStats.getPercentile(25) + " / " + perfStats.getPercentile(50) + " / " + perfStats.getPercentile(75) + " / " + perfStats.getMax());
}

Example 47 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project pyramid by cheng-li.

the class MultiLabelSynthesizer method gaussianNoise.

/**
 * 2 labels, 3 features, multi-variate gaussian noise
 * y0: w=(0,1,0)
 * y1: w=(1,0,0)
 * y2: w=(0,0,1)
 * @return
 */
public static MultiLabelClfDataSet gaussianNoise(int numData) {
    int numClass = 3;
    int numFeature = 3;
    MultiLabelClfDataSet dataSet = MLClfDataSetBuilder.getBuilder().numFeatures(numFeature).numClasses(numClass).numDataPoints(numData).build();
    // generate weights
    Vector[] weights = new Vector[numClass];
    for (int k = 0; k < numClass; k++) {
        Vector vector = new DenseVector(numFeature);
        weights[k] = vector;
    }
    weights[0].set(1, 1);
    weights[1].set(0, 1);
    weights[2].set(2, 1);
    // generate features
    for (int i = 0; i < numData; i++) {
        for (int j = 0; j < numFeature; j++) {
            dataSet.setFeatureValue(i, j, Sampling.doubleUniform(-1, 1));
        }
    }
    double[] means = new double[numClass];
    double[][] covars = new double[numClass][numClass];
    covars[0][0] = 0.5;
    covars[0][1] = 0.02;
    covars[1][0] = 0.02;
    covars[0][2] = -0.03;
    covars[2][0] = -0.03;
    covars[1][1] = 0.2;
    covars[1][2] = -0.03;
    covars[2][1] = -0.03;
    covars[2][2] = 0.3;
    MultivariateNormalDistribution distribution = new MultivariateNormalDistribution(means, covars);
    // assign labels
    int numFlipped = 0;
    for (int i = 0; i < numData; i++) {
        double[] noises = distribution.sample();
        for (int k = 0; k < numClass; k++) {
            double dot = weights[k].dot(dataSet.getRow(i));
            double score = dot + noises[k];
            if (score >= 0) {
                dataSet.addLabel(i, k);
            }
            if (dot * score < 0) {
                numFlipped += 1;
            }
        }
    }
    System.out.println("number of flipped bits = " + numFlipped);
    return dataSet;
}

Example 48 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian 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;
}

Example 49 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class SpotInspector method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        IJ.showStatus("");
        return;
    }
    // Check if the original image is open
    ImageSource source = results.getSource();
    if (source == null) {
        IJ.error(TITLE, "Unknown original source image");
        return;
    }
    source = source.getOriginal();
    source.setReadHint(ReadHint.NONSEQUENTIAL);
    if (!source.open()) {
        IJ.error(TITLE, "Cannot open original source image: " + source.toString());
        return;
    }
    final float stdDevMax = getStandardDeviation(results);
    if (stdDevMax < 0) {
        // TODO - Add dialog to get the initial peak width
        IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
        return;
    }
    // Rank spots
    rankedResults = new ArrayList<>(results.size());
    // Data for the sorting
    final PrecisionResultProcedure pp;
    if (settings.sortOrderIndex == 1) {
        pp = new PrecisionResultProcedure(results);
        pp.getPrecision();
    } else {
        pp = null;
    }
    // Build procedures to get:
    // Shift = position in pixels - originXY
    final StandardResultProcedure sp;
    if (settings.sortOrderIndex == 9) {
        sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
        sp.getXyr();
    } else {
        sp = null;
    }
    // SD = gaussian widths only for Gaussian PSFs
    final WidthResultProcedure wp;
    if (settings.sortOrderIndex >= 6 && settings.sortOrderIndex <= 8) {
        wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
        wp.getWxWy();
    } else {
        wp = null;
    }
    // Amplitude for Gaussian PSFs
    final HeightResultProcedure hp;
    if (settings.sortOrderIndex == 2) {
        hp = new HeightResultProcedure(results, IntensityUnit.PHOTON);
        hp.getH();
    } else {
        hp = null;
    }
    final Counter c = new Counter();
    results.forEach((PeakResultProcedure) result -> {
        final float[] score = getScore(result, c.getAndIncrement(), pp, sp, wp, hp, stdDevMax);
        rankedResults.add(new PeakResultRank(result, score[0], score[1]));
    });
    Collections.sort(rankedResults, PeakResultRank::compare);
    // Prepare results table
    final ImageJTablePeakResults table = new ImageJTablePeakResults(false, results.getName(), true);
    table.copySettings(results);
    table.setTableTitle(TITLE);
    table.setAddCounter(true);
    table.setShowZ(results.is3D());
    // TODO - Add to settings
    table.setShowFittingData(true);
    table.setShowNoiseData(true);
    if (settings.showCalibratedValues) {
        table.setDistanceUnit(DistanceUnit.NM);
        table.setIntensityUnit(IntensityUnit.PHOTON);
    } else {
        table.setDistanceUnit(DistanceUnit.PIXEL);
        table.setIntensityUnit(IntensityUnit.COUNT);
    }
    table.begin();
    // Add a mouse listener to jump to the frame for the clicked line
    textPanel = table.getResultsWindow().getTextPanel();
    // We must ignore old instances of this class from the mouse listeners
    id = currentId.incrementAndGet();
    textPanel.addMouseListener(new MouseAdapter() {

        @Override
        public void mouseClicked(MouseEvent event) {
            SpotInspector.this.mouseClicked(event);
        }
    });
    // Add results to the table
    int count = 0;
    for (final PeakResultRank rank : rankedResults) {
        rank.rank = count++;
        table.add(rank.peakResult);
    }
    table.end();
    if (settings.plotScore || settings.plotHistogram) {
        // Get values for the plots
        float[] xValues = null;
        float[] yValues = null;
        double yMin;
        double yMax;
        int spotNumber = 0;
        xValues = new float[rankedResults.size()];
        yValues = new float[xValues.length];
        for (final PeakResultRank rank : rankedResults) {
            xValues[spotNumber] = spotNumber + 1;
            yValues[spotNumber++] = recoverScore(rank.score);
        }
        // Set the min and max y-values using 1.5 x IQR
        final DescriptiveStatistics stats = new DescriptiveStatistics();
        for (final float v : yValues) {
            stats.addValue(v);
        }
        if (settings.removeOutliers) {
            final double lower = stats.getPercentile(25);
            final double upper = stats.getPercentile(75);
            final double iqr = upper - lower;
            yMin = Math.max(lower - iqr, stats.getMin());
            yMax = Math.min(upper + iqr, stats.getMax());
            IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
        } else {
            yMin = stats.getMin();
            yMax = stats.getMax();
            IJ.log(String.format("Data range: %f - %f", yMin, yMax));
        }
        plotScore(xValues, yValues, yMin, yMax);
        plotHistogram(yValues);
    }
    // Extract spots into a stack
    final int w = source.getWidth();
    final int h = source.getHeight();
    final int size = 2 * settings.radius + 1;
    final ImageStack spots = new ImageStack(size, size, rankedResults.size());
    // To assist the extraction of data from the image source, process them in time order to allow
    // frame caching. Then set the appropriate slice in the result stack
    Collections.sort(rankedResults, (o1, o2) -> Integer.compare(o1.peakResult.getFrame(), o2.peakResult.getFrame()));
    for (final PeakResultRank rank : rankedResults) {
        final PeakResult r = rank.peakResult;
        // Extract image
        // Note that the coordinates are relative to the middle of the pixel (0.5 offset)
        // so do not round but simply convert to int
        final int x = (int) (r.getXPosition());
        final int y = (int) (r.getYPosition());
        // Extract a region but crop to the image bounds
        int minX = x - settings.radius;
        int minY = y - settings.radius;
        final int maxX = Math.min(x + settings.radius + 1, w);
        final int maxY = Math.min(y + settings.radius + 1, h);
        int padX = 0;
        int padY = 0;
        if (minX < 0) {
            padX = -minX;
            minX = 0;
        }
        if (minY < 0) {
            padY = -minY;
            minY = 0;
        }
        final int sizeX = maxX - minX;
        final int sizeY = maxY - minY;
        float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
        // Prevent errors with missing data
        if (data == null) {
            data = new float[sizeX * sizeY];
        }
        ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
        // Pad if necessary, i.e. the crop is too small for the stack
        if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
            final ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
            spotIp2.insert(spotIp, padX, padY);
            spotIp = spotIp2;
        }
        final int slice = rank.rank + 1;
        spots.setPixels(spotIp.getPixels(), slice);
        spots.setSliceLabel(MathUtils.rounded(rank.originalScore), slice);
    }
    source.close();
    // Reset to the rank order
    Collections.sort(rankedResults, PeakResultRank::compare);
    final ImagePlus imp = ImageJUtils.display(TITLE, spots);
    imp.setRoi((PointRoi) null);
    // Make bigger
    for (int i = 10; i-- > 0; ) {
        imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
    }
}

Example 50 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    settings = Settings.load();
    // Saved by reference so just save now
    settings.save();
    // Read in multiple traced datasets
    // All datasets must have the same pixel pitch and exposure time
    // Get parameters
    // Convert datasets to tracks
    // For each track compute the 4 local track features using the configured window
    // 
    // Optional:
    // Fit a multi-variate Gaussian mixture model to the data
    // (using the configured number of components/populations)
    // Assign each point in the track using the model.
    // Smooth the assignments.
    // 
    // The alternative is to use the localisation category to assign populations.
    // 
    // Plot histograms of each track parameter, coloured by component
    final List<MemoryPeakResults> combinedResults = new LocalList<>();
    if (!showInputDialog(combinedResults)) {
        return;
    }
    final boolean hasCategory = showHasCategoryDialog(combinedResults);
    if (!showDialog(hasCategory)) {
        return;
    }
    ImageJUtils.log(TITLE + "...");
    final List<Trace> tracks = getTracks(combinedResults, settings.window, settings.minTrackLength);
    if (tracks.isEmpty()) {
        IJ.error(TITLE, "No tracks. Please check the input data and min track length setting.");
        return;
    }
    final Calibration cal = combinedResults.get(0).getCalibration();
    final CalibrationReader cr = new CalibrationReader(cal);
    // Use micrometer / second
    final TypeConverter<DistanceUnit> distanceConverter = cr.getDistanceConverter(DistanceUnit.UM);
    final double exposureTime = cr.getExposureTime() / 1000.0;
    final Pair<int[], double[][]> trackData = extractTrackData(tracks, distanceConverter, exposureTime, hasCategory);
    final double[][] data = trackData.getValue();
    // Histogram the raw data.
    final Array2DRowRealMatrix raw = new Array2DRowRealMatrix(data, false);
    final WindowOrganiser wo = new WindowOrganiser();
    // Store the histogram data for plotting the components
    final double[][] columns = new double[FEATURE_NAMES.length][];
    final double[][] limits = new double[FEATURE_NAMES.length][];
    // Get column data
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        columns[i] = raw.getColumn(i);
        if (i == FEATURE_D) {
            // Plot using a logarithmic scale
            SimpleArrayUtils.apply(columns[i], Math::log10);
        }
        limits[i] = MathUtils.limits(columns[i]);
    }
    // Compute histogram bins
    final int[] bins = new int[FEATURE_NAMES.length];
    if (settings.histogramBins > 0) {
        Arrays.fill(bins, settings.histogramBins);
    } else {
        for (int i = 0; i < FEATURE_NAMES.length; i++) {
            bins[i] = HistogramPlot.getBins(StoredData.create(columns[i]), BinMethod.FD);
        }
        // Use the maximum so all histograms look the same
        Arrays.fill(bins, MathUtils.max(bins));
    }
    // Compute plots
    final Plot[] plots = new Plot[FEATURE_NAMES.length];
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        final double[][] hist = HistogramPlot.calcHistogram(columns[i], limits[i][0], limits[i][1], bins[i]);
        plots[i] = new Plot(TITLE + " " + FEATURE_NAMES[i], getFeatureLabel(i, i == FEATURE_D), "Frequency");
        plots[i].addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(plots[i].getTitle(), plots[i], ImageJUtils.NO_TO_FRONT, wo);
    }
    wo.tile();
    // The component for each data point
    int[] component;
    // The number of components
    int numComponents;
    // Data used to fit the Gaussian mixture model
    double[][] fitData;
    // The fitted model
    MixtureMultivariateGaussianDistribution model;
    if (hasCategory) {
        // Use the category as the component.
        // No fit data and no output model
        fitData = null;
        model = null;
        // The component is stored at the end of the raw track data.
        final int end = data[0].length - 1;
        component = Arrays.stream(data).mapToInt(d -> (int) d[end]).toArray();
        numComponents = MathUtils.max(component) + 1;
        // In the EM algorithm the probability of each data point is computed and normalised to
        // sum to 1. The normalised probabilities are averaged to create the weights.
        // Note the probability of each data point uses the previous weight and the algorithm
        // iterates.
        // This is not a fitted model but the input model so use
        // zero weights to indicate no fitting was performed.
        final double[] weights = new double[numComponents];
        // Remove the trailing component to show the 'model' in a table.
        createModelTable(Arrays.stream(data).map(d -> Arrays.copyOf(d, end)).toArray(double[][]::new), weights, component);
    } else {
        // Multivariate Gaussian mixture EM
        // Provide option to not use the anomalous exponent in the population mix.
        int sortDimension = SORT_DIMENSION;
        if (settings.ignoreAlpha) {
            // Remove index 0. This shifts the sort dimension.
            sortDimension--;
            fitData = Arrays.stream(data).map(d -> Arrays.copyOfRange(d, 1, d.length)).toArray(double[][]::new);
        } else {
            fitData = SimpleArrayUtils.deepCopy(data);
        }
        final MultivariateGaussianMixtureExpectationMaximization mixed = fitGaussianMixture(fitData, sortDimension);
        if (mixed == null) {
            IJ.error(TITLE, "Failed to fit a mixture model");
            return;
        }
        model = sortComponents(mixed.getFittedModel(), sortDimension);
        // For the best model, assign to the most likely population.
        component = assignData(fitData, model);
        // Table of the final model using the original data (i.e. not normalised)
        final double[] weights = model.getWeights();
        numComponents = weights.length;
        createModelTable(data, weights, component);
    }
    // Output coloured histograms of the populations.
    final LUT lut = LutHelper.createLut(settings.lutIndex);
    IntFunction<Color> colourMap;
    if (LutHelper.getColour(lut, 0).equals(Color.BLACK)) {
        colourMap = i -> LutHelper.getNonZeroColour(lut, i, 0, numComponents - 1);
    } else {
        colourMap = i -> LutHelper.getColour(lut, i, 0, numComponents - 1);
    }
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        // Extract the data for each component
        final double[] col = columns[i];
        final Plot plot = plots[i];
        for (int n = 0; n < numComponents; n++) {
            final StoredData feature = new StoredData();
            for (int j = 0; j < component.length; j++) {
                if (component[j] == n) {
                    feature.add(col[j]);
                }
            }
            if (feature.size() == 0) {
                continue;
            }
            final double[][] hist = HistogramPlot.calcHistogram(feature.values(), limits[i][0], limits[i][1], bins[i]);
            // Colour the points
            plot.setColor(colourMap.apply(n));
            plot.addPoints(hist[0], hist[1], Plot.BAR);
        }
        plot.updateImage();
    }
    createTrackDataTable(tracks, trackData, fitData, model, component, cal, colourMap);
// Analysis.
// Assign the original localisations to their track component.
// Q. What about the start/end not covered by the window?
// Save tracks as a dataset labelled with the sub-track ID?
// Output for the bound component and free components track parameters.
// Compute dwell times.
// Other ...
// Track analysis plugin:
// Extract all continuous segments of the same component.
// Produce MSD plot with error bars.
// Fit using FBM model.
}