override fun inverseCumulativeProbability()

in plot-base/src/commonMain/kotlin/org/jetbrains/letsPlot/core/plot/base/stat/math3/AbstractRealDistribution.kt [74:177]
64 lines of code
18 McCabe index (conditional complexity)

    override fun inverseCumulativeProbability(p: Double): Double {
        /*
         * 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) {
            error("OutOfRange [0, 1] - p$p")
        }

        var lowerBound = supportLowerBound
        if (p == 0.0) {
            return lowerBound
        }

        var upperBound = supportUpperBound
        if (p == 1.0) {
            return upperBound
        }

        val mu = numericalMean
        val sig = sqrt(numericalVariance)
        val chebyshevApplies: Boolean
        chebyshevApplies = !(mu.isInfinite() || mu.isNaN() || sig.isInfinite() || sig.isNaN())

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

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

        val toSolve = object : UnivariateFunction {
            override fun value(x: Double): Double {
                return cumulativeProbability(x) - p
            }
        }

        val x = UnivariateSolverUtils.solve(
            toSolve,
            lowerBound,
            upperBound,
            solverAbsoluteAccuracy
        )

        if (!isSupportConnected) {
            /* Test for plateau. */
            val dx = solverAbsoluteAccuracy
            if (x - dx >= supportLowerBound) {
                val px = cumulativeProbability(x)
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x
                    while (upperBound - lowerBound > dx) {
                        val midPoint = 0.5 * (lowerBound + upperBound)
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint
                        } else {
                            upperBound = midPoint
                        }
                    }
                    return upperBound
                }
            }
        }
        return x
    }