Examples with OneVariableFunction dr.math.interfaces.OneVariableFunction used on opensource projects

Example 1 with OneVariableFunction

use of dr.math.interfaces.OneVariableFunction in project beast-mcmc by beast-dev.

the class InverseGaussianDistribution method quantile.

/**
     * quantiles (=inverse cumulative density function)
     * <p/>
     *
     * Same implementation as SuppleDists in R.
     *
     * Using Whitmore and Yalovsky for an initial guess. Works well for
	 * large t=lambda/mu > 2 perhaps
	 * Whitmore, G.A. and Yalovsky, M. (1978). A normalizing logarithmic
	 * transformation for inverse Gaussian random variables,
	 * Technometrics 20-2, 207-208
	 * For small t, with x<0.5 mu, use gamma approx to 1/x -- alpha=1/2 and beta =2/lambda and 1-p
	 * When x>0.5mu, approx x with gamma for p, and exponentiate -- don't know why this works.
     *
     * There are cases  which even this method produces inaccurate results (e.g. when shape = 351). Therefore,
     * we have a catch that will determine whether or not the result is accurate enough and if not,
     * then a zerofinder will be used to find a more accurate approximation. 
     *
     * @param z     argument
     * @param m     mean
     * @param shape shape parameter
     * @return icdf at z
     */
public static double quantile(double z, double m, double shape) {
    if (z < 0.01 || z > 0.99) {
        System.err.print("Quantile is " + z);
        throw new RuntimeException("Quantile is too low/high to calculate (numerical estimation for extreme values is incomplete)");
    }
    /* Approximation method used by Mudholkar GS, Natarajan R (1999)
         * Approximations for the inverse gaussian probabilities and percentiles.
         * Communications in Statistics - Simulation and Computation 28: 1051 - 1071.
         */
    double initialGuess;
    if (shape / m > 2.0) {
        initialGuess = (NormalDistribution.quantile(z, 0.0, 1.0) - 0.5 * Math.sqrt(m / shape)) / Math.sqrt(shape / m);
        initialGuess = m * Math.exp(initialGuess);
    } else {
        initialGuess = shape / (GammaDistribution.quantile(1.0 - z, 0.5, 1.0) * 2.0);
        if (initialGuess > m / 2.0) {
            // too large for the gamma approx
            // this seems to work for the upper tail ???
            initialGuess = m * Math.exp(GammaDistribution.quantile(z, 0.5, 1.0) * 0.1);
        }
    }
    //        double phi = shape / m;
    //        if(phi>50.0) {
    // Use Normal Distribution
    //            initialGuess = (NormalDistribution.quantile(z, m,Math.sqrt(m*m*m/shape)));//-0.5*Math.sqrt(m/shape))/Math.sqrt(m*m*m/shape);
    //        }
    final InverseGaussianDistribution f = new InverseGaussianDistribution(m, shape);
    final double y = z;
    NewtonZeroFinder zeroFinder = new NewtonZeroFinder(new OneVariableFunction() {

        public double value(double x) {
            return f.cdf(x) - y;
        }
    }, initialGuess);
    zeroFinder.evaluate();
    if (Double.isNaN(zeroFinder.getResult()) || zeroFinder.getPrecision() > 0.000005) {
        zeroFinder = new NewtonZeroFinder(new OneVariableFunction() {

            public double value(double x) {
                return f.cdf(x) - y;
            }
        }, initialGuess);
        zeroFinder.initializeIterations();
        int i;
        double previousPrecision = 0.0, previousResult = Double.NaN;
        double max = 10000.0, min = 0.00001;
        for (i = 0; i < 50; i++) {
            zeroFinder.evaluateIteration();
            double precision = f.cdf(zeroFinder.getResult()) - z;
            if ((previousPrecision > 0 && precision < 0) || (previousPrecision < 0 && precision > 0)) {
                max = Math.max(previousResult, zeroFinder.getResult());
                min = Math.min(previousResult, zeroFinder.getResult());
                max = Math.min(10000.0, max);
                break;
            }
            previousPrecision = precision;
            previousResult = zeroFinder.getResult();
        }
        return calculateZeroFinderApproximation(z, m, shape, min, max, initialGuess);
    }
    return zeroFinder.getResult();
}

Example 2 with OneVariableFunction

use of dr.math.interfaces.OneVariableFunction in project beast-mcmc by beast-dev.

the class InverseGaussianDistribution method calculateZeroFinderApproximation.

/** Finds the approximation of the inverse Gaussian quantile using a zero finder
     * until it converges
     *
     * @param z            quantile
     * @param m            mean
     * @param shape        shape
     * @param min          min search value
     * @param max          max search value
     * @param initialGuess first guess of the quantile
     * @return estimated x value at quantile z
     */
private static double calculateZeroFinderApproximation(double z, double m, double shape, double min, double max, double initialGuess) {
    final InverseGaussianDistribution f = new InverseGaussianDistribution(m, shape);
    final double y = z;
    BisectionZeroFinder bisectionZeroFinder = new BisectionZeroFinder(new OneVariableFunction() {

        public double value(double x) {
            return f.cdf(x) - y;
        }
    }, min, max);
    bisectionZeroFinder.setInitialValue(initialGuess);
    bisectionZeroFinder.initializeIterations();
    double bestValue = Double.NaN;
    /* I found that the converged value is not necesssarily the best */
    double bestPrecision = 10;
    double precision = 10;
    double previousPrecision = 10;
    int count = 0;
    while (precision > 0.001 && count < 10) {
        bisectionZeroFinder.evaluateIteration();
        precision = Math.abs(f.cdf(bisectionZeroFinder.getResult()) - z);
        if (precision < bestPrecision) {
            bestPrecision = precision;
            bestValue = bisectionZeroFinder.getResult();
        } else if (previousPrecision == precision) {
            count++;
        }
        previousPrecision = precision;
    }
    bisectionZeroFinder.finalizeIterations();
    //return bisectionZeroFinder.getResult();
    return bestValue;
}

Example 3 with OneVariableFunction

use of dr.math.interfaces.OneVariableFunction in project beast-mcmc by beast-dev.

the class InverseGaussianDistributionTest method testCDFAndQuantile2.

public void testCDFAndQuantile2() {
    final InverseGaussianDistribution f = new InverseGaussianDistribution(1, 1);
    for (double i = 0.01; i < 0.95; i += 0.01) {
        final double y = i;
        BisectionZeroFinder zeroFinder = new BisectionZeroFinder(new OneVariableFunction() {

            public double value(double x) {
                return f.cdf(x) - y;
            }
        }, 0.01, 100);
        zeroFinder.setMaximumIterations(100);
        zeroFinder.evaluate();
        assertEquals(f.quantile(i), zeroFinder.getResult(), 1e-3);
    }
}

Example 4 with OneVariableFunction

use of dr.math.interfaces.OneVariableFunction in project beast-mcmc by beast-dev.

the class LogNormalDistributionTest method testCDFAndQuantile2.

public void testCDFAndQuantile2() {
    final LogNormalDistribution f = new LogNormalDistribution(1, 1);
    for (double i = 0.01; i < 0.95; i += 0.01) {
        final double y = i;
        BisectionZeroFinder zeroFinder = new BisectionZeroFinder(new OneVariableFunction() {

            public double value(double x) {
                return f.cdf(x) - y;
            }
        }, 0.01, 100);
        zeroFinder.evaluate();
        assertEquals(f.quantile(i), zeroFinder.getResult(), 1e-6);
    }
}

Example 5 with OneVariableFunction

use of dr.math.interfaces.OneVariableFunction in project beast-mcmc by beast-dev.

the class InverseGaussianDistribution method calculateZeroFinderApproximation.

/** Finds the approximation of the inverse Gaussian quantile using a zero finder
     * given a set number of maximum iterations
     * UNUSED METHOD
     *
     * @param z             quantile
     * @param m             mean
     * @param shape         shape
     * @param numIterations number of iterations used to approximate value
     * @param min           min search value
     * @param max           max search value
     * @return estimated x value at quantile z
     */
private static double calculateZeroFinderApproximation(double z, double m, double shape, int numIterations, double min, double max) {
    final InverseGaussianDistribution f = new InverseGaussianDistribution(m, shape);
    final double y = z;
    BisectionZeroFinder bisectionZeroFinder = new BisectionZeroFinder(new OneVariableFunction() {

        public double value(double x) {
            return f.cdf(x) - y;
        }
    }, min, max);
    //}, 0.0001, 100000);
    bisectionZeroFinder.setMaximumIterations(numIterations);
    bisectionZeroFinder.evaluate();
    return bisectionZeroFinder.getResult();
}