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

Example 11 with Gaussian

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

the class OnlineStatisticsProviderTest method testNormallyDistributedRandomDataShifted.

@Test
public void testNormallyDistributedRandomDataShifted() {
    List<Double> values = new ArrayList<>();
    GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(0L));
    for (int i = 0; i < 1000000; ++i) {
        double d = gaussian.nextNormalizedDouble() + 10;
        values.add(d);
    }
    validateEquality(values);
}

Example 12 with Gaussian

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

the class StellarStatisticsFunctionsTest method testMergeProviders.

@Test
public void testMergeProviders() throws Exception {
    List<StatisticsProvider> providers = new ArrayList<>();
    /*
    Create 10 providers, each with a sample drawn from a gaussian distribution.
    Update the reference stats from commons math to ensure we are
     */
    GaussianRandomGenerator gaussian = new GaussianRandomGenerator(new MersenneTwister(1L));
    SummaryStatistics sStatistics = new SummaryStatistics();
    DescriptiveStatistics dStatistics = new DescriptiveStatistics();
    for (int i = 0; i < 10; ++i) {
        List<Double> sample = new ArrayList<>();
        for (int j = 0; j < 100; ++j) {
            double s = gaussian.nextNormalizedDouble();
            sample.add(s);
            sStatistics.addValue(s);
            dStatistics.addValue(s);
        }
        StatisticsProvider provider = (StatisticsProvider) run("STATS_ADD(STATS_INIT(), " + Joiner.on(",").join(sample) + ")", new HashMap<>());
        providers.add(provider);
    }
    /*
    Merge the providers and validate
     */
    Map<String, Object> providerVariables = new HashMap<>();
    for (int i = 0; i < providers.size(); ++i) {
        providerVariables.put("provider_" + i, providers.get(i));
    }
    StatisticsProvider mergedProvider = (StatisticsProvider) run("STATS_MERGE([" + Joiner.on(",").join(providerVariables.keySet()) + "])", providerVariables);
    OnlineStatisticsProviderTest.validateStatisticsProvider(mergedProvider, sStatistics, dStatistics);
}

Example 13 with Gaussian

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

the class FitTimeSeries method fitToGaussian.

static GaussianFitResults fitToGaussian(double init_center_i, double init_center_j, double init_radius2, FloatImage image) {
    // 
    // do some optimization on the image (fitting to a Gaussian)
    // set initial guesses from ROI operation.
    // 
    ISize imageSize = image.getISize();
    final int num_i = imageSize.getX();
    final int num_j = imageSize.getY();
    final float[] floatPixels = image.getFloatPixels();
    // 
    // initial guess based on previous fit of ROI
    // do gaussian fit in index space for center and standard deviation (later to translate it back to world coordinates)
    // 
    final int window_size = (int) Math.sqrt(init_radius2) * 4;
    // final int window_min_i = 0;       // (int) Math.max(0, Math.floor(init_center_i - window_size/2));
    // final int window_max_i = num_i-1; // (int) Math.min(num_i-1, Math.ceil(init_center_i + window_size/2));
    // final int window_min_j = 0;       // (int) Math.max(0, Math.floor(init_center_j - window_size/2));
    // final int window_max_j = num_j-1; // (int) Math.min(num_j-1, Math.ceil(init_center_j + window_size/2));
    final int window_min_i = (int) Math.max(0, Math.floor(init_center_i - window_size / 2));
    final int window_max_i = (int) Math.min(num_i - 1, Math.ceil(init_center_i + window_size / 2));
    final int window_min_j = (int) Math.max(0, Math.floor(init_center_j - window_size / 2));
    final int window_max_j = (int) Math.min(num_j - 1, Math.ceil(init_center_j + window_size / 2));
    final int PARAM_INDEX_CENTER_I = 0;
    final int PARAM_INDEX_CENTER_J = 1;
    final int PARAM_INDEX_K = 2;
    final int PARAM_INDEX_HIGH = 3;
    final int PARAM_INDEX_RADIUS_SQUARED = 4;
    final int NUM_PARAMETERS = 5;
    double[] initParameters = new double[NUM_PARAMETERS];
    initParameters[PARAM_INDEX_CENTER_I] = init_center_i;
    initParameters[PARAM_INDEX_CENTER_J] = init_center_j;
    initParameters[PARAM_INDEX_HIGH] = 1.0;
    initParameters[PARAM_INDEX_K] = 10;
    initParameters[PARAM_INDEX_RADIUS_SQUARED] = init_radius2;
    PowellOptimizer optimizer = new PowellOptimizer(1e-4, 1e-1);
    MultivariateFunction func = new MultivariateFunction() {

        @Override
        public double value(double[] point) {
            double center_i = point[PARAM_INDEX_CENTER_I];
            double center_j = point[PARAM_INDEX_CENTER_J];
            double high = point[PARAM_INDEX_HIGH];
            double K = point[PARAM_INDEX_K];
            double radius2 = point[PARAM_INDEX_RADIUS_SQUARED];
            double error2 = 0;
            for (int j = window_min_j; j <= window_max_j; j++) {
                // double y = j - center_j;
                double y = j;
                for (int i = window_min_i; i <= window_max_i; i++) {
                    // double x = i - center_i;
                    double x = i;
                    double modelValue = high - FastMath.exp(-K * FastMath.exp(-2 * (x * x + y * y) / radius2));
                    double imageValue = floatPixels[j * num_i + i];
                    double error = modelValue - imageValue;
                    error2 += error * error;
                }
            }
            System.out.println(new GaussianFitResults(center_i, center_j, radius2, K, high, error2));
            return error2;
        }
    };
    PointValuePair pvp = optimizer.optimize(new ObjectiveFunction(func), new InitialGuess(initParameters), new MaxEval(100000), GoalType.MINIMIZE);
    double[] fittedParamValues = pvp.getPoint();
    double fitted_center_i = fittedParamValues[PARAM_INDEX_CENTER_I];
    double fitted_center_j = fittedParamValues[PARAM_INDEX_CENTER_J];
    double fitted_radius2 = fittedParamValues[PARAM_INDEX_RADIUS_SQUARED];
    double fitted_K = fittedParamValues[PARAM_INDEX_K];
    double fitted_high = fittedParamValues[PARAM_INDEX_HIGH];
    double objectiveFunctionValue = pvp.getValue();
    return new GaussianFitResults(fitted_center_i, fitted_center_j, fitted_radius2, fitted_K, fitted_high, objectiveFunctionValue);
}

Example 14 with Gaussian

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

the class Slice method findSupportPoints.

public void findSupportPoints() {
    supportPoints.clear();
    supportMax = 0;
    // Map<Line, Double> totalLineSupport = new HashedMap();
    double distSigma = 0.2;
    Gaussian distanceG = new Gaussian(0, distSigma), angleG = new Gaussian(0, 0.3);
    for (Line l : gis.allLines()) {
        double length = l.length();
        int noPoints = (int) Math.max(2, length / 0.1);
        for (int n = 0; n <= noPoints; n++) {
            SupportPoint sp = new SupportPoint(l, l.fromPPram(n / (double) noPoints));
            // Double.MAX_VALUE;
            sp.support = 0;
            for (Line nearL : slice.getNear(sp.pt, distSigma * 3)) {
                double angle = nearL.absAngle(sp.line);
                // Math.min  (sp.pt.distance(line.start), sp.pt.distance(line.end)) ;
                double dist = nearL.distance(sp.pt, true);
                sp.support += distanceG.value(dist) * angleG.value(angle) * nearL.length();
            }
            supportPoints.add(sp);
        }
    }
    foundLines = null;
    repaint();
}

Example 15 with Gaussian

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