Examples with Frequency org.apache.commons.math3.stat.Frequency used on opensource projects

Example 21 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method performDistanceAnalysis.

private void performDistanceAnalysis(double[][] intraHist, int p99) {
    // We want to know the fraction of distances between molecules at the 99th percentile
    // that are intra- rather than inter-molecule.
    // Do single linkage clustering of closest pair at this distance and count the number of
    // links that are inter and intra.
    // Convert molecules for clustering
    final ArrayList<ClusterPoint> points = new ArrayList<>(settings.molecules.size());
    for (final Molecule m : settings.molecules) {
        // Precision was used to store the molecule ID
        points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons));
    }
    final ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, SimpleImageJTrackProgress.getInstance());
    IJ.showStatus("Clustering to check inter-molecule distances");
    engine.setTrackJoins(true);
    final List<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]);
    IJ.showStatus("");
    if (clusters != null) {
        final double[] intraIdDistances = engine.getIntraIdDistances();
        final double[] interIdDistances = engine.getInterIdDistances();
        final int all = interIdDistances.length + intraIdDistances.length;
        log("  * Fraction of inter-molecule particle linkage @ %s nm = %s %%", MathUtils.rounded(intraHist[0][p99], 4), (all > 0) ? MathUtils.rounded(100.0 * interIdDistances.length / all, 4) : "0");
        // Show a double cumulative histogram plot
        final double[][] intraIdHist = MathUtils.cumulativeHistogram(intraIdDistances, false);
        final double[][] interIdHist = MathUtils.cumulativeHistogram(interIdDistances, false);
        // Plot
        final String title = TITLE + " molecule linkage distance";
        final Plot plot = new Plot(title, "Distance", "Frequency");
        plot.addPoints(intraIdHist[0], intraIdHist[1], Plot.LINE);
        double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0;
        if (interIdHist[1].length > 0) {
            max = Math.max(max, interIdHist[1][interIdHist[1].length - 1]);
        }
        plot.setLimits(0, intraIdHist[0][intraIdHist[0].length - 1], 0, max);
        plot.setColor(Color.blue);
        plot.addPoints(interIdHist[0], interIdHist[1], Plot.LINE);
        plot.setColor(Color.black);
        plot.addLegend("Intra-molecule\nInter-molecule");
        ImageJUtils.display(title, plot);
    } else {
        log("Aborted clustering to check inter-molecule distances");
    }
}

Example 22 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class PcPalmClusters method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (!showDialog()) {
        return;
    }
    PcPalmMolecules.logSpacer();
    ImageJUtils.log(TITLE);
    PcPalmMolecules.logSpacer();
    final long start = System.currentTimeMillis();
    HistogramData histogramData;
    if (fileInput) {
        histogramData = loadHistogram(settings.histogramFile);
    } else {
        histogramData = doClustering();
    }
    if (histogramData == null) {
        return;
    }
    final float[][] hist = histogramData.histogram;
    // Create a histogram of the cluster sizes
    String title = TITLE + " Molecules/cluster";
    final String xTitle = "Molecules/cluster";
    final String yTitle = "Frequency";
    // Plot the histogram
    Plot plot = new Plot(title, xTitle, yTitle);
    plot.addPoints(hist[0], hist[1], Plot.BAR);
    ImageJUtils.display(title, plot);
    final HistogramData noiseData = loadNoiseHistogram(histogramData);
    if (noiseData != null && subtractNoise(histogramData, noiseData)) {
        // Update the histogram
        title += " (noise subtracted)";
        plot = new Plot(title, xTitle, yTitle);
        plot.addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(title, plot);
        // Automatically save
        if (autoSave) {
            final String newFilename = FileUtils.replaceExtension(histogramData.filename, ".noise.tsv");
            if (saveHistogram(histogramData, newFilename)) {
                ImageJUtils.log("Saved noise-subtracted histogram to " + newFilename);
            }
        }
    }
    // Fit the histogram
    final double[] fitParameters = fitBinomial(histogramData);
    if (fitParameters != null) {
        // Add the binomial to the histogram
        final int n = (int) fitParameters[0];
        final double p = fitParameters[1];
        ImageJUtils.log("Optimal fit : N=%d, p=%s", n, MathUtils.rounded(p));
        final BinomialDistribution dist = new BinomialDistribution(n, p);
        // A zero-truncated binomial was fitted.
        // pi is the adjustment factor for the probability density.
        final double pi = 1 / (1 - dist.probability(0));
        if (!fileInput) {
            // Calculate the estimated number of clusters from the observed molecules:
            // Actual = (Observed / p-value) / N
            final double actual = (numberOfMolecules / p) / n;
            ImageJUtils.log("Estimated number of clusters : (%d / %s) / %d = %s", numberOfMolecules, MathUtils.rounded(p), n, MathUtils.rounded(actual));
        }
        final double[] x = new double[n + 2];
        final double[] y = new double[n + 2];
        // Scale the values to match those on the histogram
        final double normalisingFactor = count * pi;
        for (int i = 0; i <= n; i++) {
            x[i] = i;
            y[i] = dist.probability(i) * normalisingFactor;
        }
        x[n + 1] = n + 1;
        y[n + 1] = 0;
        // Redraw the plot since the limits may have changed
        plot.setColor(Color.magenta);
        plot.addPoints(x, y, Plot.LINE);
        plot.addPoints(x, y, Plot.CIRCLE);
        plot.setColor(Color.black);
        plot.updateImage();
        ImageJUtils.display(title, plot);
    }
    final double seconds = (System.currentTimeMillis() - start) / 1000.0;
    final String msg = TITLE + " complete : " + seconds + "s";
    IJ.showStatus(msg);
    ImageJUtils.log(msg);
}

Example 23 with Frequency

use of org.apache.commons.math3.stat.Frequency in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method calculateAveragePrecision.

/**
 * Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision
 * distribution.
 *
 * <p>A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does
 * not fit within 3 SDs of the simple mean then the simple mean is returned.
 *
 * @param molecules the molecules
 * @param title the plot title (null if no plot should be displayed)
 * @param histogramBins the histogram bins
 * @param logFitParameters Record the fit parameters to the ImageJ log
 * @param removeOutliers The distribution is created using all values within 1.5x the
 *        inter-quartile range (IQR) of the data
 * @return The average precision
 */
public double calculateAveragePrecision(List<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
    // Plot histogram of the precision
    final float[] data = new float[molecules.size()];
    final DescriptiveStatistics stats = new DescriptiveStatistics();
    double yMin = Double.NEGATIVE_INFINITY;
    double yMax = 0;
    for (int i = 0; i < data.length; i++) {
        data[i] = (float) molecules.get(i).precision;
        stats.addValue(data[i]);
    }
    // Set the min and max y-values using 1.5 x IQR
    if (removeOutliers) {
        final double lower = stats.getPercentile(25);
        final double upper = stats.getPercentile(75);
        if (Double.isNaN(lower) || Double.isNaN(upper)) {
            if (logFitParameters) {
                ImageJUtils.log("Error computing IQR: %f - %f", lower, upper);
            }
        } else {
            final double iqr = upper - lower;
            yMin = Math.max(lower - iqr, stats.getMin());
            yMax = Math.min(upper + iqr, stats.getMax());
            if (logFitParameters) {
                ImageJUtils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
            }
        }
    }
    if (yMin == Double.NEGATIVE_INFINITY) {
        yMin = stats.getMin();
        yMax = stats.getMax();
        if (logFitParameters) {
            ImageJUtils.log("  Data range: %f - %f", yMin, yMax);
        }
    }
    int bins;
    if (histogramBins <= 0) {
        bins = (int) Math.ceil((stats.getMax() - stats.getMin()) / HistogramPlot.getBinWidthScottsRule(stats.getStandardDeviation(), (int) stats.getN()));
    } else {
        bins = histogramBins;
    }
    final float[][] hist = HistogramPlot.calcHistogram(data, yMin, yMax, bins);
    Plot plot = null;
    if (title != null) {
        plot = new Plot(title, "Precision", "Frequency");
        final float[] xValues = hist[0];
        final float[] yValues = hist[1];
        plot.addPoints(xValues, yValues, Plot.BAR);
        ImageJUtils.display(title, plot);
    }
    // Extract non-zero data
    float[] x = Arrays.copyOf(hist[0], hist[0].length);
    float[] y = Arrays.copyOf(hist[1], hist[1].length);
    int count = 0;
    for (int i = 0; i < y.length; i++) {
        if (y[i] > 0) {
            x[count] = x[i];
            y[count] = y[i];
            count++;
        }
    }
    x = Arrays.copyOf(x, count);
    y = Arrays.copyOf(y, count);
    // Sense check to fitted data. Get mean and SD of histogram
    final double[] stats2 = HistogramPlot.getHistogramStatistics(x, y);
    double mean = stats2[0];
    if (logFitParameters) {
        log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    }
    // Standard Gaussian fit
    final double[] parameters = fitGaussian(x, y);
    if (parameters == null) {
        log("  Failed to fit initial Gaussian");
        return mean;
    }
    double newMean = parameters[1];
    double error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    mean = newMean;
    if (logFitParameters) {
        log("  Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
    }
    final double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
    // Fit to a skewed Gaussian (or appropriate function)
    final double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
    if (skewParameters == null) {
        log("  Failed to fit Skewed Gaussian");
        return mean;
    }
    final SkewNormalFunction sn = new SkewNormalFunction(skewParameters);
    if (logFitParameters) {
        log("  Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance()));
    }
    newMean = sn.getMean();
    error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    // Use original histogram x-axis to maintain all the bins
    if (plot != null) {
        plot.setColor(Color.red);
        addToPlot(plot, hist[0], skewParameters, Plot.LINE);
        plot.setColor(Color.black);
        ImageJUtils.display(title, plot);
    }
    // Return the average precision from the fitted curve
    return newMean;
}