Examples with NewtonZeroFinder dr.math.iterations.NewtonZeroFinder used on opensource projects

Example 1 with NewtonZeroFinder

use of dr.math.iterations.NewtonZeroFinder 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 NewtonZeroFinder

use of dr.math.iterations.NewtonZeroFinder in project beast-mcmc by beast-dev.

the class PolynomialFunction method roots.

/**
 *
 * @param desiredPrecision double
 * @return double[]
 */
public double[] roots(double desiredPrecision) {
    PolynomialFunction dp = derivative();
    double start = 0;
    while (Math.abs(dp.value(start)) < desiredPrecision) start = Math.random();
    PolynomialFunction p = this;
    NewtonZeroFinder rootFinder = new NewtonZeroFinder(this, dp, start);
    rootFinder.setDesiredPrecision(desiredPrecision);
    Vector rootCollection = new Vector(degree());
    while (true) {
        rootFinder.evaluate();
        if (!rootFinder.hasConverged())
            break;
        double r = rootFinder.getResult();
        rootCollection.addElement(r);
        p = p.deflate(r);
        if (p.degree() == 0)
            break;
        rootFinder.setFunction(p);
        try {
            rootFinder.setDerivative(p.derivative());
        } catch (IllegalArgumentException e) {
        }
    }
    double[] roots = new double[rootCollection.size()];
    Enumeration e = rootCollection.elements();
    int n = 0;
    while (e.hasMoreElements()) {
        roots[n++] = (Double) e.nextElement();
    }
    return roots;
}