private double inverseProbability()

in commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java [162:255]
52 lines of code
16 McCabe index (conditional complexity)

    private double inverseProbability(final double p, final double q, boolean complement) {
        /* 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.
         *
         * In the case of the survival probability the bracket can be set using the same
         * bound given that the argument p = 1 - q, with q the survival probability.
         */

        double lowerBound = getSupportLowerBound();
        if (p == 0) {
            return lowerBound;
        }
        double upperBound = getSupportUpperBound();
        if (q == 0) {
            return upperBound;
        }

        final double mu = getMean();
        final double sig = Math.sqrt(getVariance());
        final boolean chebyshevApplies = Double.isFinite(mu) &&
                                         ArgumentUtils.isFiniteStrictlyPositive(sig);

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            lowerBound = createFiniteLowerBound(p, q, complement, upperBound, mu, sig, chebyshevApplies);
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            upperBound = createFiniteUpperBound(p, q, complement, lowerBound, mu, sig, chebyshevApplies);
        }

        // Here the bracket [lower, upper] uses finite values. If the support
        // is infinite the bracket can truncate the distribution and the target
        // probability can be outside the range of [lower, upper].
        if (upperBound == Double.MAX_VALUE) {
            if (complement) {
                if (survivalProbability(upperBound) > q) {
                    return getSupportUpperBound();
                }
            } else if (cumulativeProbability(upperBound) < p) {
                return getSupportUpperBound();
            }
        }
        if (lowerBound == -Double.MAX_VALUE) {
            if (complement) {
                if (survivalProbability(lowerBound) < q) {
                    return getSupportLowerBound();
                }
            } else if (cumulativeProbability(lowerBound) > p) {
                return getSupportLowerBound();
            }
        }

        final DoubleUnaryOperator fun = complement ?
            arg -> survivalProbability(arg) - q :
            arg -> cumulativeProbability(arg) - p;
        // Note the initial value is robust to overflow.
        // Do not use 0.5 * (lowerBound + upperBound).
        final double x = new BrentSolver(SOLVER_RELATIVE_ACCURACY,
                                         SOLVER_ABSOLUTE_ACCURACY,
                                         SOLVER_FUNCTION_VALUE_ACCURACY)
            .findRoot(fun,
                      lowerBound,
                      lowerBound + 0.5 * (upperBound - lowerBound),
                      upperBound);

        if (!isSupportConnected()) {
            return searchPlateau(complement, lowerBound, x);
        }
        return x;
    }