org.apache.commons.math3.optim.univariate.BrentOptimizer Java Examples

/**
 * 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();
}
 

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();
}
 

/**
 * 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();
}
 

public static double argmax(final Function<Double, Double> function, final double min, final double max, final double guess,
                            final double relativeTolerance, final double absoluteTolerance, final int maxEvaluations) {
    final BrentOptimizer optimizer = new BrentOptimizer(relativeTolerance, absoluteTolerance);
    final SearchInterval interval = new SearchInterval(min, max, guess);
    return optimizer.optimize(new UnivariateObjectiveFunction(function::apply), GoalType.MAXIMIZE, interval, new MaxEval(maxEvaluations)).getPoint();
}
 

/**
 * The {@code BrentOptimizer} default stopping criterion uses the
 * tolerances to check the domain (point) values, not the function
 * values.
 * The {@code relativeTolerance} and {@code absoluteTolerance}
 * arguments are thus passed to a {@link SimpleUnivariateValueChecker
 * custom checker} that will use the function values.
 *
 * @param optimizer Optimizer on behalf of which the line search
 * be performed.
 * Its {@link MultivariateOptimizer#computeObjectiveValue(double[])
 * computeObjectiveValue} method will be called by the
 * {@link #search(double[],double[]) search} method.
 * @param relativeTolerance Search will stop when the function relative
 * difference between successive iterations is below this value.
 * @param absoluteTolerance Search will stop when the function absolute
 * difference between successive iterations is below this value.
 * @param initialBracketingRange Extent of the initial interval used to
 * find an interval that brackets the optimum.
 * If the optimized function varies a lot in the vicinity of the optimum,
 * it may be necessary to provide a value lower than the distance between
 * successive local minima.
 */
public LineSearch(MultivariateOptimizer optimizer,
                  double relativeTolerance,
                  double absoluteTolerance,
                  double initialBracketingRange) {
    mainOptimizer = optimizer;
    lineOptimizer = new BrentOptimizer(REL_TOL_UNUSED,
                                       ABS_TOL_UNUSED,
                                       new SimpleUnivariateValueChecker(relativeTolerance,
                                                                        absoluteTolerance));
    this.initialBracketingRange = initialBracketingRange;
}
 

/**
 * The {@code BrentOptimizer} default stopping criterion uses the
 * tolerances to check the domain (point) values, not the function
 * values.
 * The {@code relativeTolerance} and {@code absoluteTolerance}
 * arguments are thus passed to a {@link SimpleUnivariateValueChecker
 * custom checker} that will use the function values.
 *
 * @param optimizer Optimizer on behalf of which the line search
 * be performed.
 * Its {@link MultivariateOptimizer#computeObjectiveValue(double[])
 * computeObjectiveValue} method will be called by the
 * {@link #search(double[],double[]) search} method.
 * @param relativeTolerance Search will stop when the function relative
 * difference between successive iterations is below this value.
 * @param absoluteTolerance Search will stop when the function absolute
 * difference between successive iterations is below this value.
 * @param initialBracketingRange Extent of the initial interval used to
 * find an interval that brackets the optimum.
 * If the optimized function varies a lot in the vicinity of the optimum,
 * it may be necessary to provide a value lower than the distance between
 * successive local minima.
 */
public LineSearch(MultivariateOptimizer optimizer,
                  double relativeTolerance,
                  double absoluteTolerance,
                  double initialBracketingRange) {
    mainOptimizer = optimizer;
    lineOptimizer = new BrentOptimizer(REL_TOL_UNUSED,
                                       ABS_TOL_UNUSED,
                                       new SimpleUnivariateValueChecker(relativeTolerance,
                                                                        absoluteTolerance));
    this.initialBracketingRange = initialBracketingRange;
}
 

/**
 * Maximizes a univariate function using a grid search followed by Brent's algorithm.
 *
 * @param fn the likelihood function to minimize
 * @param gridStart the lower bound for the grid search
 * @param gridEnd the upper bound for the grid search
 * @param gridStep step size for the grid search
 * @param relErr relative error tolerance for Brent's algorithm
 * @param absErr absolute error tolerance for Brent's algorithm
 * @param maxIter maximum # of iterations to perform in Brent's algorithm
 * @param maxEval maximum # of Likelihood function evaluations in Brent's algorithm
 *
 * @return the value of the parameter that maximizes the function
 */
public static double maximize(UnivariateFunction fn, double gridStart, double gridEnd,
    double gridStep, double relErr, double absErr, int maxIter, int maxEval) {
  Interval interval = gridSearch(fn, gridStart, gridEnd, gridStep);
  BrentOptimizer bo = new BrentOptimizer(relErr, absErr);
  UnivariatePointValuePair max = bo.optimize(
      new MaxIter(maxIter),
      new MaxEval(maxEval),
      new SearchInterval(interval.getInf(), interval.getSup()),
      new UnivariateObjectiveFunction(fn),
      GoalType.MAXIMIZE);
  return max.getPoint();
}