Examples with Gamma org.apache.commons.math3.special.Gamma used on opensource projects

Example 21 with Gamma

use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method fit.

/**
 * Fit the EM-gain distribution (Gaussian * Gamma).
 *
 * @param histogram The distribution
 */
private void fit(int[] histogram) {
    final int[] limits = limits(histogram);
    final double[] x = getX(limits);
    final double[] y = getY(histogram, limits);
    Plot plot = new Plot(TITLE, "ADU", "Frequency");
    double yMax = MathUtils.max(y);
    plot.setLimits(limits[0], limits[1], 0, yMax);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    ImageJUtils.display(TITLE, plot);
    // Estimate remaining parameters.
    // Assuming a gamma_distribution(shape,scale) then mean = shape * scale
    // scale = gain
    // shape = Photons = mean / gain
    double mean = getMean(histogram) - settings.bias;
    // Note: if the bias is too high then the mean will be negative. Just move the bias.
    while (mean < 0) {
        settings.bias -= 1;
        mean += 1;
    }
    double photons = mean / settings.gain;
    if (settings.settingSimulate) {
        ImageJUtils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.settingBias, MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    }
    ImageJUtils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    final int max = (int) x[x.length - 1];
    double[] pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot.setColor(Color.blue);
    plot.addPoints(x, pg, Plot.LINE);
    ImageJUtils.display(TITLE, plot);
    // Perform a fit
    final CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16);
    final double[] startPoint = new double[] { photons, settings.gain, settings.noise, settings.bias };
    int maxEval = 3000;
    final String[] paramNames = { "Photons", "Gain", "Noise", "Bias" };
    // Set bounds
    final double[] lower = new double[] { 0, 0.5 * settings.gain, 0, settings.bias - settings.noise };
    final double[] upper = new double[] { 2 * photons, 2 * settings.gain, settings.gain, settings.bias + settings.noise };
    // Restart until converged.
    // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data.
    PointValuePair solution = null;
    for (int iter = 0; iter < 3; iter++) {
        IJ.showStatus("Fitting histogram ... Iteration " + iter);
        try {
            // Basic Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef()));
            if (solution == null || optimum.getValue() < solution.getValue()) {
                final double[] point = optimum.getPointRef();
                // Check the bounds
                for (int i = 0; i < point.length; i++) {
                    if (point[i] < lower[i] || point[i] > upper[i]) {
                        throw new ComputationException(String.format("Fit out of of estimated range: %s %f", paramNames[i], point[i]));
                    }
                }
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
        try {
            // Bounded Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef())));
            final double[] point = adapter.unboundedToBounded(optimum.getPointRef());
            optimum = new PointValuePair(point, optimum.getValue());
            if (solution == null || optimum.getValue() < solution.getValue()) {
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Bounded Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
    }
    ImageJUtils.finished();
    if (solution == null) {
        ImageJUtils.log("Failed to fit the distribution");
        return;
    }
    final double[] point = solution.getPointRef();
    photons = point[0];
    settings.gain = point[1];
    settings.noise = point[2];
    settings.bias = (int) Math.round(point[3]);
    final String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    ImageJUtils.log(label);
    if (settings.settingSimulate) {
        ImageJUtils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", MathUtils.rounded(relativeError(settings.bias, settings.settingBias)), MathUtils.rounded(relativeError(settings.gain, settings.settingGain)), MathUtils.rounded(relativeError(settings.noise, settings.settingNoise)), MathUtils.rounded(relativeError(photons, settings.settingPhotons)));
    }
    // Show the PoissonGammaGaussian approximation
    double[] approxValues = null;
    if (settings.showApproximation) {
        approxValues = new double[x.length];
        final PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / settings.gain, settings.noise);
        final double expected = photons * settings.gain;
        for (int i = 0; i < approxValues.length; i++) {
            approxValues[i] = fun.likelihood(x[i] - settings.bias, expected);
        }
        yMax = MathUtils.maxDefault(yMax, approxValues);
    }
    // Replot
    pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot = new Plot(TITLE, "ADU", "Frequency");
    plot.setLimits(limits[0], limits[1], 0, yMax * 1.05);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    plot.setColor(Color.red);
    plot.addPoints(x, pg, Plot.LINE);
    plot.addLabel(0, 0, label);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, approxValues, Plot.LINE);
    }
    ImageJUtils.display(TITLE, plot);
}

Example 22 with Gamma

use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method plotJumpDistances.

/**
 * Plot a cumulative histogram and standard histogram of the jump distances.
 *
 * @param title the title
 * @param jumpDistances the jump distances
 * @param dimensions the number of dimensions for the jumps
 */
private void plotJumpDistances(String title, DoubleData jumpDistances, int dimensions) {
    // Cumulative histogram
    // --------------------
    final double[] values = jumpDistances.values();
    String title2 = title + " Cumulative Jump Distance " + dimensions + "D";
    final double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(values);
    Plot jdPlot = new Plot(title2, "Distance (um^2)", "Cumulative Probability");
    jdPlot.addPoints(jdHistogram[0], jdHistogram[1], Plot.LINE);
    ImageJUtils.display(title2, jdPlot, windowOrganiser);
    // Plot the expected function
    // This is the Chi-squared distribution: The sum of the squares of k independent
    // standard normal random variables with k = dimensions. It is a special case of
    // the gamma distribution. If the normals have non-unit variance the distribution
    // is scaled.
    // Chi ~ Gamma(k/2, 2) // using the scale parameterisation of the gamma
    // s^2 * Chi ~ Gamma(k/2, 2*s^2)
    // So if s^2 = 2D:
    // 2D * Chi ~ Gamma(k/2, 4D)
    final double estimatedD = pluginSettings.simpleD * pluginSettings.simpleSteps;
    final double max = MathUtils.max(values);
    final double[] x = SimpleArrayUtils.newArray(1000, 0, max / 1000);
    final double k = dimensions / 2.0;
    final double mean = 4 * estimatedD;
    final GammaDistribution dist = new GammaDistribution(null, k, mean);
    final double[] y = new double[x.length];
    for (int i = 0; i < x.length; i++) {
        y[i] = dist.cumulativeProbability(x[i]);
    }
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(title2, jdPlot);
    // Histogram
    // ---------
    title2 = title + " Jump " + dimensions + "D";
    final HistogramPlot histogramPlot = new HistogramPlotBuilder(title2, jumpDistances, "Distance (um^2)").build();
    // Assume the plot works
    histogramPlot.show(windowOrganiser);
    // Recompute the expected function
    for (int i = 0; i < x.length; i++) {
        y[i] = dist.density(x[i]);
    }
    // Scale to have the same area
    final double[] xvalues = histogramPlot.getPlotXValues();
    if (xvalues.length > 1) {
        final double area1 = jumpDistances.size() * (xvalues[1] - xvalues[0]);
        final double area2 = dist.cumulativeProbability(x[x.length - 1]);
        final double scale = area1 / area2;
        for (int i = 0; i < y.length; i++) {
            y[i] *= scale;
        }
    }
    jdPlot = histogramPlot.getPlot();
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(histogramPlot.getPlotTitle(), jdPlot);
}

Example 23 with Gamma

use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFisherInformation method findUpperLimit.

/**
 * Find the upper limit of the integrated function A^2/P. P is the Poisson-Gamma convolution, A is
 * the partial gradient.
 *
 * <p>When both A and P are convolved with a Gaussian kernel, the integral of this function - 1 is
 * the Fisher information.
 *
 * <p>This method is used to determine the upper limit of the function using a binary search.
 *
 * @param theta the Poisson mean
 * @param max the max of the function (returned from {@link #findMaximum(double, double)})
 * @param rel Relative threshold
 * @return [upper,upper value,evaluations]
 */
public double[] findUpperLimit(final double theta, double[] max, double rel) {
    if (rel < MIN_RELATIVE_TOLERANCE) {
        throw new IllegalArgumentException("Relative tolerance too small: " + rel);
    }
    final UnivariateFunction f = new UnivariateFunction() {

        double[] dgDp = new double[1];

        @Override
        public double value(double x) {
            final double G = PoissonGammaFunction.unscaledPoissonGammaPartial(x, theta, gain, dgDp);
            return getF(G, dgDp[0]);
        }
    };
    int eval = 0;
    // Increase from the max until the tolerance is achieved.
    // Use the mean to get a rough initial step size
    final double mean = theta * gain;
    double step = Math.max(mean, 1);
    double upper = max[0];
    double upperValue = max[1];
    final double threshold = upperValue * rel;
    double lower = upper;
    while (upperValue > threshold) {
        lower = upper;
        upper += step;
        step *= 2;
        eval++;
        upperValue = f.value(upper);
    }
    // Binary search the bracket between lower and upper
    while (lower + 1 < upper) {
        final double mid = (lower + upper) * 0.5;
        eval++;
        final double midg = f.value(mid);
        if (midg > threshold) {
            lower = mid;
        } else {
            upper = mid;
            upperValue = midg;
        }
    }
    return new double[] { upper, upperValue, eval };
}

Example 24 with Gamma

use of org.apache.commons.math3.special.Gamma in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFunction method createIntegrator.

private UnivariateIntegrator createIntegrator() {
    UnivariateIntegrator in = integrator;
    if (in == null) {
        // This is the integrator for the Poisson-Gamma when observed count x>=1
        // i.e. not at the boundary x=0.
        final double relativeAccuracy = 1e-4;
        final double absoluteAccuracy = 1e-16;
        int minimalIterationCount;
        switch(convolutionMode) {
            case SIMPSON_PDF:
                // This is a CustomSimpsonIntegrator that computes 1 refinement
                // on the first iteration.
                // Number of function evaluations = 2^(iteration+1) + 1
                // => 5 for 1 iterations
                // => 9 for 2 iterations
                minimalIterationCount = 1;
                in = new CustomSimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minimalIterationCount, CustomSimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
                break;
            case LEGENDRE_GAUSS_PDF:
                // Not sure how to configure this.
                // The integration points are used for each sub-interval.
                // Function evaluations = integrationpoints * intervals.
                // The intervals start at 1,2 and increase by at least 4 at each stage after that.
                // At least 1 stage is done thus 3 * integrationpoints functions evaluations
                // will be done for minimalIterationCount=1.
                minimalIterationCount = 1;
                final int maximalIterationCount = 32;
                final int integrationpoints = 8;
                in = new IterativeLegendreGaussIntegrator(integrationpoints, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
                break;
            default:
                throw new IllegalStateException();
        }
        integrator = in;
    }
    return in;
}