Example usage for org.apache.commons.math3.analysis.solvers UnivariateSolverUtils solve

Introduction

In this page you can find the example usage for org.apache.commons.math3.analysis.solvers UnivariateSolverUtils solve.

Prototype

public static double solve(UnivariateFunction function, double x0, double x1, double absoluteAccuracy)
        throws NullArgumentException, NoBracketingException 

Source Link

Document

Convenience method to find a zero of a univariate real function.

Usage

From source file:org.nd4j.linalg.api.rng.distribution.BaseDistribution.java

/**
 * {@inheritDoc}//  ww w.  j a  v  a 2  s  .  c om
 * <p/>
 * The default implementation returns
 * <ul>
 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
 * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
 * </ul>
 */
@Override
public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
    /*
     * IMPLEMENTATION NOTES
     * --------------------
     * Where applicable, use is made of the one-sided Chebyshev inequality
     * to bracket the root. This inequality states that
     * P(X - mu >= k * sig) <= 1 / (1 + k^2),
     * mu: mean, sig: standard deviation. Equivalently
     * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
     * F(mu + k * sig) >= k^2 / (1 + k^2).
     *
     * For k = sqrt(p / (1 - p)), we find
     * F(mu + k * sig) >= p,
     * and (mu + k * sig) is an upper-bound for the root.
     *
     * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
     * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
     * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
     * P(X <= mu - k * sig) <= 1 / (1 + k^2),
     * F(mu - k * sig) <= 1 / (1 + k^2).
     *
     * For k = sqrt((1 - p) / p), we find
     * F(mu - k * sig) <= p,
     * and (mu - k * sig) is a lower-bound for the root.
     *
     * In cases where the Chebyshev inequality does not apply, geometric
     * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
     * the root.
     */
    if (p < 0.0 || p > 1.0) {
        throw new OutOfRangeException(p, 0, 1);
    }

    double lowerBound = getSupportLowerBound();
    if (p == 0.0) {
        return lowerBound;
    }

    double upperBound = getSupportUpperBound();
    if (p == 1.0) {
        return upperBound;
    }

    final double mu = getNumericalMean();
    final double sig = FastMath.sqrt(getNumericalVariance());
    final boolean chebyshevApplies;
    chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig)
            || Double.isNaN(sig));

    if (lowerBound == Double.NEGATIVE_INFINITY) {
        if (chebyshevApplies) {
            lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
        } else {
            lowerBound = -1.0;
            while (cumulativeProbability(lowerBound) >= p) {
                lowerBound *= 2.0;
            }
        }
    }

    if (upperBound == Double.POSITIVE_INFINITY) {
        if (chebyshevApplies) {
            upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
        } else {
            upperBound = 1.0;
            while (cumulativeProbability(upperBound) < p) {
                upperBound *= 2.0;
            }
        }
    }

    final UnivariateFunction toSolve = new UnivariateFunction() {

        public double value(final double x) {
            return cumulativeProbability(x) - p;
        }
    };

    double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy());

    if (!isSupportConnected()) {
        /* Test for plateau. */
        final double dx = getSolverAbsoluteAccuracy();
        if (x - dx >= getSupportLowerBound()) {
            double px = cumulativeProbability(x);
            if (cumulativeProbability(x - dx) == px) {
                upperBound = x;
                while (upperBound - lowerBound > dx) {
                    final double midPoint = 0.5 * (lowerBound + upperBound);
                    if (cumulativeProbability(midPoint) < px) {
                        lowerBound = midPoint;
                    } else {
                        upperBound = midPoint;
                    }
                }
                return upperBound;
            }
        }
    }
    return x;
}

---