Search in sources :

Example 6 with BrentOptimizer

use of org.apache.commons.math3.optim.univariate.BrentOptimizer in project gatk by broadinstitute.

the class PosteriorSummaryUtils method calculatePosteriorMode.

/**
     * Given a list of posterior samples, returns an estimate of the posterior mode (using
     * mllib kernel density estimation in {@link KernelDensity} and {@link BrentOptimizer}).
     * Note that estimate may be poor if number of samples is small (resulting in poor kernel density estimation),
     * or if posterior is not unimodal (or is sufficiently pathological otherwise). If the samples contain
     * {@link Double#NaN}, {@link Double#NaN} will be returned.
     * @param samples   posterior samples, cannot be {@code null} and number of samples must be greater than 0
     * @param ctx       {@link JavaSparkContext} used by {@link KernelDensity} for mllib kernel density estimation
     */
public static double calculatePosteriorMode(final List<Double> samples, final JavaSparkContext ctx) {
    Utils.nonNull(samples);
    Utils.validateArg(samples.size() > 0, "Number of samples must be greater than zero.");
    //calculate sample min, max, mean, and standard deviation
    final double sampleMin = Collections.min(samples);
    final double sampleMax = Collections.max(samples);
    final double sampleMean = new Mean().evaluate(Doubles.toArray(samples));
    final double sampleStandardDeviation = new StandardDeviation().evaluate(Doubles.toArray(samples));
    //if samples are all the same or contain NaN, can simply return mean
    if (sampleStandardDeviation == 0. || Double.isNaN(sampleMean)) {
        return sampleMean;
    }
    //use Silverman's rule to set bandwidth for kernel density estimation from sample standard deviation
    //see https://en.wikipedia.org/wiki/Kernel_density_estimation#Practical_estimation_of_the_bandwidth
    final double bandwidth = SILVERMANS_RULE_CONSTANT * sampleStandardDeviation * Math.pow(samples.size(), SILVERMANS_RULE_EXPONENT);
    //use kernel density estimation to approximate posterior from samples
    final KernelDensity pdf = new KernelDensity().setSample(ctx.parallelize(samples, 1)).setBandwidth(bandwidth);
    //use Brent optimization to find mode (i.e., maximum) of kernel-density-estimated posterior
    final BrentOptimizer optimizer = new BrentOptimizer(RELATIVE_TOLERANCE, RELATIVE_TOLERANCE * (sampleMax - sampleMin));
    final UnivariateObjectiveFunction objective = new UnivariateObjectiveFunction(f -> pdf.estimate(new double[] { f })[0]);
    //search for mode within sample range, start near sample mean
    final SearchInterval searchInterval = new SearchInterval(sampleMin, sampleMax, sampleMean);
    return optimizer.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
Also used : Mean(org.apache.commons.math3.stat.descriptive.moment.Mean) SearchInterval(org.apache.commons.math3.optim.univariate.SearchInterval) UnivariateObjectiveFunction(org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction) BrentOptimizer(org.apache.commons.math3.optim.univariate.BrentOptimizer) KernelDensity(org.apache.spark.mllib.stat.KernelDensity) StandardDeviation(org.apache.commons.math3.stat.descriptive.moment.StandardDeviation)

Example 7 with BrentOptimizer

use of org.apache.commons.math3.optim.univariate.BrentOptimizer in project gatk by broadinstitute.

the class CNLOHCaller method optimizeIt.

private double optimizeIt(final Function<Double, Double> objectiveFxn, final SearchInterval searchInterval) {
    final MaxEval BRENT_MAX_EVAL = new MaxEval(1000);
    final double RELATIVE_TOLERANCE = 0.001;
    final double ABSOLUTE_TOLERANCE = 0.001;
    final BrentOptimizer OPTIMIZER = new BrentOptimizer(RELATIVE_TOLERANCE, ABSOLUTE_TOLERANCE);
    final UnivariateObjectiveFunction objective = new UnivariateObjectiveFunction(x -> objectiveFxn.apply(x));
    return OPTIMIZER.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
Also used : MaxEval(org.apache.commons.math3.optim.MaxEval) UnivariateObjectiveFunction(org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction) BrentOptimizer(org.apache.commons.math3.optim.univariate.BrentOptimizer)

Aggregations

BrentOptimizer (org.apache.commons.math3.optim.univariate.BrentOptimizer)7 UnivariateObjectiveFunction (org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction)6 MaxEval (org.apache.commons.math3.optim.MaxEval)4 SearchInterval (org.apache.commons.math3.optim.univariate.SearchInterval)4 Mean (org.apache.commons.math3.stat.descriptive.moment.Mean)2 StandardDeviation (org.apache.commons.math3.stat.descriptive.moment.StandardDeviation)2 KernelDensity (org.apache.spark.mllib.stat.KernelDensity)2 MedianWindow (gdsc.core.utils.MedianWindow)1 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)1 FRC (gdsc.smlm.ij.frc.FRC)1 FRCCurve (gdsc.smlm.ij.frc.FRC.FRCCurve)1 FRCCurveResult (gdsc.smlm.ij.frc.FRC.FRCCurveResult)1 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)1 Plot2 (ij.gui.Plot2)1 LoessInterpolator (org.apache.commons.math3.analysis.interpolation.LoessInterpolator)1 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)1 SimpleCurveFitter (org.apache.commons.math3.fitting.SimpleCurveFitter)1 WeightedObservedPoint (org.apache.commons.math3.fitting.WeightedObservedPoint)1 WeightedObservedPoints (org.apache.commons.math3.fitting.WeightedObservedPoints)1 PointValuePair (org.apache.commons.math3.optim.PointValuePair)1