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