/* * 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; } }