Search in sources :

Example 26 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project ffx by mjschnie.

the class ReflectionList method ordinal.

/**
 * <p>
 * ordinal</p>
 *
 * @param s a double.
 * @return a double.
 */
public double ordinal(double s) {
    double r = (s - minres) / (maxres - minres);
    r = min(r, 0.999) * 1000.0;
    int i = (int) r;
    r -= floor(r);
    return ((1.0 - r) * hist[i] + r * hist[i + 1]);
}
Also used : FastMath.rint(org.apache.commons.math3.util.FastMath.rint)

Example 27 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.

the class CreateData method createUniformDistribution.

/**
	 * Create distribution within an XY border
	 * 
	 * @param border
	 * @return
	 */
private UniformDistribution createUniformDistribution(double border) {
    double depth = (settings.fixedDepth) ? settings.depth / settings.pixelPitch : settings.depth / (2 * settings.pixelPitch);
    // Ensure the focal plane is in the middle of the zDepth
    double[] max = new double[] { settings.size / 2 - border, settings.size / 2 - border, depth };
    double[] min = new double[3];
    for (int i = 0; i < 3; i++) min[i] = -max[i];
    if (settings.fixedDepth)
        min[2] = max[2];
    // Try using different distributions:
    final RandomGenerator rand1 = createRandomGenerator();
    if (settings.distribution.equals(DISTRIBUTION[UNIFORM_HALTON])) {
        return new UniformDistribution(min, max, rand1.nextInt());
    }
    if (settings.distribution.equals(DISTRIBUTION[UNIFORM_SOBOL])) {
        SobolSequenceGenerator rvg = new SobolSequenceGenerator(3);
        rvg.skipTo(rand1.nextInt());
        return new UniformDistribution(min, max, rvg);
    }
    // Create a distribution using random generators for each dimension 
    UniformDistribution distribution = new UniformDistribution(min, max, this);
    return distribution;
}
Also used : UniformDistribution(gdsc.smlm.model.UniformDistribution) SobolSequenceGenerator(org.apache.commons.math3.random.SobolSequenceGenerator) RandomGenerator(org.apache.commons.math3.random.RandomGenerator)

Example 28 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min 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
	 * @param title
	 *            the plot title (null if no plot should be displayed)
	 * @param histogramBins
	 * @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(ArrayList<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
    // Plot histogram of the precision
    float[] data = new float[molecules.size()];
    DescriptiveStatistics stats = new DescriptiveStatistics();
    double yMin = Double.NEGATIVE_INFINITY, 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) {
        double lower = stats.getPercentile(25);
        double upper = stats.getPercentile(75);
        if (Double.isNaN(lower) || Double.isNaN(upper)) {
            if (logFitParameters)
                Utils.log("Error computing IQR: %f - %f", lower, upper);
        } else {
            double iqr = upper - lower;
            yMin = FastMath.max(lower - iqr, stats.getMin());
            yMax = FastMath.min(upper + iqr, stats.getMax());
            if (logFitParameters)
                Utils.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)
            Utils.log("  Data range: %f - %f", yMin, yMax);
    }
    float[][] hist = Utils.calcHistogram(data, yMin, yMax, histogramBins);
    Plot2 plot = null;
    if (title != null) {
        plot = new Plot2(title, "Precision", "Frequency");
        float[] xValues = hist[0];
        float[] yValues = hist[1];
        if (xValues.length > 0) {
            double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
            plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.max(yValues) * 1.05);
        }
        plot.addPoints(xValues, yValues, Plot2.BAR);
        Utils.display(title, plot);
    }
    // Extract non-zero data
    float[] x = Arrays.copyOf(hist[0], hist[0].length);
    float[] y = hist[1];
    int count = 0;
    float dx = (x[1] - x[0]) * 0.5f;
    for (int i = 0; i < y.length; i++) if (y[i] > 0) {
        x[count] = x[i] + dx;
        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
    double[] stats2 = Utils.getHistogramStatistics(x, y);
    double mean = stats2[0];
    if (logFitParameters)
        log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    // Standard Gaussian fit
    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]);
    double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
    // Fit to a skewed Gaussian (or appropriate function)
    double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
    if (skewParameters == null) {
        log("  Failed to fit Skewed Gaussian");
        return mean;
    }
    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) {
        x = hist[0];
        for (int i = 0; i < y.length; i++) x[i] += dx;
        plot.setColor(Color.red);
        addToPlot(plot, x, skewParameters, Plot2.LINE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    }
    // Return the average precision from the fitted curve
    return newMean;
}
Also used : DescriptiveStatistics(org.apache.commons.math3.stat.descriptive.DescriptiveStatistics) Plot2(ij.gui.Plot2) SkewNormalFunction(gdsc.smlm.function.SkewNormalFunction) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Example 29 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.

the class PoissonFunctionTest method cumulativeProbabilityIsOneWithRealAbove4.

private void cumulativeProbabilityIsOneWithRealAbove4(final double gain, final double mu, int min, int max) {
    final double o = mu * gain;
    final PoissonFunction f = new PoissonFunction(1.0 / gain, true);
    double p = 0;
    SimpsonIntegrator in = new SimpsonIntegrator(3, 30);
    try {
        p = in.integrate(Integer.MAX_VALUE, new UnivariateFunction() {

            public double value(double x) {
                return f.likelihood(x, o);
            }
        }, min, max);
        System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertEquals(String.format("g=%f, mu=%f", gain, mu), 1, p, 0.02);
    } catch (TooManyEvaluationsException e) {
        //double inc = max / 20000.0;
        //for (double x = 0; x <= max; x += inc)
        //{
        //	final double pp = f.likelihood(x, o);
        //	//System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, pp);
        //	p += pp;
        //}
        //System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertFalse(e.getMessage(), true);
    }
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) TooManyEvaluationsException(org.apache.commons.math3.exception.TooManyEvaluationsException) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Example 30 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method cumulativeProbabilityIsOneWithRealData.

private void cumulativeProbabilityIsOneWithRealData(final double mu, double min, double max, boolean test) {
    double p = 0;
    UnivariateIntegrator in = new SimpsonIntegrator();
    p = in.integrate(20000, new UnivariateFunction() {

        public double value(double x) {
            double v;
            v = PoissonCalculator.likelihood(mu, x);
            //System.out.printf("x=%f, v=%f\n", x, v);
            return v;
        }
    }, min, max);
    System.out.printf("mu=%f, p=%f\n", mu, p);
    if (test) {
        Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
    }
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) UnivariateIntegrator(org.apache.commons.math3.analysis.integration.UnivariateIntegrator) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Aggregations

ArrayList (java.util.ArrayList)16 List (java.util.List)10 DescriptiveStatistics (org.apache.commons.math3.stat.descriptive.DescriptiveStatistics)8 SummaryStatistics (org.apache.commons.math3.stat.descriptive.SummaryStatistics)7 Map (java.util.Map)6 UnivariateFunction (org.apache.commons.math3.analysis.UnivariateFunction)6 MaxEval (org.apache.commons.math3.optim.MaxEval)6 Collectors (java.util.stream.Collectors)5 ExpressionException (cbit.vcell.parser.ExpressionException)4 Plot2 (ij.gui.Plot2)4 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)4 InitialGuess (org.apache.commons.math3.optim.InitialGuess)4 PointValuePair (org.apache.commons.math3.optim.PointValuePair)4 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)4 VisibleForTesting (com.google.common.annotations.VisibleForTesting)3 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)3 HashMap (java.util.HashMap)3 SimpsonIntegrator (org.apache.commons.math3.analysis.integration.SimpsonIntegrator)3 BrentOptimizer (org.apache.commons.math3.optim.univariate.BrentOptimizer)3 SearchInterval (org.apache.commons.math3.optim.univariate.SearchInterval)3