PoissonDistribution.java

/*
 *      Copyright (C) 2015  higherfrequencytrading.com
 *      Copyright 2001-2015 The Apache Software Foundation
 *      Copyright 2010-2012 CS Systèmes d'Information
 *
 *      This program is free software: you can redistribute it and/or modify
 *      it under the terms of the GNU Lesser General Public License as published by
 *      the Free Software Foundation, either version 3 of the License.
 *
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU Lesser General Public License for more details.
 *
 *      You should have received a copy of the GNU Lesser General Public License
 *      along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package net.openhft.chronicle.hash.impl.util.math;

public class PoissonDistribution {

    private static final double EPSILON = 1e-12;
    /**
     * Poisson distribution is used to estimate segment fillings. Segments are not bound with
     * Integer.MAX_VALUE, but it's not clear if algorithms from Commons Math could work with such
     * big values as Long.MAX_VALUES
     */
    private static final long UPPER_BOUND = 1L << 36;
    private static final int MAX_ITERATIONS = 10000000;

    public static double cumulativeProbability(double mean, long x) {
        if (x < 0) {
            return 0;
        }
        if (x >= UPPER_BOUND) {
            return 1;
        }
        return Gamma.regularizedGammaQ((double) x + 1, mean, EPSILON, MAX_ITERATIONS);
    }

    public static long inverseCumulativeProbability(double mean, double p) {
        checkProbability(p);

        long lower = 0;
        if (p == 0.0) {
            return lower;
        }
        lower -= 1; // this ensures cumulativeProbability(lower) < p, which
        // is important for the solving step

        long upper = UPPER_BOUND;
        if (p == 1.0) {
            return upper;
        }

        // use the one-sided Chebyshev inequality to narrow the bracket
        // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
        final double mu = mean;
        // in Poisson distribution, variance == mean
        double variance = mean;
        final double sigma = Math.sqrt(variance);
        final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
                Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0);
        if (chebyshevApplies) {
            double k = Math.sqrt((1.0 - p) / p);
            double tmp = mu - k * sigma;
            if (tmp > lower) {
                lower = ((int) Math.ceil(tmp)) - 1;
            }
            k = 1.0 / k;
            tmp = mu + k * sigma;
            if (tmp < upper) {
                upper = ((int) Math.ceil(tmp)) - 1;
            }
        }

        return solveInverseCumulativeProbability(mean, p, lower, upper);
    }

    private static void checkProbability(double p) {
        if (p < 0.0 || p > 1.0) {
            throw new IllegalArgumentException("probability should be in [0.0, 1.0] bounds, " + p +
                    " given");
        }
    }

    /**
     * This is a utility function used by {@link
     * #inverseCumulativeProbability(double, double)}. It assumes {@code 0 < p < 1} and
     * that the inverse cumulative probability lies in the bracket {@code
     * (lower, upper]}. The implementation does simple bisection to find the
     * smallest {@code p}-quantile <code>inf{x in Z | P(X<=x) >= p}</code>.
     *
     * @param p the cumulative probability
     * @param lower a value satisfying {@code cumulativeProbability(lower) < p}
     * @param upper a value satisfying {@code p <= cumulativeProbability(upper)}
     * @return the smallest {@code p}-quantile of this distribution
     */
    private static long solveInverseCumulativeProbability(double mean, final double p,
                                                         long lower, long upper) {

        while (lower + 1 < upper) {
            long xm = (lower + upper) / 2;
            if (xm < lower || xm > upper) {
                /*
                 * Overflow.
                 * There will never be an overflow in both calculation methods
                 * for xm at the same time
                 */
                xm = lower + (upper - lower) / 2;
            }

            double pm = checkedCumulativeProbability(mean, xm);
            if (pm >= p) {
                upper = xm;
            } else {
                lower = xm;
            }
        }
        return upper;
    }

    public static double meanByCumulativeProbabilityAndValue(double p, long x, double precision) {
        checkProbability(p);
        assert x > 0 && x < UPPER_BOUND;

        double lower = 0;
        double upper = UPPER_BOUND;
        while (lower + precision < upper) {
            double m = (lower + upper) / 2;
            double pm = checkedCumulativeProbability(m, x);
            if (pm < p) {
                upper = m;
            } else {
                lower = m;
            }
        }
        return lower;
    }

    /**
     * Computes the cumulative probability function and checks for {@code NaN}
     * values returned. Throws {@code MathInternalError} if the value is
     * {@code NaN}. Rethrows any exception encountered evaluating the cumulative
     * probability function. Throws {@code MathInternalError} if the cumulative
     * probability function returns {@code NaN}.
     *
     * @param argument input value
     * @return the cumulative probability
     * @throws AssertionError if the cumulative probability is {@code NaN}
     */
    private static double checkedCumulativeProbability(double mean, long argument) {
        double result = cumulativeProbability(mean, argument);
        if (Double.isNaN(result)) {
            throw new AssertionError("Discrete cumulative probability function returned NaN " +
                    "for argument " + argument);
        }
        return result;
    }
}